1r.gwflow(1)                 GRASS GIS User's Manual                r.gwflow(1)
2
3
4

NAME

6       r.gwflow   -  Numerical calculation program for transient, confined and
7       unconfined groundwater flow in two dimensions.
8

KEYWORDS

10       raster, groundwater flow, hydrology
11

SYNOPSIS

13       r.gwflow
14       r.gwflow --help
15       r.gwflow [-f]  phead=name  status=name  hc_x=name  hc_y=name   [q=name]
16       s=name   [recharge=name]   top=name  bottom=name output=name  [vx=name]
17       [vy=name]        [budget=name]       type=string       [river_bed=name]
18       [river_head=name]           [river_leak=name]          [drain_bed=name]
19       [drain_leak=name]    dtime=float    [maxit=integer]     [maxit=integer]
20       [error=float]    [solver=name]    [--overwrite]   [--help]  [--verbose]
21       [--quiet]  [--ui]
22
23   Flags:
24       -f
25           Allocate a full quadratic linear  equation  system,  default  is  a
26           sparse linear equation system.
27
28       --overwrite
29           Allow output files to overwrite existing files
30
31       --help
32           Print usage summary
33
34       --verbose
35           Verbose module output
36
37       --quiet
38           Quiet module output
39
40       --ui
41           Force launching GUI dialog
42
43   Parameters:
44       phead=name [required]
45           Name of input raster map with initial piezometric head in [m]
46
47       status=name [required]
48           Name  of  input  raster  map  providing  boundary condition status:
49           0-inactive, 1-active, 2-dirichlet
50
51       hc_x=name [required]
52           Name of input raster map with x-part of the hydraulic  conductivity
53           tensor in [m/s]
54
55       hc_y=name [required]
56           Name  of input raster map with y-part of the hydraulic conductivity
57           tensor in [m/s]
58
59       q=name
60           Name of input raster map with water sources and sinks in [m^3/s]
61
62       s=name [required]
63           Name of input raster map with storativity for confined or effective
64           porosity for unconfined groundwater flow booth in [-]
65
66       recharge=name
67           Recharge input raster map e.g: 6*10^-9 per cell in [m^3/s*m^2]
68
69       top=name [required]
70           Name  of input raster map describing the top surface of the aquifer
71           in [m]
72
73       bottom=name [required]
74           Name of input raster map  describing  the  bottom  surface  of  the
75           aquifer in [m]
76
77       output=name [required]
78           Output raster map storing the numerical result [m]
79
80       vx=name
81           Output  raster  map to store the groundwater filter velocity vector
82           part in x direction [m/s]
83
84       vy=name
85           Output raster map to store the groundwater filter  velocity  vector
86           part in y direction [m/s]
87
88       budget=name
89           Output  raster  map  to  store the groundwater budget for each cell
90           [m^3/s]
91
92       type=string [required]
93           The type of groundwater flow
94           Options: confined, unconfined
95           Default: confined
96
97       river_bed=name
98           Name of input raster map providing the height of the river  bed  in
99           [m]
100
101       river_head=name
102           Name  of  input  raster map providing the water level (head) of the
103           river with leakage connection in [m]
104
105       river_leak=name
106           Name of input raster map providing the leakage coefficient  of  the
107           river bed in [1/s].
108
109       drain_bed=name
110           Name  of  input raster map providing the height of the drainage bed
111           in [m]
112
113       drain_leak=name
114           Name of input raster map providing the leakage coefficient  of  the
115           drainage bed in [1/s]
116
117       dtime=float [required]
118           The calculation time in seconds
119           Default: 86400
120
121       maxit=integer
122           Maximum  number of iteration used to solve the linear equation sys‐
123           tem
124           Default: 10000
125
126       maxit=integer
127           The maximum number of iterations in the linearization approach
128           Default: 25
129
130       error=float
131           Error break criteria for iterative solver
132           Default: 0.000001
133
134       solver=name
135           The type of solver which should solve the symmetric linear equation
136           system
137           Options: cg, pcg, cholesky
138           Default: cg
139

