Tutorial 0 – Part B

Modelling groundwater flow using

FREEWAT
v.2.9 released on 20/09/2019

Giovanna De Filippis1, Iacopo Borsi2, Matteo Ghetta1, Rudy Rossetto1

1. Scuola Superiore Sant’Anna – Pisa (Italy)

2. TEA SISTEMI S.p.A. – Pisa (Italy)
 DOCUMENTATION LICENSES

Please attribute FREEWAT

with a link to http://www.freewat.eu/

Except where otherwise noted,

this slides are licensed under a
Creative Commons Attribution 4.0 International License:
http://creativecommons.org/licenses/by/4.0/

FREEWAT - Free and Open Source Software Tools for Water Resource Management
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant
agreement No 642224
 DATA LICENSES

Please attribute FREEWAT with a link to

http://www.freewat.eu/

Except where otherwise noted,

data used in this tutorial are licensed under a
Open Database License:
http://opendatacommons.org/licenses/odbl/1.0/

Any rights in individual contents of the databases are

licensed under the Database Contents License:
http://opendatacommons.org/licenses/dbcl/1.0/
FREEWAT - Free and Open Source Software Tools for Water Resource Management
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant
agreement No 642224
 Introduction

Before starting with this second part of the tutorial, please recall that its
general aim is to introduce some capabilities of the GIS-integrated
FREEWAT platform, by developing a numerical model inspired to a
synthetic case study.

To include also the geographical component, the synthetic model was

located in the context of a natural alluvial plain north of Grosseto (Italy)
and crossed by the Bruna river along a roughly straight line from north
southward (Fig.1 in Tutorial 0 - Part A).

The estimated time required to perform this part of the tutorial is about
three hours.

1
 Previous steps/1

In the first part of the tutorial (Part A) the model was run over two Stress
Periods (SPs) to achieve the following specific objectives.

 The simulation of the natural groundwater flow in steady-state

condition and the investigation of the relationships between ground-
and surface-water in terms of water budget.

 The assessment and quantification of the hydraulic head decrease in

each model layer, due to the presence of two pumping wells
penetrating the deeper aquifer, still focusing on the ground-/surface-
water exchanges.

2
 Previous steps/2

The two steady-state SPs were defined as following:

 The 1st SP lasted one day, during which no stresses (recharge,

pumping wells, etc.) affected the hydrogeological system. At the end
of this SP, the hydraulic heads in each model layer were computed
and represented as color and contour maps (Fig.11 in Tutorial 0 - Part
A). The groundwater balance (Fig.12 in Tutorial 0 - Part A) was then
read to check consistency with the conceptual model (Fig.3 and Fig.4
in Tutorial 0 - Part A).

 The 2nd SP lasted 30 days, during which two pumping wells

penetrating the deepest aquifer were activated. At the end of this SP,
the same analysis of the model outputs was carried on (Fig.15 and
Fig.16 in Tutorial 0 - Part A).

3
 Previous steps/3

 A raster analysis was finally performed to quantify the decrease of

the hydraulic head in each model layer, due to the presence of the
pumping wells (Fig.17 in Tutorial 0 - Part A).

While building the model, some MODFLOW Packages were activated to

simulate specific boundary conditions and source/sink terms.

Several QGIS tools were used as well to handle and represent data and
features. In this regard, coupling the use of relevant free and open
source codes for hydrological modelling and the potentialities of QGIS
makes the FREEWAT platform powerful and easily accessible.

4
 Objectives/1
This second part of the tutorial has three specific objectives.

 The numerical model will be further modified to simulate the impact

of a Managed Aquifer Recharge scheme (an Infiltration Pond, IP) at
the center of the domain to recharge at 2*10-7 m/s the superficial
aquifer. This analysis will be performed without pumping.

AUSTRALIAN GUIDELINES
FOR WATER RECYCLING

 The effects of such infiltration pond will be evaluated by computing,

