Groundwater flow simulation through the application of the FREEWAT modeling platform

Groundwater modeling is an efficient tool for education, hydrologic research, water management, groundwater protection, and remediation. Numerical models are a simplified representation of the operation of a real ground-water system with mathematical equations solved by a computer program (Reilly & Harbaugh 2004; Harbaugh 2005).

Over the years, many numerical groundwater modeling software based on different methods has been developed and used widely. In relation to code availability, there are two categories: public domain (open-source) and proprietary or commercial (closed-source). Public domain codes are free to use and can be readily accessed, reviewed, and modified while proprietary codes are only available for use after purchase. Examples of public domain software are: MODFLOW (USGS), SUTRA (USGS), VS2DI (USGS), MT3DMS (University of Alabama), PHAST (USGS), RT3D (Pacific Northwest Laboratory), TOUGH2 (Lawrence Berkeley National Laboratory), ParFLOW (Lawrence Livermore National Laboratory), and STOMP (Pacific Northwest National Laboratory). Examples of commercial codes are: MODHMS (Hydrogeologic Inc.), MODFLOW-SURFACT (Waterloo Hydrologic), FEFlow (Wasy Ltd), and MikeSHE (DHI Water and Environment).

MODFLOW is a physically based, spatially distributed code developed by USGS, which simulates steady and non-steady flow in an irregularly shaped flow system in which aquifer layers can be confined, unconfined, or a combination of confined and unconfined (Harbaugh 2005). Flow from external stresses, such as flow to wells, areal recharge, evapotranspiration, flow to drains, and flow through river beds, can be simulated. MODFLOW includes the main program and independent subroutines (packages) that define a specific characteristic of the hydrological system that is to be simulated. Separating the program into packages allows the user to investigate a particular hydrologic feature of the model independently.

Several GUIs (Graphical User Interface) for MODFLOW have been developed. These are classified as commercial GUIs such as Groundwater Modeling System, GMS (U.S. Department of Defense), Visual MODFLOW (Waterloo Hydrogeologic Inc.), Groundwater Vistas (Environmental Simulation, Inc.), and open source GUIs such as ModelMuse (USGS), ModelMate (USGS), Processing MODFLOW for Window, and FREEWAT (Rossetto et al. 2018).

Integrating geographic information system (GIS) and hydrological codes facilitates the use of complex modeling environments, allowing to store, manage, and visualize large spatial datasets as stated by relevant EU regulations and recommendations.

FREEWAT (Rossetto et al. 2018) aims at promoting the application of EU water-related directives by creating a public domain, QGIS-integrated platform, developed to simulate several hydrological processes in order to address decision-making in water resources management (De Filippis et al. 2016).

The FREEWAT platform is a QGIS (QGIS Development Team 2009) plug-in, which combines the abilities of GIS for geo-processing and post-processing tools with several codes of the MODFLOW, USGS family for the simulation of the hydrological cycle, hydrochemical or economic-social processes. Input and output data are managed through the SpatiaLite (SpatiaLite Development Team 2011) Data Base Management System (De Filippis et al. 2016). FloPy is used as reference Python library to connect with hydrological codes.

The FREEWAT platform includes the following modules for pre-processing and model implementation:

• A pre-processing tool allowing to import, analyze, and visualize time series data which can be used for model calibration.

• MODFLOW_2005 (Harbaugh 2005) to perform groundwater flow simulation in the saturated and unsaturated zones.

• A module for sensitivity analysis, calibration, and parameter estimation using the UCODE_2014 (Poeter et al. 2014) which improves the model fit, by reducing the difference between model-simulated heads and flows and the observed data.

• Tools for the analysis, interpretation, and visualization of hydrogeological and hydrochemical data (Criollo et al. 2019).

• A module to simulate solute transport in groundwater flow systems using MT3DMS (Zheng & Wang 1999).

• A module for water management and planning of rural environments by integrating MODFLOW-OWHM (One-Water Hydrologic Flow Model; Hanson et al. 2014).

