Abstract

FREEWAT is a free and open source QGIS-integrated platform, developed to simulate several hydrological processes by combining the capabilities of geographic information system (GIS) for geo-processing and post-processing tools with several codes of the well-known USGS MODFLOW ‘family’. FREEWAT platform was applied for the groundwater flow simulation of a coastal aquifer system, located in northern Greece. The simulation was conducted using the MODFLOW_2005 code, the Observation Analysis Tool (a FREEWAT module facilitating the integration of time series observations into modeling), while the UCODE_2014 code was used as the main module for the sensitivity analysis and parameter estimation. The statistics used include composite scaled sensitivities, parameter correlation coefficients, and leverage. The simulation of the investigated aquifer system was found to be satisfactory, indicating that the simulated level values were slightly greater than the observed values after the optimization.

INTRODUCTION

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).

MATERIALS AND METHODS

Characterization of the study area and conceptual model

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 km2, 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).

Figure 1

Geographical location of the study area and monitoring points for groundwater level measurements and estimation of aquifer hydraulic parameters.

Figure 1

Geographical location of the study area and monitoring points for groundwater level measurements and estimation of aquifer hydraulic parameters.

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).

Figure 2

Geological map of the study area (adapted from Petalas 1997).

Figure 2

Geological map of the study area (adapted from Petalas 1997).

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−5 and 2.22 × 10−3 m/s, with an average value of 4.67 × 10−4 m/s. The values of transmissivity, T, range between 1.15 × 10−3 and 9.971 × 10−2 m2/s, with an average value of 2.395 × 10−2 m2/s and the values of storativity, S, range between 1.3 × 10−4 and 2.9 × 10−2 with an average value of 5.06 × 10−3. 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.

Figure 3

Distribution of hydraulic conductivity K (×10−5 m/s; monitoring points of Figure 1).

Figure 3

Distribution of hydraulic conductivity K (×10−5 m/s; monitoring points of Figure 1).

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).

Figure 4

Annual precipitation for the period 1954–2005.

Figure 4

Annual precipitation for the period 1954–2005.

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.

Figure 5

Piezometric map April 2003 (monitoring points of Figure 1).

Figure 5

Piezometric map April 2003 (monitoring points of Figure 1).

Figure 6

Piezometric map October 2003 (monitoring points of Figure 1).

Figure 6

Piezometric map October 2003 (monitoring points of Figure 1).

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.

Model development and implementation

Model grid network

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).

Figure 7

Conceptual model and grid network of the study area.

Figure 7

Conceptual model and grid network of the study area.

Boundary conditions and source/sink terms

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).

Observations

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 observation package (HOB)

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.

Observation analysis tool (OAT)

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).

Model stress periods

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.

Table 1

Stress periods and model time steps

Stress period (SP) Time step (TS) Calendar period 
SP1 TS1 15/04/2003 
SP2 TS1 16/04/2003–16/05/2003 
TS2 17/05/2003–16/06/2003 
TS3 17/06/2003–16/07/2003 
TS4 17/07/2003–16/08/2003 
TS5 17/08/2003–16/09/2003 
TS6 17/09/2003–13/10/2003 
Stress period (SP) Time step (TS) Calendar period 
SP1 TS1 15/04/2003 
SP2 TS1 16/04/2003–16/05/2003 
TS2 17/05/2003–16/06/2003 
TS3 17/06/2003–16/07/2003 
TS4 17/07/2003–16/08/2003 
TS5 17/08/2003–16/09/2003 
TS6 17/09/2003–13/10/2003 

Calibration results and model balance

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 m3 and 39,937,388 m3, 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.

Figure 8

Head observation wells created with the observation analysis tool (monitoring points of Figure 2).

Figure 8

Head observation wells created with the observation analysis tool (monitoring points of Figure 2).

The Nash–Sutcliffe efficiency index, NS, is calculated as (Nash & Sutcliffe 1970):  
formula
(1)
where qsim and qobs are the simulated and the observed values. An efficiency of NS = 1 corresponds to a perfect match between model and observations. Essentially, the closer the model efficiency is to 1, the more accurate the model is.
The mass balance error (m) is calculated as:  
formula
(2)
where qsim and qobs are the simulated and the observed values. The value of m might be positive, indicating that the simulated values are larger than the observed values or negative, indicating that the observed values are larger than the simulated values. When the observed values match the simulated ones, m = 0.

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.