for each model layer, the difference between the hydraulic head
simulated under stressed conditions (with managed recharge
applied) and that simulated in steady-state conditions (with no
recharge applied).
AUSTRALIAN GUIDELINES FOR WATER RECYCLING: MANAGING HEALTH AND ENVIRONMENTAL RISKS (PHASE 2). Managed Aquifer Recharge (Natural Resource
Management Ministerial Council + Environment Protection and Heritage Council + National Health and Medical Research Council, 2009). 5
 Objectives/2
 Then, after activating pumping, we will evaluate the impact
of the Managed Aquifer Recharge scheme on the whole
water budget.
One more technical objective is the following.
 Use of the MODFLOW Recharge (RCH) Package to simulate areally-
distributed recharge to the groundwater system.

ASSUMPTIONS!!!
There are a couple of assumptions underlying this tutorial.
1) During the simulations we will not be changing the specified-head values at the
western boundary. However, by means of the IP we will create a groundwater mound.
In the real world, this would change the head at the western boundary. We will then
assume that no change in head is caused by the presence of such IP.
2) We set two no-flow boundaries at the northern and southern boundaries of the model
domain for all the simulations. In reality, this will not hold true and in fact it will prevent
flow from the groundwater mound to move out through these two boundaries.
Because of this, the analysis on head distributions and of the water budget have to be
considered cautiously. We set here these considerations just to show the kind of analysis
that can be performed.
6
In a next release, the tutorial will be expanded in space in order to skip these issues.
 Time discretisation
Taking step from the previously implemented model, we will add
and run a 3rd SP, deactivating the abstraction wells and activating
the IP.

Then, a 4th SP will be run to represent a further configuration of the

model, with both the wells and the IP activated.

The following scheme is so adopted for the entire tutorial:

Stress Period From (sec) To (sec) Length State Time step Stresses
(days) involved
1* 0 86400 1 Steady-state 1 None

2* 86400 2678400 30 Steady-state 1 Wells

3 2678400 5270400 30 Steady-state 1 Recharge

4 5270400 7862400 30 Steady-state 1 Recharge +
Wells

* These SPs were already performed in Tutorial 0 - Part A

7
 Preliminary steps/1

• Run QGIS and open the project ‘t0.qgs ’

• Just to recall what is the state of your model, use the following menu
to check which Model Data Objects (MDOs) are stored in the Database
(DB):
Database -> DB Manager -> DB Manager
- Expand SpatiaLite
- Expand t0.sqlite

• The following MDOs should be stored within the DB ‘t0.sqlite’ :

- ‘grid_grid’ (the model grid)
- ‘upper_aquifer’ (model layer 1)
- ‘aquitard’ (model layer 2)
- ‘lower_aquifer’ (model layer 3)

8
 Preliminary steps/2

- ‘river_1sp_riv’ (RIV boundary condition activated during the

1st SP) and the related table ‘river_1sp_riv_table’
- ‘river_2sp_riv’ (RIV boundary condition activated during the
2nd SP) and the related table ‘river_2sp_riv_table’
- ‘pumping_2sp_well’ (WEL sink term activated during the 2nd SP)

• All these MDOs should also be listed in the Layers Panel (LP) and
(some of them) eventually displayed in the Map canvas.

9
 The infiltration pond
We will now simulate the presence of an infiltration pond at the
center of the domain (Fig.18) as Managed Aquifer Recharge
scheme infiltrating 12 l/s:

Fig.18 Location of the infiltration pond. GH is the profile along which

the hydrostratigraphic section shown in Fig.19 is drawn. 10
 Add stress period/1

The infiltration pond will be active for a period of 30 days. So we will

now define a 3rd SP 30 days long.
While defining the SP length, remember that time unit of this model is
seconds!

• To define the 3rd SP use the menu:

