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

NAME

6       r.solute.transport  - Numerical calculation program for transient, con‐
7       fined and unconfined solute transport in two dimensions
8

KEYWORDS

10       raster, hydrology, solute transport
11

SYNOPSIS

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
22   Flags:
23       -f
24           Use  a  full  filled quadratic linear equation system, default is a
25           sparse linear equation system.
26
27       -c
28           Use the Courant-Friedrichs-Lewy criteria for time step calculation
29
30       --overwrite
31           Allow output files to overwrite existing files
32
33       --help
34           Print usage summary
35
36       --verbose
37           Verbose module output
38
39       --quiet
40           Quiet module output
41
42       --ui
43           Force launching GUI dialog
44
45   Parameters:
46       c=name [required]
47           The initial concentration in [kg/m^3]
48
49       phead=name [required]
50           The piezometric head in [m]
51
52       hc_x=name [required]
53           The x-part of the hydraulic conductivity tensor in [m/s]
54
55       hc_y=name [required]
56           The y-part of the hydraulic conductivity tensor in [m/s]
57
58       status=name [required]
59           The status for each cell, = 0 - inactive cell, 1 - active cell, 2 -
60           dirichlet- and 3 - transfer boundary condition
61
62       diff_x=name [required]
63           The x-part of the diffusion tensor in [m^2/s]
64
65       diff_y=name [required]
66           The y-part of the diffusion tensor in [m^2/s]
67
68       q=name
69           Groundwater sources and sinks in [m^3/s]
70
71       cin=name
72           Concentration  sources  and sinks bounded to a water source or sink
73           in [kg/s]
74
75       cs=name [required]
76           Concentration of inner sources and inner sinks in  [kg/s]  (i.e.  a
77           chemical reaction)
78
79       rd=name [required]
80           Retardation factor [-]
81
82       nf=name [required]
83           Effective porosity [-]
84
85       top=name [required]
86           Top surface of the aquifer in [m]
87
88       bottom=name [required]
89           Bottom surface of the aquifer in [m]
90
91       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
95       vx=name
96           Calculate and store the groundwater filter velocity vector part  in
97           x direction [m/s]
98
99       vy=name
100           Calculate  and store the groundwater filter velocity vector part in
101           y direction [m/s]
102
103       dtime=float [required]
104           The calculation time in seconds
105           Default: 86400
106
107       maxit=integer
108           Maximum number of iteration used to solve the linear equation  sys‐
109           tem
110           Default: 10000
111
112       error=float
113           Error break criteria for iterative solver
114           Default: 0.000001
115
116       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
121       relax=float
122           The  relaxation  parameter  used  by  the jacobi and sor solver for
123           speedup or stabilizing
124           Default: 1
125
126       al=float
127           The longditudinal dispersivity length. [m]
128           Default: 0.0
129
130       at=float
131           The transversal dispersivity length. [m]
132           Default: 0.0
133
134       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
139       stab=string
140           Set the flow stabilizing scheme (full or exponential upwinding).
141           Options: full, exp
142           Default: full
143

DESCRIPTION

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
151       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

NOTES

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
177       (dc/dt)*R = div ( D grad c - uc) + cs -q/nf(c - c_in)
178
179           •   c -- the concentration [kg/m^3]
180
181           •   u -- vector of mean groundwater flow velocity
182
183           •   dt -- the time step for transient calculation in seconds [s]
184
185           •   R -- the linear retardation coefficient [-]
186
187           •   D -- the diffusion and dispersion tensor [m^2/s]
188
189           •   cs -- inner concentration sources/sinks [kg/m^3]
190
191           •   c_in -- the solute concentration of influent water [kg/m^3]
192
193           •   q -- inner well sources/sinks [m^3/s]
194
195           •   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
201           •   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
204           •   1 == active - this cell is used for sloute  transport  calcula‐
205               tion, inner sources can be defined for those cells
206
207           •   2  ==  Dirichlet - cells of this type will have a fixed concen‐
208               tration value which do not change over time
209
210           •   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

EXAMPLE

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

SEE ALSO

273       r.gwflow
274       r3.gwflow
275       r.out.vtk
276

AUTHOR

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

SOURCE CODE

284       Available at: r.solute.transport source code (history)
285
286       Accessed: Saturday Jan 21 21:15:10 2023
287
288       Main  index  | Raster index | Topics index | Keywords index | Graphical
289       index | Full index
290
291       © 2003-2023 GRASS Development Team, GRASS GIS 8.2.1 Reference Manual
292
293
294
295GRASS 8.2.1                                              r.solute.transport(1)
Impressum