The FREEWAT platform was tested in 14 case studies in EU and non-EU countries aiming at specific water resources issues (Dadaser-Celik & Celik 2017; Grodzynskyi & Svidzinska 2017; Kopač & Vremec 2017; Panteleit et al. 2017; Perdikaki et al. 2017; Positano & Nannucci 2017; Cannata et al. 2018). These case studies are divided into two categories: (i) nine case studies (eight in EU countries and one in Switzerland) are specifically referred to meet the requirements of the Water Framework Directive, Ground Water Directive, and other water-related Directives for water resource management; (ii) five case studies (two in EU countries, one in Ukraine, one in Turkey, and one in Africa) are devoted to address water management issues mostly related to the rural environment.

The objective of the paper is to demonstrate the application of the FREEWAT platform for the groundwater flow simulation of a sedimentary aquifer, located in northern Greece, at the south-west coastal part of the Prefecture of Rhodope. The study area was selected due to the existence of a large amount of hydrogeologic, hydrologic, geological, and meteorological data from previous studies (Diamantis et al. 1994; Petalas 1997; Kallioras et al. 2006a, 2006b, 2006c, 2010; Kallioras 2008).

The study area is located in northern Greece, at the south-west coastal part of the Prefecture of Rhodope and extends between the lakes Ismarida and Vistonida from east to west (Figure 1). The total area of investigation covers approximately 160 km², including the populated areas which occupy a very small part (approximately 1.6% of the total area). Most of the area is almost exclusively used for agricultural activities and is being extensively irrigated from groundwater resources pumped from the investigated aquifer system through a dense groundwater wells network (Figure 1) (Petalas 1997; Kallioras et al. 2006a, 2006b, 2006c; Kallioras 2008).

The surface relief of the study area is characterized as hilly to semi-hilly (70 m being the highest altitude) and most of the study area presents slope ranges which do not exceed 20% (Diamantis et al. 1994). The geomorphological relief contains a significant hydrographic network which includes Rivers Vosvozis and Aspropotamos and contributes to the recharge conditions of the study aquifer system (Diamantis et al. 1994). At the southern part, a series of coastal lakes appear which are hydraulically disconnected from the study aquifer system (Kallioras et al. 2006b).

The geological and hydrogeological conditions of the study area have been investigated by previous studies (Stamatakis 1992; Diamantis et al. 1994; Petalas 1997; Kallioras et al. 2006a, 2006b, 2006c, 2010; Kallioras 2008). The study area is part of the Rhodope massif and includes Quaternary, Pliocene, Upper Miocene, and Paleocene deposits. The Quaternary deposits appear at the southern part and at the northern parts of the area of investigation, and due to their high permeability they are of great importance for the recharge conditions of the study area. Pliocene deposits generally include coarse-grained materials, which extend throughout the study area and Upper Miocene formations consist of high permeability materials. The aforementioned materials are interrupted by the presence of gray-green clay, which negatively influences the vertical movement of groundwater, as it frequently appears in different layers throughout the entire column of the area. The Paleocene formations include mainly clays and sandstones and are located at the northern and eastern boundaries of the study area. Two major normal faults of NE–SW direction appear in the Upper Miocene deposits from the east to the west, which form the boundary conditions of the study area (Figure 2).

The hydraulic parameters of the study aquifer system have been investigated in previous studies (Stamatakis 1992; Diamantis et al. 1994; Petalas 1997). The points where the pumping tests were carried out are shown in Figure 1. The values of hydraulic conductivity, K, range between 2.49 × 10⁻⁵ and 2.22 × 10⁻³ m/s, with an average value of 4.67 × 10⁻⁴ m/s. The values of transmissivity, T, range between 1.15 × 10⁻³ and 9.971 × 10⁻² m²/s, with an average value of 2.395 × 10⁻² m²/s and the values of storativity, S, range between 1.3 × 10⁻⁴ and 2.9 × 10⁻² with an average value of 5.06 × 10⁻³. The spatial distribution of the values of hydraulic conductivity is shown in Figure 3, where it is observed that the highest values of hydraulic conductivity appear in the central part of the study area, which is correlated to the intensive groundwater exploitation of that certain region.