FREEWAT -> Model Setup -> Add Stress Period
- Time Table : timetable_t0
- SP Length : 2592000
- Time Steps : 1
- Multiplier : 1.0
- State: Steady State
- Click OK

11
 Add stress period/2

You will see that the table ‘timetable_t0 ’ is updated and contains three
records with information about the SPs defined so far.

12
 Boundary conditions and source term
We will now specify boundary conditions and the source term
defined during SP 3 (Fig.19):

- Specified-head along the western boundary

- No-flow at the northern and southern boundaries and at the
bottom of the system (no particular package/setting is required)
- River along the eastern boundary (MODFLOW RIV Package is used)
- Recharge to the uppermost model layer, through the infiltration pond
(MODFLOW RCH Package is used)
Specified head
(m msl)

Fig.19 Representation of boundary conditions and the source term during SP 3 (vertical
cross-section along the profile GH; modified after McDonald et al., 1992). 13
 MODFLOW RCH Package

MODFLOW Recharge (RCH) Package allows to simulate areally

distributed, direct recharge to groundwater, usually used to
represent rainfall recharge.

The User must define, for each SP, the recharge flux to be applied to the
map area, in units of length per time [L/T]. This recharge flux is then
internally multiplied by the area of each cell, to get the recharge flow rate
at each cell, then expressed as a fluid volume per unit time [L3/T].

The User must also specify if the areal recharge has to be applied (irch
parameter):
(1) to the upper model layer,
(2) to the uppermost variable-head cell in each vertical column, or
(3) to the model layer specified by irch .

14
 RCH source term/1

• First of all, define the RCH source term by activating the RCH
Package with the menu:
FREEWAT -> MODFLOW Boundary Conditions -> Create RCH Layer

- Model Name : t0
- Grid Layer : grid_grid
- Name of new layer : recharge_3sp
- Click OK

A new MDO is created (‘recharge_3sp_rch ’), loaded in the LP and

displayed in the Map canvas as a set of grid cells.
One rate value is assigned at the center of all grid cells and for each SP by
default.
15
 RCH source term/2

The Attribute Table of this MDO contains several records (one for each
grid cell) and data for each SP:

- sp_1_rech : recharge rate values assigned at each grid cell during

the 1st SP;
- sp_1_irch. This field refers to the layer number to which recharge
is applied during the 1st SP. For the aims of this tutorial
its value is not important, because recharge will be applied
to the top layer and this will be stated during the Run phase;
- sp_2_rech and sp_3_rech , sp_2_irch and sp_3_irch : defined as for SP 1.

16
 RCH source term/3
• Open the Attribute Table of the MDO ‘recharge_3sp_rch ’,
activate the editing mode and use the expression bar to set:
- sp_1_rech = 0 m/s
- sp_2_rech = 0 m/s
- sp_3_rech = 0 m/s
- Each time click on Update All (this means that a null recharge flux
is applied at each grid cell and for each SP)

• Save your edits and deactivate the editing mode

17
 RCH source term/4
• To apply the recharge just in the area where the infiltration
pond is located, by using the MODFLOW RCH Package, load
the shapefile inf_pond.shp from the folder
input_data/infiltration_pond.
The area where the recharge will be applied is displayed in the Map
canvas:

View of the Map area

For the aims of the numerical model, the cells involved in the application
of the MODFLOW RCH Package will be those intersected by this polygon. 18
 RCH source term/5
• You must first select grid cells intersected by this polygon.
To do so, select the ‘recharge_3sp_rch ’ layer in the LP and use
Select Features by Polygon [see 1.]. You must draw a polygon similar to
that contained in the inf_pond.shp file (single click on each vertex and
right-click when finished)
• You should see the selected (yellow) cells in the Map canvas. If not,
make sure the first visible layer is ‘recharge_3sp_rch ’, by unchecking all
the other overlaid layers in the LP

1.

19
 RCH source term/6
