1r.gwflow(1) GRASS GIS User's Manual r.gwflow(1)
2
6 r.gwflow - Numerical calculation program for transient, confined and
7 unconfined groundwater flow in two dimensions.
8
10 raster, groundwater flow, hydrology
11
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 Flags:
24 -f
25 Allocate a full quadratic linear equation system, default is a
26 sparse linear equation system.
27 --overwrite
29 Allow output files to overwrite existing files
30 --help
32 Print usage summary
33 --verbose
35 Verbose module output
36 --quiet
38 Quiet module output
39 --ui
41 Force launching GUI dialog
42 Parameters:
44 phead=name [required]
45 Name of input raster map with initial piezometric head in [m]
46 status=name [required]
48 Name of input raster map providing boundary condition status:
49 0-inactive, 1-active, 2-dirichlet
50 hc_x=name [required]
52 Name of input raster map with x-part of the hydraulic conductivity
53 tensor in [m/s]
54 hc_y=name [required]
56 Name of input raster map with y-part of the hydraulic conductivity
57 tensor in [m/s]
58 q=name
60 Name of input raster map with water sources and sinks in [m^3/s]
61 s=name [required]
63 Name of input raster map with storativity for confined or effective
64 porosity for unconfined groundwater flow booth in [-]
65 recharge=name
67 Recharge input raster map e.g: 6*10^-9 per cell in [m^3/s*m^2]
68 top=name [required]
70 Name of input raster map describing the top surface of the aquifer
71 in [m]
72 bottom=name [required]
74 Name of input raster map describing the bottom surface of the
75 aquifer in [m]
76 output=name [required]
78 Output raster map storing the numerical result [m]
79 vx=name
81 Output raster map to store the groundwater filter velocity vector
82 part in x direction [m/s]
83 vy=name
85 Output raster map to store the groundwater filter velocity vector
86 part in y direction [m/s]
87 budget=name
89 Output raster map to store the groundwater budget for each cell
90 [m^3/s]
91 type=string [required]
93 The type of groundwater flow
94 Options: confined, unconfined
95 Default: confined
96 river_bed=name
98 Name of input raster map providing the height of the river bed in
99 [m]
100 river_head=name
102 Name of input raster map providing the water level (head) of the
103 river with leakage connection in [m]
104 river_leak=name
106 Name of input raster map providing the leakage coefficient of the
107 river bed in [1/s].
108 drain_bed=name
110 Name of input raster map providing the height of the drainage bed
111 in [m]
112 drain_leak=name
114 Name of input raster map providing the leakage coefficient of the
115 drainage bed in [1/s]
116 dtime=float [required]
118 The calculation time in seconds
119 Default: 86400
120 maxit=integer
122 Maximum number of iteration used to solve the linear equation sys‐
123 tem
124 Default: 10000
125 maxit=integer
127 The maximum number of iterations in the linearization approach
128 Default: 25
129 error=float
131 Error break criteria for iterative solver
132 Default: 0.000001
133 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
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 This module is sensitive to mask settings. All cells which are outside
147 the mask are ignored and handled as no flow boundaries.
148 Workflow of r.gwflow
150 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
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 (dh/dt)*S = div (K grad h) + q
171 In detail for 2 dimensions:
173 (dh/dt)*S = Kxx * (d^2h/dx^2) + Kyy * (d^2h/dy^2) + q
175 · h -- the piezometric head im [m]
177 · dt -- the time step for transient calculation in [s]
179 · S -- the specific storage [1/m]
181 · Kxx -- the hydraulic conductivity tensor part in x direction in
183 [m/s]
184 · Kyy -- the hydraulic conductivity tensor part in y direction in
186 [m/s]
187 · q - inner source/sink in meter per second [1/s]
189 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 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 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 · 0 == inactive - the cell with status 0 will not be calculated,
206 active cells will have a no flow boundary to this cell
207 · 1 == active - this cell is used for groundwater floaw calcula‐
209 tion, inner sources and recharge can be defined for those cells
210 · 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
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
273 r.solute.transport, r3.gwflow, r.out.vtk
274
276 Sören Gebbert
277 This work is based on the Diploma Thesis of Sören Gebbert available
279 here at Technical University Berlin in Germany.
280
282 Available at: r.gwflow source code (history)
283 Main index | Raster index | Topics index | Keywords index | Graphical
285 index | Full index
286 © 2003-2020 GRASS Development Team, GRASS GIS 7.8.5 Reference Manual
288GRASS 7.8.5 r.gwflow(1)