1r.solute.transport(1) GRASS GIS User's Manual r.solute.transport(1)
2
6 r.solute.transport - Numerical calculation program for transient, con‐
7 fined and unconfined solute transport in two dimensions
8
10 raster, hydrology, solute transport
11
13 r.solute.transport
14 r.solute.transport --help
15 r.solute.transport [-fc] c=name phead=name hc_x=name hc_y=name sta‐
16 tus=name diff_x=name diff_y=name [q=name] [cin=name] cs=name
17 rd=name nf=name top=name bottom=name output=name [vx=name] [vy=name]
18 dtime=float [maxit=integer] [error=float] [solver=name] [re‐
19 lax=float] [al=float] [at=float] [loops=float] [stab=string]
20 [--overwrite] [--help] [--verbose] [--quiet] [--ui]
21 Flags:
23 -f
24 Use a full filled quadratic linear equation system, default is a
25 sparse linear equation system.
26 -c
28 Use the Courant-Friedrichs-Lewy criteria for time step calculation
29 --overwrite
31 Allow output files to overwrite existing files
32 --help
34 Print usage summary
35 --verbose
37 Verbose module output
38 --quiet
40 Quiet module output
41 --ui
43 Force launching GUI dialog
44 Parameters:
46 c=name [required]
47 The initial concentration in [kg/m^3]
48 phead=name [required]
50 The piezometric head in [m]
51 hc_x=name [required]
53 The x-part of the hydraulic conductivity tensor in [m/s]
54 hc_y=name [required]
56 The y-part of the hydraulic conductivity tensor in [m/s]
57 status=name [required]
59 The status for each cell, = 0 - inactive cell, 1 - active cell, 2 -
60 dirichlet- and 3 - transfer boundary condition
61 diff_x=name [required]
63 The x-part of the diffusion tensor in [m^2/s]
64 diff_y=name [required]
66 The y-part of the diffusion tensor in [m^2/s]
67 q=name
69 Groundwater sources and sinks in [m^3/s]
70 cin=name
72 Concentration sources and sinks bounded to a water source or sink
73 in [kg/s]
74 cs=name [required]
76 Concentration of inner sources and inner sinks in [kg/s] (i.e. a
77 chemical reaction)
78 rd=name [required]
80 Retardation factor [-]
81 nf=name [required]
83 Effective porosity [-]
84 top=name [required]
86 Top surface of the aquifer in [m]
87 bottom=name [required]
89 Bottom surface of the aquifer in [m]
90 output=name [required]
92 The resulting concentration of the numerical solute transport cal‐
93 culation will be written to this map. [kg/m^3]
94 vx=name
96 Calculate and store the groundwater filter velocity vector part in
97 x direction [m/s]
98 vy=name
100 Calculate and store the groundwater filter velocity vector part in
101 y direction [m/s]
102 dtime=float [required]
104 The calculation time in seconds
105 Default: 86400
106 maxit=integer
108 Maximum number of iteration used to solve the linear equation sys‐
109 tem
110 Default: 10000
111 error=float
113 Error break criteria for iterative solver
114 Default: 0.000001
115 solver=name
117 The type of solver which should solve the linear equation system
118 Options: gauss, lu, jacobi, sor, bicgstab
119 Default: bicgstab
120 relax=float
122 The relaxation parameter used by the jacobi and sor solver for
123 speedup or stabilizing
124 Default: 1
125 al=float
127 The longditudinal dispersivity length. [m]
128 Default: 0.0
129 at=float
131 The transversal dispersivity length. [m]
132 Default: 0.0
133 loops=float
135 Use this number of time loops if the CFL flag is off. The timestep
136 will become dt/loops.
137 Default: 1
138 stab=string
140 Set the flow stabilizing scheme (full or exponential upwinding).
141 Options: full, exp
142 Default: full
143
145 This numerical program calculates numerical implicit transient and
146 steady state solute transport in porous media in the saturated zone of
147 an aquifer. The computation is based on raster maps and the current re‐
148 gion settings. All initial- and boundary-conditions must be provided as
149 raster maps. The unit in the location must be meters.
150 This module is sensitive to mask settings. All cells which are outside
152 the mask are ignored and handled as no flow boundaries.
153 This module calculates the concentration of the solution and optional
154 the velocity field, based on the hydraulic conductivity, the effective
155 porosity and the initial piecometric heads. The vector components can
156 be visualized with paraview if they are exported with r.out.vtk.
157 Use r.gwflow to compute the piezometric heights of the aquifer. The
158 piezometric heights and the hydraulic conductivity are used to compute
159 the flow direction and the mean velocity of the groundwater. This is
160 the base of the solute transport computation.
161 The solute transport will always be calculated transient. For stady
162 state computation set the timestep to a large number (billions of sec‐
163 onds).
164 To reduce the numerical dispersion, which is a consequence of the con‐
165 vection term and the finite volume discretization, you can use small
166 time steps and choose between full and exponential upwinding.
167
169 The solute transport calculation is based on a diffusion/convection
170 partial differential equation and a numerical implicit finite volume
171 discretization. Specific for this kind of differential equation is the
172 combination of a diffusion/dispersion term and a convection term. The
173 discretization results in an unsymmetric linear equation system in form
174 of Ax = b, which must be solved. The solute transport partial differen‐
175 tial equation is of the following form:
176 (dc/dt)*R = div ( D grad c - uc) + cs -q/nf(c - c_in)
178 • c -- the concentration [kg/m^3]
180 • u -- vector of mean groundwater flow velocity
182 • dt -- the time step for transient calculation in seconds [s]
184 • R -- the linear retardation coefficient [-]
186 • D -- the diffusion and dispersion tensor [m^2/s]
188 • cs -- inner concentration sources/sinks [kg/m^3]
190 • c_in -- the solute concentration of influent water [kg/m^3]
192 • q -- inner well sources/sinks [m^3/s]
194 • nf -- the effective porosity [-]
196 Three different boundary conditions are implemented, the Dirichlet,
197 Transmission and Neumann conditions. The calculation and boundary sta‐
198 tus of single cells can be set with the status map. The following
199 states are supportet:
200 • 0 == inactive - the cell with status 0 will not be calculated,
202 active cells will have a no flow boundary to an inactive cell
203 • 1 == active - this cell is used for sloute transport calcula‐
205 tion, inner sources can be defined for those cells
206 • 2 == Dirichlet - cells of this type will have a fixed concen‐
208 tration value which do not change over time
209 • 3 == Transmission - cells of this type should be placed on
211 out-flow boundaries to assure the flow of the solute stream out
212 Note that all required raster maps are read into main memory. Addition‐
213 ally the linear equation system will be allocated, so the memory con‐
214 sumption of this module rapidely grow with the size of the input maps.
215 The resulting linear equation system Ax = b can be solved with several
216 solvers. Several iterative solvers with unsymmetric sparse and qua‐
217 dratic matrices support are implemented. The jacobi method, the
218 Gauss-Seidel method and the biconjugate gradients-stabilized (bicgstab)
219 method. Additionally a direct Gauss solver and LU solver are avail‐
220 able. Those direct solvers only work with quadratic matrices, so be
221 careful using them with large maps (maps of size 10.000 cells will need
222 more than one gigabyte of ram). Always prefer a sparse matrix solver.
223
225 Use this small python script to create a working groundwater flow /
226 solute transport area and data. Make sure you are not in a lat/lon
227 projection.
228 #!/usr/bin/env python3
229 # This is an example script how groundwater flow and solute transport are
230 # computed within GRASS GIS
231 import sys
232 import os
233 import grass.script as gs
234 # Overwrite existing maps
235 gs.run_command("g.gisenv", set="OVERWRITE=1")
236 gs.message(_("Set the region"))
237 # The area is 200m x 100m with a cell size of 1m x 1m
238 gs.run_command("g.region", res=1, res3=1, t=10, b=0, n=100, s=0, w=0, e=200)
239 gs.run_command("r.mapcalc", expression="phead = if(col() == 1 , 50, 40)")
240 gs.run_command("r.mapcalc", expression="phead = if(col() ==200 , 45 + row()/40, phead)")
241 gs.run_command("r.mapcalc", expression="status = if(col() == 1 || col() == 200 , 2, 1)")
242 gs.run_command("r.mapcalc", expression="well = if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.001, 0)")
243 gs.run_command("r.mapcalc", expression="hydcond = 0.00005")
244 gs.run_command("r.mapcalc", expression="recharge = 0")
245 gs.run_command("r.mapcalc", expression="top_conf = 20")
246 gs.run_command("r.mapcalc", expression="bottom = 0")
247 gs.run_command("r.mapcalc", expression="poros = 0.17")
248 gs.run_command("r.mapcalc", expression="syield = 0.0001")
249 gs.run_command("r.mapcalc", expression="null = 0.0")
250 gs.message(_("Compute a steady state groundwater flow"))
251 gs.run_command("r.gwflow", solver="cg", top="top_conf", bottom="bottom", phead="phead",\
252 status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
253 recharge="recharge", output="gwresult_conf", dt=8640000000000, type="confined")
254 gs.message(_("generate the transport data"))
255 gs.run_command("r.mapcalc", expression="c = if(col() == 15 && row() == 75 , 500.0, 0.0)")
256 gs.run_command("r.mapcalc", expression="cs = if(col() == 15 && row() == 75 , 0.0, 0.0)")
257 gs.run_command("r.mapcalc", expression="tstatus = if(col() == 1 || col() == 200 , 3, 1)")
258 gs.run_command("r.mapcalc", expression="diff = 0.0000001")
259 gs.run_command("r.mapcalc", expression="R = 1.0")
260 # Compute the initial state
261 gs.run_command("r.solute.transport", solver="bicgstab", top="top_conf",\
262 bottom="bottom", phead="gwresult_conf", status="tstatus", hc_x="hydcond", hc_y="hydcond",\
263 rd="R", cs="cs", q="well", nf="poros", output="stresult_conf_0", dt=3600, diff_x="diff",\
264 diff_y="diff", c="c", al=0.1, at=0.01)
265 # Compute the solute transport for 300 days in 10 day steps
266 for dt in range(30):
267 gs.run_command("r.solute.transport", solver="bicgstab", top="top_conf",\
268 bottom="bottom", phead="gwresult_conf", status="tstatus", hc_x="hydcond", hc_y="hydcond",\
269 rd="R", cs="cs", q="well", nf="poros", output="stresult_conf_" + str(dt + 1), dt=864000, diff_x="diff",\
270 diff_y="diff", c="stresult_conf_" + str(dt), al=0.1, at=0.01)
271
273 r.gwflow
274 r3.gwflow
275 r.out.vtk
276
278 Sören Gebbert
279 This work is based on the Diploma Thesis of Sören Gebbert available
281 here at Technical University Berlin in Germany.
282
284 Available at: r.solute.transport source code (history)
285 Accessed: Saturday Jan 21 20:39:01 2023
287 Main index | Raster index | Topics index | Keywords index | Graphical
289 index | Full index
290 © 2003-2023 GRASS Development Team, GRASS GIS 8.2.1 Reference Manual
292GRASS 8.2.1 r.solute.transport(1)