DESCRIPTION

141       This  numerical  program  calculates  implicit  transient, confined and
142       unconfined groundwater flow in two dimensions based on raster maps  and
143       the  current region settings.  All initial and boundary conditions must
144       be provided as raster maps. The unit in the location must be meters.
145
146       This module is sensitive to mask settings. All cells which are  outside
147       the mask are ignored and handled as no flow boundaries.
148
149       Workflow of r.gwflow
150
151
152       r.gwflow  calculates the piezometric head and optionally the water bud‐
153       get and the filter velocity field, based on the hydraulic  conductivity
154       and the piezometric head.  The vector components can be visualized with
155       paraview if they are exported with r.out.vtk.
156       The groundwater flow will always be calculated  transient.   For  stady
157       state  computation set the timestep to a large number (billions of sec‐
158       onds) or set the storativity/ effective porosity raster map to zero.
159       The water budget is calculated for each non inactive cell. The  sum  of
160       the  budget  for  each non inactive cell must be near zero.  This is an
161       indicator of the quality of the numerical result.
162

NOTES

164       The groundwater flow calculation is based on Darcy’s law and a  numeri‐
165       cal  implicit  finite volume discretization. The discretization results
166       in a symmetric and positive definite linear equation system in form  of
167       Ax = b, which must be solved. The groundwater flow partial differential
168       equation is of the following form:
169
170       (dh/dt)*S = div (K grad h) + q
171
172       In detail for 2 dimensions:
173
174       (dh/dt)*S = Kxx * (d^2h/dx^2) + Kyy * (d^2h/dy^2) + q
175
176           ·   h -- the piezometric head im [m]
177
178           ·   dt -- the time step for transient calculation in [s]
179
180           ·   S -- the specific storage [1/m]
181
182           ·   Kxx -- the hydraulic conductivity tensor part in x direction in
183               [m/s]
184
185           ·   Kyy -- the hydraulic conductivity tensor part in y direction in
186               [m/s]
187
188           ·   q - inner source/sink in meter per second [1/s]
189
190       Confined and unconfined groundwater flow is supported.  Be  aware  that
191       the  storativity  input  parameter  is  handled  differently in case of
192       unconfined flow. Instead of the storativity, the effective porosity  is
193       expected.
194
195       To  compute  unconfined  groundwater  flow,  a simple Picard based lin‐
196       earization scheme is used to solve the  resulting  non-linear  equation
197       system.
198
199       Two  different  boundary  conditions are implemented, the Dirichlet and
200       Neumann conditions. By default the calculation area  is  surrounded  by
201       homogeneous  Neumann boundary conditions.  The calculation and boundary
202       status of single cells must be set with a  status  map,  the  following
203       states are supportet:
204
205           ·   0  == inactive - the cell with status 0 will not be calculated,
206               active cells will have a no flow boundary to this cell
207
208           ·   1 == active - this cell is used for groundwater floaw  calcula‐
209               tion, inner sources and recharge can be defined for those cells
210
211           ·   2 == Dirichlet - cells of this type will have a fixed piezomet‐
212               ric head value which do not change over the time
213       Note that all required raster maps are read into main memory. Addition‐
214       ally  the  linear equation system will be allocated, so the memory con‐
215       sumption of this module rapidely grow with the size of the input maps.
216       The resulting linear equation system Ax = b can be solved with  several
217       solvers.   An iterative solvers with sparse and quadratic matrices sup‐
218       port is implemented.  The conjugate gradients  method  with  (pcg)  and
219       without  (cg)  precondition.   Additionally a direct Cholesky solver is
220       available. This direct solver only work with normal quadratic matrices,
221       so  be  careful  using  them with large maps (maps of size 10.000 cells
222       will need more than one gigabyte  of  RAM).   Always  prefer  a  sparse
223       matrix solver.
224

EXAMPLE