The climate of the study area is characterized as typical Mediterranean, with warm and dry summers as well as mild and wet winters. According to data retrieved by the rainfall gauge of Porpi, located at the center of the study area, for the period between 1954 and 2005, the average annual precipitation is 555 mm (Figure 4), November and December being the months where highest precipitation values have been recorded and July and August the months with highest temperatures and lowest precipitation heights (Kallioras et al. 2006b).

The potentiometric conditions prior to and after the irrigation period of the study area are shown in the piezometric maps, Figures 5 and 6, respectively (April and October 2003). The main recharge axis from the adjacent northern aquifer system of the alluvial cone of Kompsatos and the formation of the deep cone of depression (−35 to −40 m groundwater levels) in the central part of the area, due to overexploitation of the groundwater resources of the aquifer, are observed.

The conceptual model was developed by the collection of a large amount of relevant hydrogeologic, hydrologic, geological, and meteorological data of the study area (Diamantis et al. 1994; Petalas 1997; Kallioras et al. 2006a, 2006b, 2006c, 2010; Kallioras 2008). The recharge of the aquifer system is achieved by direct infiltration from precipitation, irrigation returns, percolation from surface water bodies of the area, and the lateral hydraulic contact with the northern alluvial cone of Kompsatos River (Kallioras 2008).

The direct infiltration from precipitation was estimated to be approximately 10% as reported by Kallioras (2008) and the precipitation data which were used as model inputs concern the period April–October 2003. The precipitation data were retrieved from the rainfall gauge of Porpi, located at the centre of the study area. According to Kallergis (2000), generally 30% of the amount of groundwater which is abstracted for irrigation returns to the unconfined aquifers as irrigation returns.

The hydrogeological boundaries of the model were determined according to the geological map of the area (Kallioras 2008):

• The upper geological layer (approximately 5 m) overlying the main aquifer includes impermeable formations (siltstones, sandstones, conglomerate, clays) and provides the confining conditions of the aquifer system.

• The middle geological formation (approximately 50 m) which is considered as the confined aquifer of investigation, consists of permeable alluvial deposits, mainly sands, gravels, and pebbles.

• The bottom of the aquifer system is composed of a gray-green clay layer, which extends up to the whole area of investigation.

The model grid network of the study area has dimensions of 35 rows and 40 columns and each cell has dimensions of 400 × 400 m (Figure 7).

The term boundary conditions represent any flow or head constraints within the flow domain, with respect to the hydrogeological conditions which occur at the lateral boundaries of the investigated aquifer system (Harbaugh 2005).

The boundary conditions and source/sink terms of the study area are represented through the following MODFLOW packages:

• General head boundaries (GHB) package for the parts of the aquifer which are in hydraulic connection with the alluvial cone of Kompsatos River, at parts of Vistonida lake inlet and at the southern part of Ismarida Lake.

• River boundaries (RIV) package across the Vozvozis River.

• Well (WELL) package to simulate wells that withdraw water from the aquifer at a constant rate during a stress period.

• Recharge (RCH) package to simulate the direct infiltration from precipitation.

• Evapotranspiration (EVT) package to simulate the conjunctive effects of plant transpiration and direct evaporation from groundwater.

• No flow boundaries at the parts of the aquifer where no flow occurs (inflow or outflow).

MODFLOW_2005 (Harbaugh 2005) provides the ability to compare simulated heads and flows to measured data values, which is called the Observation Process (OBS), and OBS in MODFLOW_2005 is derived from OBS in MODFLOW_2000 (Hill et al. 2000). The OBS is the first step for the calibration of a model and the application of the sensitivity analysis for FREEWAT. The types of observations which could be used in the OBS process include: hydraulic heads; changes in hydraulic head over time; flows to or from surface-water bodies represented using the General head boundary (GHB), Drain (DRN), River (RIV), or Streamflow Routing (STR) Packages; and flow to or from constant-head cells (Foglia et al. 2017).

