The use of the pilot points method on groundwater modelling for a degraded aquifer with limited field data: the case of Lake Karla aquifer Open Access

Groundwater resources are vital in providing water for domestic, agricultural, and industrial use globally. Over two billion people rely on groundwater as their primary water source, while half or more of the irrigation water for agriculture is supplied from underground sources (Pisinaras et al. 2013; Famiglietti 2014). However, despite its critical importance to global water supplies, inadequate management practices are applied for monitoring groundwater resources. Furthermore, it is crucial to estimate the pressure on surface and groundwater systems due to groundwater depletion. Climate change and related changes on the components of the hydrological cycle complicate the task of sustaining groundwater supplies for the future. Hence, about 20% of the world's aquifers systems are estimated to be under over exploitation status and with increasing trends due to the climate change phenomenon and population growth (Upton et al. 2020).

Understanding why groundwater depletion has occurred may be the first step towards sustainable groundwater resources management. Numerical simulation models are effective and valuable tools for computing water budgets of regional groundwater systems and the planning and policy making for groundwater resources management (Konikow & Kendy 2005). The analysis of groundwater system through simulation is the main process before a researcher/modeller deals with the management of groundwater resources (Doherty & Hunt 2010). The simulation problems in groundwater hydrology are part of the complex environmental issues that are difficult to solve. These problems in groundwater modelling become extremely difficult when in addition to the lack of field data problems also meet a lack of quality data (Singh 2014).

One of the most important issues for water resources management is developing strategies for groundwater modelling that are adaptable to data scarcity and field measurements and so begins the difficulty of the inverse problem. According to Zhou and her associates (2014) many approaches have been proposed to resolve the inverse problem in hydrogeological modelling. Inverse methods are categorized as direct and indirect ones. The inverse problem in hydrogeological fields is difficult to solve because it is ill posed and the calibration parameters are more than the observations (Anderson et al. 2015). Common indirect methods used to solve the inverse problem are the geostatistical approach, the maximum likelihood method where zonation is employed, and pilot points. The most advanced methods where the stochastic approach is introduced are the self-calibrated method, the Markov chain Monte Carlo method, the gradual deformation method, the probability perturbation method, and the ensemble Kalman filter. The geostatistical approach has several advantages, however it is not appropriate for use in the inverse model where the hydraulic conductivity variable poses spatial distributions with significant heterogeneity, forcing the transition from linearized to non-linearized approaches (Zhou et al. 2014).

Zonation assumes an abrupt spatial change in parameter values, whereas the pilot points method produces smoothly distributed parameters and flexibility compared to the zonation approach (Anderson et al. 2015; Rajabi et al. 2018). These strategies are particularly important in arid and semi-arid areas where access to data is poor and data collection is difficult, such as the Lake Karla aquifer, Greece. This study analyzes the groundwater resources status of the Lake Karla aquifer focusing on a precise simulation–validation procedure of the groundwater model. A surface hydrological simulation tool was applied to calculate the groundwater recharge while the MODFLOW with a Newton formulation was employed to simulate groundwater flow. The groundwater model was calibrated for spatial distributed hydraulic conductivity based on the pilot points method, since there are limited data for this parameter.

This paper presents a new framework for calibrating and validating the groundwater flow simulation applying two automated tools with the use of pilot points. The first relates to calibration post-process analysis and is called PLPROC and the second one deals with the post-process analysis of the results and is called mod2obs (Doherty 2016, 2018). Both tools are part of the PEST: Model-Independent Parameter Estimation and Uncertainty Analysis, an open-source, public-domain software suite (Doherty 2020a, 2020b). The PLPROC tool has been applied in the Precipice Sandstone aquifer system of the Surat Basin, Great Artesian Basin, in Australia, by Hayes et al. (2020). The overall purpose of this paper is to provide a useful framework in surface and groundwater modelling dealing with scarce hydrogeological information in rural areas.