• Open the Attribute Table of the MDO ‘recharge_3sp_rch ’
• Select Show Selected Features (just 6 blue records should
be displayed) to apply the recharge flux just to the selected cells
• Then, in editing mode, select the field sp_3_rech and:
- In the expression bar, type 2e-7 m/s (such recharge flux is applied
at each of these six cells during SP 3)
- Click Update Selected
- Click Save Edits
- Deactivate the Toggle editing mode
- Close the Attribute Table and deselect all cells with

20
 Update RIV boundary condition
• Before running the model
you must update the RIV
boundary condition, defining
river parameters for the 3rd SP.
To do so, you must create a
new MDO (‘river_3sp_riv ’),
exactly as done before, but
using the file riv_par_3sp.csv
in the folder
input_data/river_parameters

Please notice that specified-head and no-flow boundary conditions will not be
updated: while the no-flow boundary is a default setting in MODFLOW, when non-
boundary is specified, the specified-head was defined by editing the fields ACTIVE
and STRT in the Attribute Table of each model layer. 21
 Update WEL sink term
• Before running the model you must update also the WEL sink
term. To do so, select the ‘grid_grid ’ layer in the LP and use
Select features by area or single click to select cells (5;6) and (10;5)
• Then create a new MDO (‘pumping_3sp_well ’) exactly as done in
Tutorial 0 - Part A and edit the Attribute Table of this new MDO as in
the figure below. Remember that the wells are active just during SP 2!

22
 Run the model/1

• Now we can run the model. In the Run Model window you must
check the RIV Package (and the corresponding MDO created
‘river_3sp_riv ’), the WEL Package (and the corresponding MDO created
‘pumping_3sp_well’ ) and the RCH Package (and the corresponding
MDO created ‘recharge_3sp_rch ’)

• Among the Rch Options, select ‘Recharge to top grid ’ in the drop-down
menu to apply recharge to model layer 1.
As already stated before, other two Rch Options are available:
- ‘Recharge layer defined in irch ’ (if you want to apply recharge to
selected model layers in selected cells, you must modify irch in
the Attribute Table of the rch MDO, as needed)
- ‘Recharge to highest active cell ’ (if you want to apply recharge to the
uppermost active cell, which is not necessarily in model layer 1)
23
Run the model/2

24
 View model output

• Now, we will display model outputs simulated at the end of the

3rd SP, edit their view style and extract contours from these raster data

25
Head and contour at the end of SP 3

Model layer 1 Model layer 2

Model layer 3
Fig.20 Hydraulic head simulated at the end of SP 3, represented
as color maps and contour lines (values in meters msl). 26
 MODFLOW input and output files
MODFLOW input and output files are updated and the input file
t0.rch is created, as you can also verify within the t0.nam file.

As done before, you can check the model balance at the end of SP 3 in
the t0.list file (Fig.21).

Fig.21 Model balance report at the end of SP 3 (content of the t0.list file). 27
 Check model balance
With reference to Fig.21, and with respect to Fig.12 in Tutorial 0 –
Part A, a new flow term is now involved in the model balance at the
end of the 3rd SP: the source term due to the presence of the IP
(RECHARGE - IN). This term is worth 1.20*10-2 m3/s (1036.80 m3/day).

Due to the head increase at the location of the IP, both inflow and outflow
occur through the western boundary (CONSTANT HEAD – IN and
CONSTANT HEAD – OUT).

The discharge at the eastern boundary (RIVER LEAKAGE - OUT) is

estimated to be about 1.62*10-2 m3/s (1399.68 m3/day). This last term is
higher than at the end of the 1st SP (it was 1.13*10-2 m3/s - 976.32 m3/day;
see Fig.12 in Tutorial 0 – Part A).
The increase of the hydraulic head where the IP is located is generating
flow which, in turn, is drained by the river.

28
 Raster analysis/1

We will now evaluate such head increase due to recharge to the top layer.