Head observations are groundwater head measurements in a monitoring well. The HOB package is used to compare calculated data values, derived from a model, to measured values acquired from field measurements. Hydraulic heads can be defined at any location and time. Temporal and spatial interpolation (horizontally and vertically) is applied when an observation does not correspond directly to a node (Hill et al. 2000; Foglia et al. 2017).

Information about the location of the observation, date, and time of the individual observation, the observation value as well as statistical weighting method and weighting value, is required for the implementation of head observations. For this purpose comma-separated value files (csv-files) or sensor data using the FREEWAT OAT tool can be used.

The OAT is a (FREEWAT module) pre-processing tool which facilitates the import analysis and visualization of time series data in order to enhance model development and model calibration (De Filippis et al. 2017).

The OAT library consists of a sensor stored in the SpatiaLite database representing a single point of observation which is composed of time series data and metadata. Sensors' metadata are used to identify the point of observation including the name, location (latitude, longitude, elevation), observed property, unit of measure, coordinate system, time zone, weight statistic, and data availability (time interval). The method represents a processing method that can be applied to a sensor. Data entry can be accomplished from data stored in local files or databases into the FREEWAT GIS environment, MODFLOW model results (from list, HOB, or gage files) and sensor data from an istSOS server (Cannata et al. 2017). Sensor objects can be plotted, exported as new CSV files, or used to create observation input files. The head observation layer (HOB) created by OAT, can be used for model calibration in FREEWAT using UCODE (Poeter et al. 2014).

The simulation involved the period between April and October 2003 and the model was divided into two stress periods. The first stress period (steady-state) lasted 1 day, during which no stresses (recharge, pumping wells, etc.) affected the hydrogeological system. The second stress period (transient) lasted 181 days, during which 411 pumping wells penetrating the deepest layer were activated. The effects on the simulated hydraulic head and the aquifer water budget were assessed at the end of the second stress period. The calendar periods between the stress periods and the time steps are shown in Table 1.

The model inflows during the simulation period include the following:

• Lateral groundwater inflows from the alluvial cone of Kompsatos River (General head boundary conditions at the NW boundaries of the conceptual model, across Aspropotamos torrent).

• Lateral groundwater inflows from the quaternary aquifer system of Neo Sidirochori (General head boundary conditions at the NE boundaries of the conceptual model).

• Saline water inflow at Vistonida inlet and at the southern part of Ismarida Lake (General head boundary conditions at the W and E boundaries of the conceptual model).

• Percolation from Vozvozis River (River boundary conditions across the Vozvozis River at the E part of the conceptual model).

• Irrigation returns (Well boundary conditions – positive flow rates – of the conceptual model).

• Distributed recharge based on rainfall rate (Recharge boundary conditions of the conceptual model).

The model outflows for the same simulation period are associated with the following:

• Groundwater pumping during the irrigation period (Well boundary conditions – negative flow rates – of the conceptual model).

• Distributed evapotranspiration loss from the water table (Evapotranspiration boundary conditions of the conceptual model).

• Lateral groundwater outflows from the alluvial cone of Kompsatos River, at Vistonida inlet and at the southern part of Ismarida Lake (General head boundary conditions at the NW, W, and E boundaries of the conceptual model).

The total amount of water, which inflows and outflows into the aquifer system accounts for 39,937,380 m³ and 39,937,388 m³, respectively. To evaluate the calibration results and compare the calculated groundwater levels to the measured ones, the HOB was used, where 41 wells were used as observation wells (Figure 8). For the evaluation of model fit performance, the Nash–Sutcliffe efficiency index and the mass balance error were used.

The Nash–Sutcliffe coefficient and the mass balance error accounted for 0.2 and 29.4, respectively, indicating that the simulated level values are greater than the observed values. Figure 9 presents the observed plotted versus the simulated groundwater levels and Figure 10 presents the difference between these values (residuals). Both figures refer to the measurement of 13 October 2003 in order to illustrate the correlation between the two different values. Table 2 shows the model fit statistics for the 41 observations used to estimate parameters for the calibration period.