The study area of about 500 km² is located at Eastern Thessaly, Greece and occupies the plain zone of Lake Karla watershed (Figure 1). The mean annual values of precipitation, temperature and potential evapotranspiration of the region are 560 mm, 14 °C and 775 mm respectively (Mylopoulos & Sidiropoulos 2016). The grainy, phreatic aquifer consists of alluvial deposits (Quaternary formations). The thickness varies in places reaching its maximum value of 200 m at the centre of study area. The presence of Mavrovounio Mountain at the east, where low permeable bedrocks prevail at the borderline with the lowland creates a no flow boundary. At the west, the permeable formations of Chalkodonio Mountain (south-west) and of the plain (west north) allow a hydraulic connection of adjacent aquifers with Lake Karla aquifer. The main crops are wheat, cotton, alfalfa, maize, vegetables, and trees. Lake Karla's aquifer is a characteristic example of an agricultural area in Greek territory, where surface water systems are limited and have relatively small capacity, in contrast to the great irrigation water demands imposed by the intensive agricultural activity of the area. As a result, the largest part of the region's water needs for irrigation is covered by groundwater with the use of pumping wells, most of them being private wells (Mylopoulos & Sidiropoulos 2016).

These agricultural practices and the absence of collective irrigation networks result in pumping non-renewable groundwater resources leading to groundwater depletion. Agriculture is the primary income source for most of the local population and for decades agricultural production in the region was based on subsidies from the Common Agriculture Policy (CAP) of the European Union. A typical example is the cultivation of cotton, a high-water-demanding crop which, for years, is based on the provided subsidies without considering sustainable management of groundwater resources (Alamanos et al. 2020). Previous studies have shown that the aquifer is in over exploitation condition and effective restoration practices are needed for sustainable water resources management (Sidiropoulos et al. 2015, 2019).

The methodology is described in Figure 2 with the use of a flowchart. The first step is running the UTHBAL surface hydrology tool (Loukas et al. 2007). One of the outcomes of UTHBAL is groundwater recharge that is used as input parameter in the MODFLOW_NWT code (Niswonger et al. 2011). The preprocessor PLPROC was used during the parameterization while the optimization process follows the PEST methodology. Because the inverse problem is an ill-posed problem the regularization mode is used to prevent overfitting. Finally, the objective function results are examined. If the objective function values show a clear decrease, then post-process analysis follows (mod2obs). When there is no observed decrease, we can conclude that the groundwater hydrology model is unstable and needs to be re-evaluated.

The most common statistical criteria for the performance of a simulation model/tool are the Mean Absolute Error (MAE), Root Mean Square Error (RMSE) and for transient conditions is the Nash–Sutcliffe Efficiency (E_f). An acceptable criterion for the successful completion of calibration is a RMSE deviance of 10% of the water levels' total range. MAE is more robust in estimating the effects of outlier values of the spatial distribution than the RMSE, and therefore it is usually smaller than RMSE. However, if the calibration procedure should be continued, depends on score of the Nash–Sutcliffe Efficiency (E_f) (Anderson et al. 2015).

The surface hydrological model (UTHBAL), is a surface hydrological simulation tool that gives a representation of the rainfall–runoff interactions of a catchment area. UTHBAL uses the components of the hydrological cycle, monthly time series of precipitation (P), mean temperature (T) and potential evapotranspiration (Ep) as inputs. The surface hydrological simulation tool disaggregates the total precipitation into rainfall and snowfall. UTHBAL distinguishes the total watershed runoff (Qsim), between the surface runoff (SR), the interflow (MR), and the baseflow (Qg) utilizing a soil moisture mechanism (Nmoist) and finally estimates the groundwater recharge to the subsurface systems (Loukas et al. 2007). In this study, UTHBAL has been applied in semi-distributed mode because the aquifer system covers the lower zone of the Karla watershed (e.g. Sidiropoulos et al. 2019). The groundwater recharge rate varied for each stress period.

The annual crop data from 2001 to 2016 come from the Integrated Administration and Control System (IACS), which is maintained by the Greek Payment Authority of Common Agricultural Policy (CAP) Aid Schemes (OPEKEPE, online). The crop pattern changes annually and the irrigation demands were estimated for every month from 2001 to 2016 with Near Irrigation Requirement (NIR). Hence, the pumping rates varied for the dry/pumping period.