Table 2

Model fit performance

Parameter Value 
Absolute residual mean 7.827 
Standard error of the estimate 0.929 
Residual RMS 7.151 
Normalized RMS 0.168 
Pearson correlation coefficient 0.893 
Parameter Value 
Absolute residual mean 7.827 
Standard error of the estimate 0.929 
Residual RMS 7.151 
Normalized RMS 0.168 
Pearson correlation coefficient 0.893 
Figure 9

Simulated versus observed groundwater levels for the measurement of 13 October 2003.

Figure 9

Simulated versus observed groundwater levels for the measurement of 13 October 2003.

Figure 10

Residuals of the second stress period (observations of Figure 8).

Figure 10

Residuals of the second stress period (observations of Figure 8).

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.

Figure 11

Simulated piezometric map October 2003 (observations of Figure 8).

Figure 11

Simulated piezometric map October 2003 (observations of Figure 8).

Sensitivity analysis

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.

Fit-independent statistics

Composite scaled sensitivities (CSS) are calculated using dimensionless scaled sensitivities (DSS). Dimensionless scaled sensitivities express the importance of an observation to the evaluation of a parameter. Larger absolute values of DSS indicate greater importance. Dimensionless scaled sensitivities for observation i and parameter j are calculated as:  
formula
(3)
where i is the observation index, is the simulated value being compared to the ith observation, is the jth estimated parameter, is the sensitivity of the ith simulated value with respect to the jth parameter, b is a vector of parameter values, is the weight for the jth observation, n is the number of observations, and p is the number of parameters.

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.

Composite scaled sensitivity for any parameter j is calculated as:  
formula
(4)
where is from Equation (3).
Parameter value non-uniqueness can be detected using parameter correlation coefficients (PCC). If the absolute value of PCC is greater than 0.95 the calibration data cannot be used to independently estimate a pair of parameters. PCC is calculated as:  
formula
(5)
where are the covariances; are the variances of the square, symmetric parameter variance covariance matrix, V(b).
V(b) is calculated as:  
formula
(6)
where is the calculated error variance (CEV); Χ is an n by p matrix of sensitivities defined after Equation (3); and ω is an n by n weight matrix.
Leverage represents the potential influence an observation has on the estimation of a unique parameter. Values of leverage range from 0.0 to 1.0. Leverage is calculated as:  
formula
(7)
where hii is the leverage for observation I; xi is the transpose of the ith row of the sensitivity matrix defined after Equation (6).

Parameter estimation

Parameter estimations are performed 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). Parameter estimation refers to the process of optimizing a selected set of parameters using a nonlinear regression which minimizes the objective function. The objective function is defined as the sum of squared weighted residuals (Hill & Tiedman 2007):  
formula
(8)
where S(b) is the objective function, b the selected model parameter, ωi the observation weight, hi and the observed and simulated head, and are the observed and simulated flow, and and are the observed and simulated prior information, respectively.

RESULTS

Sensitivity analysis

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.

Table 3

Initial values of the model parameters (source of initial estimate: Stamatakis 1992; Diamantis et al. 1994; Petalas 1997; Kallioras 2008)