Figure 11 shows a piezometric map with the calculated values for the measurement of 13 October 2003. From Figures 7 and 10, it is concluded that at NE and central boundaries of the study area a coincidence of observed and simulated groundwater levels occurs. Declination between the simulated and the observed groundwater levels is observed at the western and NW parts of the investigated aquifer system.

The UCODE code is a public domain, open-source software supported by U.S. Geological Survey. UCODE is one of a set of inverse modeling codes, developed for models in which the number of parameters is less than the number of observations. FREEWAT incorporates UCODE_2014 for sensitivity and uncertainty analysis (Poeter et al. 2005; Poeter et al. 2014).

FREEWAT provides integration of UCODE_2014 (Hill & Tiedman 2007) for the sensitivity analysis and parameter estimation. A series of fit independent statistics was used for this paper: composite scaled sensitivities (CSS), parameter correlation coefficients (PCC), and leverage.

Composite scaled sensitivities (CSS) are an averaged measure of all DSS for a single parameter and indicate the importance of observations as a whole to a single parameter. Larger CSS values point out parameters for which more information is provided by all observations for each parameter.

After the first run of the model, a sensitivity analysis was performed using the UCODE_2014 code (Poeter et al. 2014) to gain information on the quantitative relationship between parameters and observations. For the sensitivity analysis, a variety of parameters representing the important processes was included to cover all possible impacts in the model setup. The parameters taken into consideration were hydraulic conductivity (HK), specific storage (SS), and direct infiltration from precipitation (RCH). Hydraulic conductivity was divided into four zones (ΗΚ1–ΗΚ4) (Figure 12). Table 3 presents the initial values of those parameters.

The statistical analysis and the evaluation of the results were carried out using the following fit-independent methods: dimensionless scaled sensitivity, composite scaled sensitivity, parameter correlation coefficient, and leverage.

Table 4 shows that the parameters HK1, RCH, and HK2 have the highest CSS value and are therefore the most sensitive parameters. This is illustrated in Figure 13. Figure 14 shows the DSS values for all parameters for each observation. Regarding the parameter correlation coefficient, the parameters HK1 and RCH have a factor of PCC = 0.91 and the parameters HK3 and HR4 have a factor of PCC = 0.90 (Figure 15). Figure 16 shows that the most influential observations in the calculation of a parameter are wel_21 and wel_35 related to parameter HK1 and wel_42 for the RCH parameter.

The parameter estimation was conducted with the same code, to obtain the set of parameters which give the best fit between observed and simulated values of the square weighted residuals regression (Foglia et al. 2017). The objective function was defined as the sum of squared weighted residuals.

In the parameter estimation process, the three parameters with the highest CSS values were used. The parameters SS, HK3, and HK4 were excluded from further estimation because of their low sensitivity. The high correlation between HK1 and RCH was addressed by imposing constraints based on existing data (Kallioras 2008). The initial and final values of the parameters are listed in Table 5.

The Nash–Sutcliffe coefficient and the mass balance error accounted for 0.74 and 0.8, respectively. Table 6 shows the model fit statistics for the 41 observations after the successful parameter estimation.

The final results of the calibrated model (observed plotted versus the simulated groundwater levels) are presented in Figure 17. Figure 18 shows the difference between the measured and calculated groundwater levels (residuals). A piezometric map with the final simulated values for the measurement of 13 October 2003 is presented in Figure 19.

The FREEWAT platform was applied to an aquifer located in northern Greece. The groundwater flow simulation was achieved by using the MODFLOW_2005 code (and packages) and the OAT. The OAT module was used to create the head observation layer where 41 wells were used as monitoring wells.