Table 1 presents the estimation of irrigation water demands covered by groundwater extraction for the period 2001–2016. Table 2 shows the pumping rates for each zone and the respective irrigation area. It is worth noting that although the reservoir exists at present, it does not participate in the current study. The former natural Lake Karla had been drained in the 1960s, causing serious environmental problems in the area. A reservoir for the restoration of Lake Karla has been recently constructed to reverse the bad environmental status. The Lake Karla reservoir project construction has been completed in 2018 (Sidiropoulos et al. 2021). Thus, its contribution to the groundwater system is not considered during the simulation.

Groundwater flow simulation is achieved with the use of the MODFLOW_NWT code through the commercial GMS graphic user interface (Aquaveo 2019). MODFLOW_NWT code is a Newton formulation for MODFLOW-2005 (Niswonger et al. 2011). A no flow boundary was used in the east part to simulate the low permeable conditions of the aquifer. While in the west, general-head boundary was used to simulate the hydraulic connection of the aquifer with the adjacent ones. Conductance was estimated by the comparison of the observation wells on either side of the border ranging from 15 to 40 m²/d. The aquifer is represented as a single layer and the resulting grid consists of 12,500 active cells with 200 m side each. Grid refinement is observed at the wells' locations. This study examines the groundwater resources status extending the time series to November 2018 and applying the MODFLOW_NWT code instead of MODFLOW 2000. A transient groundwater flow model was developed and simulated from 1987 to 2018 with monthly changing irrigation needs from 2001 to 2016 also taking under consideration the annually different crop pattern application. Additionally, the groundwater flow was simulated with varying timesteps in accordance with the exact dates of water level observations. MODFLOW_NWT was chosen because the implemented method does not deactivate dry cells, but instead uses an upstream weighting approach to limit flow from each dry cell to an adjacent, partially saturated cell. It also uses additional smoothing for conductance and storage functions which allows the Newton numerical method application.

The development and application of a simulation model relies heavily on knowledge of the dominant hydrogeology conditions in the study area, the initial piezometry, the boundary conditions, the grid geometry, the trends in the groundwater levels and the assessment of hydraulic properties (Sundell et al. 2019). Once a groundwater model is developed and well calibrated using a technically sound hydrogeological structure for pre-process and post-process, its outcome estimates the groundwater budget on a regional scale and can be used to analyze future water management scenarios (Konikow & Kendy 2005).

The calibration process attempts to bring the model response as close to the observed ones as possible. In the first stage, a trial-and-error method is employed whereas, in the next stage, the pilot points approach was chosen. The pilot point approach is followed due to the lack of field measurements of hydraulic conductivity. Figure 1 indicates that there are only 32 measurements of hydraulic conductivity in the study area which is not enough to correctly express the heterogeneity of such a complex aquifer system. Pilot points are used to express the heterogeneity of hydraulic variables in groundwater systems naturally without the problem of numerical instability in fields with limited measurements data (Doherty & Hunt 2010; Sundell et al. 2019).

As a general approach in groundwater modelling pilot points are located between the observation boreholes and the pumping wells in the direction of groundwater flow (Doherty & Hunt 2010; Baalousha et al. 2019). In this paper, spatial variations of hydraulic conductivity parameter were estimated using pilot points. The pilot points placement can be seen in detail in Figure 3(a) whereas Figure 3(b) is an illustration of groundwater flow. In our simulations 254 pilot points were not randomly chosen but manually set between the observation boreholes and pumping wells. However, it is worth noting that increasing the number of pilot points in the model lengthens the simulation duration. Specific yield value was set as 0.001 corresponding to the average value between fine sand and clay, according to (Johnson 1967). Hydraulic head values were calculated by following the inverse calibration process based on hydraulic conductivity variations.

In the presented experiments the ill-posed inverse problem was solved with the application of Tikhonov regularization method. The study uses 32 observation wells monitoring the groundwater status during 1987 (Initial Conditions), while the calibration and the validation process was conducted using datasets from 10 observation wells that monitored the groundwater status during 2015 and 2018, respectively.