226       Use  this  small  script  to create a working groundwater flow area and
227       data. Make sure you are not  in  a  lat/lon  projection.   It  includes
228       drainage and river input as well.
229       # set the region accordingly
230       g.region res=25 res3=25 t=100 b=0 n=1000 s=0 w=0 e=1000 -p3
231       #now create the input raster maps for confined and unconfined aquifers
232       r.mapcalc expression="phead = if(row() == 1 , 50, 40)"
233       r.mapcalc expression="status = if(row() == 1 , 2, 1)"
234       r.mapcalc expression="well = if(row() == 20 && col() == 20 , -0.01, 0)"
235       r.mapcalc expression="hydcond = 0.00025"
236       r.mapcalc expression="recharge = 0"
237       r.mapcalc expression="top_conf = 20.0"
238       r.mapcalc expression="top_unconf = 70.0"
239       r.mapcalc expression="bottom = 0.0"
240       r.mapcalc expression="null = 0.0"
241       r.mapcalc expression="poros = 0.15"
242       r.mapcalc expression="s = 0.0001"
243       # The maps of the river
244       r.mapcalc expression="river_bed = if(col() == 35 , 48, null())"
245       r.mapcalc expression="river_head = if(col() == 35 , 49, null())"
246       r.mapcalc expression="river_leak = if(col() == 35 , 0.0001, null())"
247       # The maps of the drainage
248       r.mapcalc expression="drain_bed = if(col() == 5 , 48, null())"
249       r.mapcalc expression="drain_leak = if(col() == 5 , 0.01, null())"
250       #confined groundwater flow with cg solver and sparse matrix, river and drain
251       #do not work with this confined aquifer (top == 20m)
252       r.gwflow solver=cg top=top_conf bottom=bottom phead=phead status=status \
253         hc_x=hydcond hc_y=hydcond q=well s=s recharge=recharge output=gwresult_conf \
254         dt=8640000 type=confined vx=gwresult_conf_velocity_x vy=gwresult_conf_velocity_y budget=budget_conf
255       #unconfined groundwater flow with cg solver and sparse matrix, river and drain are enabled
256       # We use the effective porosity as storativity parameter
257       r.gwflow solver=cg top=top_unconf bottom=bottom phead=phead \
258         status=status hc_x=hydcond hc_y=hydcond q=well s=poros recharge=recharge \
259         river_bed=river_bed river_head=river_head river_leak=river_leak \
260         drain_bed=drain_bed drain_leak=drain_leak \
261         output=gwresult_unconf dt=8640000 type=unconfined vx=gwresult_unconf_velocity_x \
262         budget=budget_unconf vy=gwresult_unconf_velocity_y
263       # The data can be visulaized with paraview when exported with r.out.vtk
264       r.out.vtk -p in=gwresult_conf,status vector=gwresult_conf_velocity_x,gwresult_conf_velocity_y,null \
265         out=/tmp/gwdata_conf2d.vtk
266       r.out.vtk -p elevation=gwresult_unconf in=gwresult_unconf,status vector=gwresult_unconf_velocity_x,gwresult_unconf_velocity_y,null \
267         out=/tmp/gwdata_unconf2d.vtk
268       #now load the data into paraview
269       paraview --data=/tmp/gwdata_conf2d.vtk &
270       paraview --data=/tmp/gwdata_unconf2d.vtk &
271

SEE ALSO

273        r.solute.transport, r3.gwflow, r.out.vtk
274

AUTHOR

276       Sören Gebbert
277
278       This  work  is  based on the Diploma Thesis of Sören Gebbert available
279       here at Technical University Berlin in Germany.
280

SOURCE CODE

282       Available at: r.gwflow source code (history)
283
284       Main index | Raster index | Topics index | Keywords index  |  Graphical
285       index | Full index
286
287       © 2003-2020 GRASS Development Team, GRASS GIS 7.8.5 Reference Manual
288
289
290
291GRASS 7.8.5                                                        r.gwflow(1)
Impressum