After the first run of the model, a sensitivity analysis was performed using the UCODE_2014 code (Poeter et al. 2014) to gain information on the quantitative relationship between parameters and observations. The parameters taken into consideration were: hydraulic conductivity (HK) which was divided into four zones, specific storage (SS), and direct infiltration from precipitation (RCH). The statistical analysis of the results was carried out using a series of fit independent statistics (Hill & Tiedman 2007) such as: dimensionless scaled sensitivity, composite scaled sensitivity, parameter correlation coefficient, and leverage. The most important parameters with the highest CSS values affecting the results were the first zone of hydraulic conductivity of the aquifer (HK1), followed by the direct infiltration from precipitation (RCH) and the second zone of hydraulic conductivity (HK2) (Figure 13). Figure 14 shows that the observations that provide most information on the parameter HK1 are wel_14, wel_12, and wel_6. Furthermore, the most influential observations in the calculation of a parameter (Figure 16) are wel_21 and wel_35 related to parameter HK1 and wel_42 for the RCH parameter.

In the parameter estimation process, which was conducted with the same code, the three parameters with the highest CSS values were used. The small CSS values of the parameters SS, HK3, and HK4 in Figure 13 and the PCC absolute values close to 1.00 in Figure 15 indicate that attempts to estimate all of the parameters may produce failed regressions. Therefore, the parameters SS, HK3, and HK4 were excluded from further estimation because of their low sensitivity. In addition, the high correlation between HK1 and RCH was addressed by imposing constraints based on existing data.

After the parameter estimation, the Nash–Sutcliffe coefficient and the mass balance error accounted for 0.74 and 0.8, respectively, indicating that the simulated level values are slightly greater than the observed values. The calculated values were compared with the measured values of 13 October 2003, and it was concluded that at the eastern, north-eastern, western, north-western, and southern boundaries of the study area a coincidence of observed and simulated groundwater levels occurred. On average, the model results were within 4.15 m of the measured heads. Small deviations between the simulated and the observed groundwater levels were observed at the central part of the area at the point of formation of the cone of depression (wel_33), at the south-west part of Vistonida Lake (wel_44, wel_11) and at the south part of Aspropotamos River (wel_43).

Although the model provided satisfactory results, model performance could be improved by further developments:

• New field measurements and pumping tests should be made to review the hydraulic characteristics and the hydrogeological conditions of the study area due to heterogeneity of the alluvial formations which compose the aquifer layer and the large fluctuation of the groundwater levels in most of the monitoring wells.

• Further investigation into the hydraulic connection of the aquifer system with the sea, Vistonida and Ismarida lakes, and the inaccurate aquifer discharge towards the southern part of the study area.

• New records of the abstraction wells and crops of the area to estimate the necessary quantities of irrigation water.

• Review of the spatial and temporal discretization of the model to represent with greater accuracy the hydraulic characteristics of the area and to assure that changes in smaller time scales can be simulated and be distinguished in model results.

• Regarding the improvement of FREEWAT capabilities, there is the need of integrating more MODFLOW solver packages, the implementation of stochastic simulations methods, and modules for improving simulation of mass exchange between underground and surface water (Rossetto et al. 2018).

The case study demonstrates that by using the FREEWAT plug-in directly in the QGIS environment the user can easily archive, pre-process, and analyze large datasets, build a set of models, and post-process results in a unique software. Furthermore, the OAT tool helps simplify the processing of time series data in order to facilitate the water quality and quantity management.

As a final remark, it should be noted that the codes integrated in FREEWAT are released under compatible free and open-source licenses; a key strategy adopted by the project in making the outcomes sustainable and easily adopted (Cannata et al. 2018). The openness of a code is an important factor in scientific analysis as it promotes the reliability of the analysis performed and fast development of the code (Dile et al. 2016; Rossetto et al. 2018). Modeling codes that are not open or free restrict the usage due to the high cost. In addition, commercial GISs are connected with high licensing costs and lead to concerns regarding rigidity in data-models and platform dependence (Bhatt et al. 2014). Therefore, free and open source GIS-integrated software tools may contribute to enhance groundwater management capabilities from a technical point of view (Dile et al. 2016; Wang et al. 2016; De Filippis et al. 2017; Rossetto et al. 2018). Nevertheless, in order to achieve a better sustainability, there is the need of supporting open data to promote water-related studies.