Post-processing data were analyzed using mod2obs (Doherty 2018), a PEST individual module that compares the simulated hydraulic head with the calibration target. Mod2obs converts the raw model outputs to the model generated observations. Mod2obs works like most programs provided by the PEST suite. The user needs to have prepared the dimension files of the model grid and the coordinates file of the observation points. Then mod2obs outputs a bore sample file of the model generated heads over time, interpolated to the sites and times in the borehole sample file (Doherty 2018).

Figures 5–7 demonstrate two types of results. The calibration and the validation using the pilot points methodology that show satisfactorily results whereas a high drawdown in hydraulic heads is observed. The above findings strongly suggest and verify the presence of groundwater depletion. The phenomenon can also be observed when examining individually boreholes AD11, SR30 and SR31 which are in the central part of the study area individually including SR63a which is located in the southwest part of the study area (Figure 1). The exact drop figures for each borehole can be seen in the Table 4.

Figure 6 shows the results achieved through calibration and validation process where the observed heads converge with the simulated ones before the pumping period and at the end of the pumping period. Large groundwater fluctuations in hydraulic heads are observed between the wet and dry periods. Statistical values for each borehole are shown in the Table 5.

Figure 7 shows the groundwater balance of the area for the period 1987–2018. The groundwater recharge variable is the sum of deep infiltration variable simulated by the UTHBAL model plus the returned irrigation flow was approximately accounted as 10% of the irrigation requirements. The groundwater balance is negative throughout the years, indicating the intensive groundwater abstraction over the region's capacity to replenish the outflows with adequate recharge to minimize the deficit water balance.

The simulation of the subsurface system of Lake Karla watershed has been presented in the past by Sidiropoulos et al. (2015, 2019). Sidiropoulos et al. (2015, 2019) used the MODFLOW 2000 code and 15 observation points for hydraulic conductivity parameter applying the geostatistical approach and the model ran from 1987 to 2009. In the current study the model was expanded forward to simulate the groundwater flow from 1987 to 2018 using an updated code, the MODFLOW_NWT mode, which resolves the problem of dry cells, uses 254 pilot points, and updates the groundwater pumping data and observation during the simulation, and utilizes the automated calibration–validation tools PLPROC and mod2obs. The methodology of pilot points deals with the data scarcity, the MODFLOW_NWT does not deactivate the model cells with low water level resulting in the simulation of groundwater flow, without the numerical error that is inserted with the assumption of dry cells by the MODFLOW 2000 code. In addition to the simulation of the groundwater flow and the evolution of all simulation and analysis tools that occur over the years, it is important to mention the difference in the groundwater levels. The simulated water levels of the historical period of 1987–2009 agree well with the simulation results of the previous studies. Furthermore, the groundwater balance indicated an improving trend during the period 2009–2018 but it remains under pressure and intensive and extensive groundwater abstractions. The varying timesteps of the simulation according to the varying dates of groundwater measurements depict clearly the actual and real conditions of the groundwater with increased timely and spatial accuracy.

This study analyzes the groundwater resources status of the Lake Karla aquifer focusing on a precise simulation–validation procedure of groundwater model based upon new and innovative methodologies. The simulated annual decrease in groundwater levels verify the groundwater depletion and shows that the groundwater system is still under intense pressure even though farming activity in recent years has turned to less water-intensive crops. Moreover, the enduring water balance deficit indicates the necessity for the construction of collective irrigation networks that will eliminate the extravagant number of private own wells, and the need for sustainable water resources management. Unless appropriate measures and policies are implemented, the Lake Karla's aquifer system may face the problem of water scarcity.

This article presents a method that uses a configurable number of pilot points to calibrate and subsequently analyze the groundwater flow system. The use of the Newton smoothing technique in solving the inverse problem with the use of pilot points ensures stability in hydraulic conductivity parameter values and the other hydrogeological variables for Lake Karla's aquifer. Furthermore, the use of tools like PLPROC and mod2obs contribute significantly to the creation of a realistic image of the hydrogeological environment of the Karla aquifer. The pilot points approach is particularly important in data sparse areas of hydrogeological information. Hence, the proposed approach could be used in conjunction with stochastic field generation techniques for assessing model uncertainty and especially in contaminant fate and transport modelling studies.

All relevant data are included in the paper or its Supplementary Information.