r.solute.transport(1grass) GRASS GIS User's Manual r.solute.transport(1grass)
NAME
r.solute.transport - Numerical calculation program for transient, con-
fined and unconfined solute transport in two dimensions
KEYWORDS
raster, hydrology, solute transport
SYNOPSIS
r.solute.transport
r.solute.transport --help
r.solute.transport [-fc] c=name phead=name hc_x=name hc_y=name sta-
tus=name diff_x=name diff_y=name [q=name] [cin=name] cs=name
rd=name nf=name top=name bottom=name output=name [vx=name] [vy=name]
dtime=float [maxit=integer] [error=float] [solver=name] [re-
lax=float] [al=float] [at=float] [loops=float] [stab=string]
[--overwrite] [--help] [--verbose] [--quiet] [--ui]
Flags:
-f
Use a full filled quadratic linear equation system, default is a
sparse linear equation system.
-c
Use the Courant-Friedrichs-Lewy criteria for time step calculation
--overwrite
Allow output files to overwrite existing files
--help
Print usage summary
--verbose
Verbose module output
--quiet
Quiet module output
--ui
Force launching GUI dialog
Parameters:
c=name [required]
The initial concentration in [kg/m^3]
phead=name [required]
The piezometric head in [m]
hc_x=name [required]
The x-part of the hydraulic conductivity tensor in [m/s]
hc_y=name [required]
The y-part of the hydraulic conductivity tensor in [m/s]
status=name [required]
The status for each cell, = 0 - inactive cell, 1 - active cell, 2 -
dirichlet- and 3 - transfer boundary condition
diff_x=name [required]
The x-part of the diffusion tensor in [m^2/s]
diff_y=name [required]
The y-part of the diffusion tensor in [m^2/s]
q=name
Groundwater sources and sinks in [m^3/s]
cin=name
Concentration sources and sinks bounded to a water source or sink
in [kg/s]
cs=name [required]
Concentration of inner sources and inner sinks in [kg/s] (i.e. a
chemical reaction)
rd=name [required]
Retardation factor [-]
nf=name [required]
Effective porosity [-]
top=name [required]
Top surface of the aquifer in [m]
bottom=name [required]
Bottom surface of the aquifer in [m]
output=name [required]
The resulting concentration of the numerical solute transport cal-
culation will be written to this map. [kg/m^3]
vx=name
Calculate and store the groundwater filter velocity vector part in
x direction [m/s]
vy=name
Calculate and store the groundwater filter velocity vector part in
y direction [m/s]
dtime=float [required]
The calculation time in seconds
Default: 86400
maxit=integer
Maximum number of iteration used to solve the linear equation sys-
tem
Default: 10000
error=float
Error break criteria for iterative solver
Default: 0.000001
solver=name
The type of solver which should solve the linear equation system
Options: gauss, lu, jacobi, sor, bicgstab
Default: bicgstab
relax=float
The relaxation parameter used by the jacobi and sor solver for
speedup or stabilizing
Default: 1
al=float
The longditudinal dispersivity length. [m]
Default: 0.0
at=float
The transversal dispersivity length. [m]
Default: 0.0
loops=float
Use this number of time loops if the CFL flag is off. The timestep
will become dt/loops.
Default: 1
stab=string
Set the flow stabilizing scheme (full or exponential upwinding).
Options: full, exp
Default: full
DESCRIPTION
This numerical program calculates numerical implicit transient and
steady state solute transport in porous media in the saturated zone of
an aquifer. The computation is based on raster maps and the current re-
gion settings. All initial- and boundary-conditions must be provided as
raster maps. The unit in the location must be meters.
This module is sensitive to mask settings. All cells which are outside
the mask are ignored and handled as no flow boundaries.
This module calculates the concentration of the solution and optional
the velocity field, based on the hydraulic conductivity, the effective
porosity and the initial piecometric heads. The vector components can
be visualized with paraview if they are exported with r.out.vtk.
Use r.gwflow to compute the piezometric heights of the aquifer. The
piezometric heights and the hydraulic conductivity are used to compute
the flow direction and the mean velocity of the groundwater. This is
the base of the solute transport computation.
The solute transport will always be calculated transient. For stady
state computation set the timestep to a large number (billions of sec-
onds).
To reduce the numerical dispersion, which is a consequence of the con-
vection term and the finite volume discretization, you can use small
time steps and choose between full and exponential upwinding.
NOTES
The solute transport calculation is based on a diffusion/convection
partial differential equation and a numerical implicit finite volume
discretization. Specific for this kind of differential equation is the
combination of a diffusion/dispersion term and a convection term. The
discretization results in an unsymmetric linear equation system in form
of Ax = b, which must be solved. The solute transport partial differen-
tial equation is of the following form:
(dc/dt)*R = div ( D grad c - uc) + cs -q/nf(c - c_in)
• c -- the concentration [kg/m^3]
• u -- vector of mean groundwater flow velocity
• dt -- the time step for transient calculation in seconds [s]
• R -- the linear retardation coefficient [-]
• D -- the diffusion and dispersion tensor [m^2/s]
• cs -- inner concentration sources/sinks [kg/m^3]
• c_in -- the solute concentration of influent water [kg/m^3]
• q -- inner well sources/sinks [m^3/s]
• nf -- the effective porosity [-]
Three different boundary conditions are implemented, the Dirichlet,
Transmission and Neumann conditions. The calculation and boundary sta-
tus of single cells can be set with the status map. The following
states are supportet:
• 0 == inactive - the cell with status 0 will not be calculated,
active cells will have a no flow boundary to an inactive cell
• 1 == active - this cell is used for sloute transport calcula-
tion, inner sources can be defined for those cells
• 2 == Dirichlet - cells of this type will have a fixed concen-
tration value which do not change over time
• 3 == Transmission - cells of this type should be placed on
out-flow boundaries to assure the flow of the solute stream out
Note that all required raster maps are read into main memory. Addition-
ally the linear equation system will be allocated, so the memory con-
sumption of this module rapidely grow with the size of the input maps.
The resulting linear equation system Ax = b can be solved with several
solvers. Several iterative solvers with unsymmetric sparse and qua-
dratic matrices support are implemented. The jacobi method, the
Gauss-Seidel method and the biconjugate gradients-stabilized (bicgstab)
method. Additionally a direct Gauss solver and LU solver are avail-
able. Those direct solvers only work with quadratic matrices, so be
careful using them with large maps (maps of size 10.000 cells will need
more than one gigabyte of ram). Always prefer a sparse matrix solver.
EXAMPLE
Use this small python script to create a working groundwater flow /
solute transport area and data. Make sure you are not in a lat/lon
projection.
#!/usr/bin/env python3
# This is an example script how groundwater flow and solute transport are
# computed within grass
import sys
import os
import grass.script as grass
# Overwrite existing maps
grass.run_command("g.gisenv", set="OVERWRITE=1")
grass.message(_("Set the region"))
# The area is 200m x 100m with a cell size of 1m x 1m
grass.run_command("g.region", res=1, res3=1, t=10, b=0, n=100, s=0, w=0, e=200)
grass.run_command("r.mapcalc", expression="phead = if(col() == 1 , 50, 40)")
grass.run_command("r.mapcalc", expression="phead = if(col() ==200 , 45 + row()/40, phead)")
grass.run_command("r.mapcalc", expression="status = if(col() == 1 || col() == 200 , 2, 1)")
grass.run_command("r.mapcalc", expression="well = if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.001, 0)")
grass.run_command("r.mapcalc", expression="hydcond = 0.00005")
grass.run_command("r.mapcalc", expression="recharge = 0")
grass.run_command("r.mapcalc", expression="top_conf = 20")
grass.run_command("r.mapcalc", expression="bottom = 0")
grass.run_command("r.mapcalc", expression="poros = 0.17")
grass.run_command("r.mapcalc", expression="syield = 0.0001")
grass.run_command("r.mapcalc", expression="null = 0.0")
grass.message(_("Compute a steady state groundwater flow"))
grass.run_command("r.gwflow", solver="cg", top="top_conf", bottom="bottom", phead="phead",\
status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
recharge="recharge", output="gwresult_conf", dt=8640000000000, type="confined")
grass.message(_("generate the transport data"))
grass.run_command("r.mapcalc", expression="c = if(col() == 15 && row() == 75 , 500.0, 0.0)")
grass.run_command("r.mapcalc", expression="cs = if(col() == 15 && row() == 75 , 0.0, 0.0)")
grass.run_command("r.mapcalc", expression="tstatus = if(col() == 1 || col() == 200 , 3, 1)")
grass.run_command("r.mapcalc", expression="diff = 0.0000001")
grass.run_command("r.mapcalc", expression="R = 1.0")
# Compute the initial state
grass.run_command("r.solute.transport", solver="bicgstab", top="top_conf",\
bottom="bottom", phead="gwresult_conf", status="tstatus", hc_x="hydcond", hc_y="hydcond",\
rd="R", cs="cs", q="well", nf="poros", output="stresult_conf_0", dt=3600, diff_x="diff",\
diff_y="diff", c="c", al=0.1, at=0.01)
# Compute the solute transport for 300 days in 10 day steps
for dt in range(30):
grass.run_command("r.solute.transport", solver="bicgstab", top="top_conf",\
bottom="bottom", phead="gwresult_conf", status="tstatus", hc_x="hydcond", hc_y="hydcond",\
rd="R", cs="cs", q="well", nf="poros", output="stresult_conf_" + str(dt + 1), dt=864000, diff_x="diff",\
diff_y="diff", c="stresult_conf_" + str(dt), al=0.1, at=0.01)
SEE ALSO
r.gwflow
r3.gwflow
r.out.vtk
AUTHOR
Sören Gebbert
This work is based on the Diploma Thesis of Sören Gebbert available
here at Technical University Berlin in Germany.
SOURCE CODE
Available at: r.solute.transport source code (history)
Accessed: unknown
Main index | Raster index | Topics index | Keywords index | Graphical
index | Full index
© 2003-2022 GRASS Development Team, GRASS GIS 7.8.7 Reference Manual
GRASS 7.8.7 r.solute.transport(1grass)
Generated by dwww version 1.14 on Sun Dec 29 18:57:15 CET 2024.