Parameter Initial value Parameter Initial value 
HK1 (m/day) HK4 (m/day) 80 
HK2 (m/day) 25 SS (m−12.6 × 10−5 
HK3 (m/day) 40 RCH (m/day) 0.0002 
Parameter Initial value Parameter Initial value 
HK1 (m/day) HK4 (m/day) 80 
HK2 (m/day) 25 SS (m−12.6 × 10−5 
HK3 (m/day) 40 RCH (m/day) 0.0002 
Figure 12

Hydraulic conductivity zones (monitoring points of Figure 1).

Figure 12

Hydraulic conductivity zones (monitoring points of Figure 1).

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.

Table 4

Composite scaled sensitivity of the model parameters

Parameter CSS Parameter CSS 
HK1 0.064 HK4 0.005 
HK2 0.031 SS 8.425 × 10−7 
HK3 0.012 RCH 0.036 
Parameter CSS Parameter CSS 
HK1 0.064 HK4 0.005 
HK2 0.031 SS 8.425 × 10−7 
HK3 0.012 RCH 0.036 
Figure 13

Parameter importance to observations, based on composite scaled sensitivity (observations of Figure 8).

Figure 13

Parameter importance to observations, based on composite scaled sensitivity (observations of Figure 8).

Figure 14

Observation importance based on dimensionless scaled sensitivity of the model parameters (observations of Figure 8).

Figure 14

Observation importance based on dimensionless scaled sensitivity of the model parameters (observations of Figure 8).

Figure 15

Parameter dependence based on parameter correlation coefficient (observations of Figure 8).

Figure 15

Parameter dependence based on parameter correlation coefficient (observations of Figure 8).

Figure 16

Observation dominance based on leverage (observations of Figure 8).

Figure 16

Observation dominance based on leverage (observations of Figure 8).

Parameter estimation

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.

Table 5

Results of the final parameter estimation (source of initial estimate: Stamatakis 1992; Diamantis et al. 1994; Petalas 1997; Kallioras 2008)

Parameter (units) Initial value Calibrated value 
HK1 (m/day) 1.3 
HK2 (m/day) 25 15 
HK3 (m/day) 40 40 
HK4 (m/day) 80 80 
SS (m−10.000026 0.000026 
RCH (m/day) 0.0002 0.0011 
Parameter (units) Initial value Calibrated value 
HK1 (m/day) 1.3 
HK2 (m/day) 25 15 
HK3 (m/day) 40 40 
HK4 (m/day) 80 80 
SS (m−10.000026 0.000026 
RCH (m/day) 0.0002 0.0011 

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.

Table 6

Model fit performance after the final parameter estimation

Parameter Value 
Absolute residual mean 4.188 
Standard error of the estimate 0.811 
Residual RMS 0.421 
Normalized RMS 0.009 
Pearson correlation coefficient 0.874 
Parameter Value 
Absolute residual mean 4.188 
Standard error of the estimate 0.811 
Residual RMS 0.421 
Normalized RMS 0.009 
Pearson correlation coefficient 0.874 

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.

Figure 17

Simulated versus observed groundwater levels for the measurement of 13 October 2003 after the parameter estimation (observations of Figure 8).

Figure 17

Simulated versus observed groundwater levels for the measurement of 13 October 2003 after the parameter estimation (observations of Figure 8).

Figure 18

Residuals of the second stress period after the final parameter estimation (observations of Figure 8).

Figure 18

Residuals of the second stress period after the final parameter estimation (observations of Figure 8).

Figure 19

Simulated piezometric map October 2003 after the final parameter estimation (observations of Figure 8).

Figure 19

Simulated piezometric map October 2003 after the final parameter estimation (observations of Figure 8).

DISCUSSION AND CONCLUSIONS

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.

REFERENCES

REFERENCES
Bhatt
G.
,
Kumar
M.
&
Duffy
C. J.
2014
A tightly coupled GIS and distributed hydrologic modeling framework
.
Environ. Model. Softw.
62
,
70
84
.
Cannata
M.
,
Neumann
J.
&
Cardoso
M.
2017
FREEWAT User Manual, Volume 5 – Observation Analysis Tool, version 0.4, March 31st, 2017
.
Cannata
M.
,
Neumann
J.
&
Rossetto
R.
2018
Open source GIS platform for water resource modelling: FREEWAT approach in the Lugano Lake
.
Spat. Inf. Res.
26
(
3
),
241
251
.
doi: 10.1007/s41324-017-0140-4
.
Criollo
R.
,
Velasco
V.
,
Nardi
A.
,
Vries
L. M.
,
Riera
C.
,
Scheiber
L.
,
Jurado
A.
,
Brouyère
S.
,
Pujades
E.
,
Rossetto
R.
&
Vázquez-Suñé
E.
2019
AkvaGIS: an open source tool for water quantity and quality management
.
Comput. Geosci.
127
,
123
132
.
doi:10.1016/j.cageo.2018.10.012
.
Dadaser-Celik
F.
&
Celik
M.
2017
Modelling surface water-groundwater interactions at the Palas Basin (Turkey) using FREEWAT
.
Acque Sotterranee – Ital. J. Groundwater
6
(
3/149
),
53
60
.
doi:10.7343/as-2017-288
.
De Filippis
G.
,
Borsi
I.
,
Foglia
L.
,
Cannata
M.
,
Velasco
V.
&
Rossetto
R.
2016
The H2020 FREEWAT Project for Developing a GIS-integrated platform for water resource management
. In:
Convegno Nazionale di Idraulica E Costruzioni Idrauliche
,
Bologna, Italy
.
De Filippis
G.
,
Borsi
I.
,
Foglia
L.
,
Cannata
M.
,
Velasco
V.
,
Vasquez-Suñe
E.
,
Ghetta
M.
&
Rossetto
R.
2017
Software tools for sustainable water resources management: the GIS-integrated FREEWAT platform
.
Rend. Online Soc. Geol. Ital.
42
,
59
61
.
doi: 10.3301/ROL.2017.14
.
Diamantis
I
,
Petalas
C.
,
Tzevelekis
T.
&
Pliakas
F.
1994
Investigation for the Possibility of Water Supply of Coastal Settlements of Thrace From Coastal Aquifers
.
Technical Report for the Region of Eastern Macedonia and Thrace
,
Greece
,
4
, p.
393
.
Dile
Y. T.
,
Daggupati
P.
,
George
C.
,
Srinivasan
R.
&
Arnold
J.
2016
Introducing a new open source GIS user interface for the SWAT model
.
Environ. Model. Softw.
85
,
129
138
.
Foglia
L.
,
Toegl
A.
&
Mehl
S.
2017
FREEWAT User Manual, Volume 6 - MODFLOW OBServation process, version 0.4, March 31st, 2017
.
Grodzynskyi
M.
&
Svidzinska
D.
2017
Modelling the impact of rural land use scenarios on water management: a FREEWAT approach to the Bakumivka catchment case study, Ukraine
.
Acque Sotterranee – Ital. J. Groundwater
6
(
3/149
),
39
50
.
doi: 10.7343/as-2017-291
.
Hanson
R. T.
,
Boyce
S. E.
,
Schmid
W.
,
Hughes
J. D.
,
Mehl
S. M.
,
Leake
S. A.
,
Maddock
T.
&
Niswonger
R. G.
2014
One-Water Hydrologic Flow Model (MODFLOW-OWHM)
.
Techniques and Methods 6-A51
,
U.S. Geological Survey, Reston
, VI, p.
134
.
Harbaugh
A. W.
2005
MODFLOW-2005, The U.S. Geological Survey Modular Ground-Water Model – the Ground-Water Flow Process
.
Techniques and Methods 6-A16
,
U.S. Geological Survey, Reston, VI
, p.
253
.
Hill
M. C.
&
Tiedman
C. R.
2007
Effective Groundwater Model Calibration
.
John Wiley & Sons Inc.
,
Hoboken, NJ
,
USA
, p.
384
.
Hill
M. C.
,
Banta
E. R.
,
Harbaugh
A. W.
&
Anderman
E. R.
2000
MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water Model – User Guide to the Observation, Sensitivity, and Parameter-Estimation Processes and Three Post-Processing Programs
.
U.S. Geological Survey Open-File Report 00-184
, p.
210
.
Kallergis
G.
2000
Applied Environmental Hydrogeology
,
2nd edn
.
Technical Chamber of Greece
,
Athens
,
Greece
, p.
345
.
Kallioras
A.
2008
Groundwater Resources Management of Aquifers Subjected to Seawater Intrusion Regime. The Case of Western Coastal Plain of the Prefecture of Rhodope
.
Doctoral dissertation
,
Department of Civil Engineering, Democritus University of Thrace
,
Xanthi
,
Greece
(in Greek)
.
Kallioras
A.
,
Pliakas
F.
,
Diamantis
I.
&
Emmanouil
M.
2006a
Application of Geographical Information Systems (GIS) for the management of coastal aquifers subjected to seawater intrusion
.
J. Environ. Sci. Health, Part A: Toxic/Hazard. Subst. Environ. Eng.
A41
(
9
),
2027
2044
.
Kallioras
A.
,
Pliakas
F.
,
Diamantis
I.
&
Kallergis
G.
2006c
Seawater intrusion and management of coastal aquifers in Greece. The case study of Rhodope western coastal aquifer
. In:
International Water Association, IWA, 10th International Conference on Diffuse Pollution
,
18–22 September 2006
,
Istanbul, Turkey
.
Kallioras
A.
,
Pliakas
F.
&
Diamantis
I.
2010
Simulation of groundwater flow in a sedimentary aquifer system subjected to overexploitation
.
Water Air Soil Pollut.
211
,
177
201
.
doi: 10.1007/s11270-009-0291-6
.
Kopač
I.
&
Vremec
M.
2017
Slovenian test case Vrbanski Plato aquifer in the EU HORIZON 2020 FREEWAT project
.
Acque Sotterranee – Ital. J. Groundwater
6
(
3/149
),
15
25
.
doi: 10.7343/as-2017-287
.
Panteleit
B.
,
Jensen
S.
,
Seiter
K.
&
Siebert
Y.
2017
Das Bremerhavener Grundwasser im Klimawandel - Eine FREEWAT-Fallstudie. (Bremerhaven's groundwater in climate change)
.
Grundwasser
23
(
3
),
233
244
(in German)
.
doi: 10.1007/s00767-017-0385-9
.
Perdikaki
M.
,
Pouliaris
C.
,
Borsi
I.
,
Rossetto
R.
&
Kallioras
A.
2017
Management of coastal hydrosystems through the application of free and open source software tool FREEWAT
.
European Water
57
,
383
388
.
Petalas
C.
1997
Analysis of Aquifer Systems in the Heterogeneous Coastal Part of Prefecture of Rhodope
.
Doctoral dissertation
,
Department of Civil Engineering, Democritus University of Thrace
,
Xanthi
,
Greece
.
Poeter
E. P.
,
Hill
M. C.
,
Banta
E. R.
,
Mehl
S.
&
Christensen
S.
2005
UCODE_2005 and Six Other Computer Codes for Universal Sensitivity Analysis, Calibration, and Uncertainty Evaluation
.
Techniques and Methods 6-A11
,
U.S. Geological Survey, Reston, VI
, p.
283
.
Poeter
E. P.
,
Hill
M. C.
,
Lu
D.
,
Tiedeman
C. R.
&
Mehl
S.
2014
UCODE_2014, with new Capabilities to Define Parameters Unique to Predictions, Calculate Weights Using Simulated Values, Estimate Parameters with SVD, Evaluate Uncertainty with MCMC, and More
.
Integrated Groundwater Modeling Center Report Number GWMI 2014-02
.
Positano
P.
&
Nannucci
M.
2017
The H2O20 FREEWAT participated approach for the Follonica-Scarlino aquifer case study. A common space to generate shared knowledge on the value of water
.
Acque Sotterranee – Ital. J. Groundwater
6
(
3/149
),
27
38
.
doi: 10.7343/as-2017-290
.
QGIS Development Team
2009
QGIS Geographic Information System
.
Open Source Geospatial Foundation Project
.
http://qgis.osgeo.org (accessed 26 February 2018)
.
Reilly
T.
&
Harbaugh
A.
2004
Guidelines for Evaluating Ground-Water Flow
.
Scientific Investigations Report 2004-5038. U.S. Department of Interior, U.S. Geological Survey
, p.
30
.
Rossetto
R.
,
De Filippis
G.
,
Borsi
I.
,
Foglia
L.
,
Cannata
M.
,
Criollo
R.
&
Vázquez-Suñé
E.
2018
Integrating free and open source tools and distributed modelling codes in GIS environment for data-based groundwater management
.
Environ. Model. Softw.
107
,
210
230
.
SpatialLite Development Team
2011
The Gaia-SINS Federated Projects Home-Page
.
http://www.gaia-gis.it/gaia-sins (accessed 20 February 2018)
.
Stamatakis
G.
1992
Investigation of Water Potential of Rhodope Prefecture
.
Undergraduate thesis
,
Civil Engineering Department, Democritus University of Thrace
,
Xanthi
,
Greece
.
Wang
L.
,
Jackson
C. R.
,
Pachocka
M.
&
Kingdon
A.
2016
A seamlessly coupled GIS and distributed groundwater flow model
.
Environ. Model. Softw.
82
,
1
6
.
Zheng
C.
&
Wang
P. P.
1999
MT3DMS, A Modular Three-Dimensional Multi-Species Transport Model for Simulation of Advection, Dispersion and Chemical Reactions of Contaminants in Groundwater Systems
.
U.S. Army Engineer Research and Development Center Contract Report SERDP-99-1
,
Vicksburg, MS
,
USA
, p.
202
.