In order to get this result, use the Raster Calculator tool to subtract
the head calculated for each model layer at the end SP 3 minus the
head calculated for each model layer at the end of SP 1.

29
 Raster analysis/2

Model layer 1 Model layer 2

Model layer 3
Fig.22 Difference between hydraulic heads simulated at the end of SPs 1 and 3,
represented as color maps and contour lines (values in meters). 30
 Add stress period/1
We will now simulate the presence of the infiltration pond and the
two pumping wells at the same time (Fig.23).

Fig.23 Location of the pumping wells and infiltration pond. 31

 Add stress period/2

The infiltration pond and pumping wells will be active at the same
time for a period of 30 days. So we will define one more SP 30 days long.
While defining the SP length, remember that time unit of this model is
seconds!

• To define the 4th SP use the menu:

FREEWAT -> Model Setup -> Add Stress Period
- Time Table : timetable_t0
- SP Length : 2592000
- Time Steps : 1
- Multiplier : 1.0
- State: Steady State
- Click OK

32
 Add stress period/3

You will see that the table ‘timetable_t0 ’ is updated and contains four
records with information about the SPs defined so far.

33
 Boundary conditions and
source/sink terms
We have now to specify boundary conditions and source/sink
terms defined during SP 4:

- Specified-head * along the western boundary

- No-flow * at the northern and southern boundaries and at the
bottom of the system (no particular package/setting is required)
- River along the eastern boundary (MODFLOW RIV Package is used)
- Recharge to the uppermost model layer, through the infiltration pond
(MODFLOW RCH Package is used)
- Pumping wells in the lower layer (MODFLOW WEL Package is used)

* Please notice that specified-head and no-flow boundary conditions will not be
updated: while the no-flow boundary is a default setting in MODFLOW, when
non-boundary is specified, the specified-head was defined by editing the fields
ACTIVE and STRT in the Attribute Table of each model layer. 34
 Update RIV boundary conditions

• Before running the model

you must update the RIV
boundary condition, defining
river parameters for the 4th SP.
To do so, you must create a
new MDO (‘river_4sp_riv ’),
exactly as done before, but
using the file riv_par_4sp.csv
in the folder
input_data/river_parameters

35
 Update WEL sink term
• Before running the model you must update also the WEL sink
term. To do so, select the ‘grid_grid ’ layer in the LP and use
Select features by area or single click to select cells (5;6) and (10;5)
• Then create a new MDO (‘pumping_4sp_well ’) exactly as done before
and edit the Attribute Table of this new MDO as in the figure below.
Remember that the wells are active just during SPs 2 and 4!

36
 Update RCH source term
• Before running the model you must update also the RCH recharge
term. To do so, create a new MDO (‘recharge_4sp_rch ’) exactly as
done before and set to zero the recharge flux at all cells and for each
SP
• Then select the cells intersected by the infiltration pond, by using
Select features by polygon , edit the selected records in the
Attribute Table of this new MDO as
in the figure and deselect all cells.
Remember that the recharge is
active just during SPs 3 and 4!

37
 Run the model/1

• Now you can run the model. In the Run Model window you must check
the RIV Package (and the corresponding MDO created ‘river_4sp_riv ’),
the WEL Package (and the corresponding MDO created
‘pumping_4sp_well ’) and the RCH Package (and the corresponding
MDO created ‘recharge_4sp_rch ’).

• Among the Rch Options, select ‘Recharge to top grid ’ in the drop-
down menu to apply recharge to model layer 1.

38
Run the model/2

39
 View model output

• Now, we will display model outputs simulated at the end of the

4th SP, edit their view style and extract contours from these raster data

40
Head and contour at the end of SP 4

Model layer 1 Model layer 2

Model layer 3
Fig.24 Hydraulic head simulated at the end of SP 4, represented
as color maps and contour lines (values in meters msl). 41
 MODFLOW input and output files

As done before, you can check the model balance at the end of SP
4 in the t0.list file (Fig.25).

Fig.25 Model balance report at the end of SP 4 (content of the t0.list file).
42
 Check model balance

With reference to Fig.25, and with respect to Fig.12 in Tutorial 0 –

Part A, both RECHARGE – IN and WELLS – OUT terms are involved in
the groundwater balance at the end of SP 4: they are worth 1.20*10-2
m3/s (1036.80 m3/day) and 8.00*10-3 m3/s (691.20 m3/day), respectively.

Both inflow and outflow occur through the western boundary (CONSTANT
HEAD – IN and CONSTANT HEAD – OUT).

The discharge at the eastern boundary (RIVER LEAKAGE - OUT) is

estimated to be about 1.51*10-2 m3/s (1304.64 m3/day), which is lower
with respect to what estimated at the end of the 3rd SP (it was 1.62*10-2
m3/s - 1399.68 m3/day), due to the decrease of the hydraulic head where
the pumping wells are located.

43
 Raster analysis/1

We will now evaluate the impact of the infiltration pond in contrasting the
head decrease at the location of pumping wells.

In order to get this result, use the Raster Calculator tool to subtract the
head calculated for each model layer at the end SP 4 minus the head
calculated for each model layer at the end of SP 2.

You can use the Identity Features tool to query the rasters at the pumping
wells location.

44
 Raster analysis/2

• Set the Transparency of the three rasters obtained to 30% and make
sure the vector layer ‘wells’ in the LP is checked, so that you can easily
identify the pumping wells location

• Now, select the raster you want to query, click on Identity Features
and then click within the cells where the wells are located, one by
one

• The Identify Results Panel will appear below the LP, where the value
of the head difference between SPs 2 and 4 is displayed. See Fig.26 for
details.

45
 Raster analysis/3
head difference ≈ 1.86 m head difference ≈ 1.84 m

head difference ≈ 1.37 m head difference ≈ 1.36 m

Model layer 1 Model layer 2

head difference ≈ 1.82 m

head difference ≈ 1.35 m

Model layer 3
Fig.26 Difference between hydraulic heads simulated at the end of SPs 4 and 2,
represented as color maps and contour lines (values in meters). 46
 If you need any assistance, please contact

Giovanna De Filippis – Scuola Superiore Sant’Anna (Pisa - Italy)

g.defilippis@santannapisa.it

FREEWAT development received funding from the following projects:

1. Hydrological part was developed starting from the project SID&GRID, funded by Regione Toscana through EU POR-FSE 2007-2013
(sidgrid.isti.cnr.it) (2010-2013)
2. Porting of SID&GRID under QGIS has been performed through funds provided by Regione Toscana to Scuola Superiore S.Anna -
Project Evoluzione del sistema open source SID&GRID di elaborazione dei dati geografici vettoriali e raster per il porting negli
ambienti QGis e Spatialite in uso presso la Regione Toscana (CIG: ZA50E4058A) (2015)
3. Saturated zone solute transport simulation capability has been developed within the EU FP7-ENV-2013-WATER-INNO-DEMO
MARSOL. MARSOL project received funding from the European Union's Seventh Framework Programme for Research, Technological
Development and Demonstration under grant agreement n. 619120 (www.marsol.eu) (2014-2017)
4. FREEWAT was developed within the EU H2020 project FREEWAT - Free and Open Source Software Tools for Water Resource
Management. FREEWAT project received funding from the European Union’s Horizon 2020 research and innovation programme
under grant agreement n. 642224 (www.freewat.eu) (2015-2017)
5. Integration of SFT (StreamFlow Transport) and LKT (Lake Transport) packages of MT3D-USGS is being performed at Scuola
Superiore Sant'Anna within the project SMAQua (SMart ICT tools per l'utilizzo efficiente dell'AcQua) - co-financed by Regione Toscana,
ASA S.p.A. and ERM Italia S.p.A. (2018-2020)