Skip site navigation (1)Skip section navigation (2)

FreeBSD Manual Pages


home | help
r.solute.transport(1)	    GRASS GIS User's Manual	 r.solute.transport(1)

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

       raster, hydrology, 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]

	   Use	a  full	 filled	quadratic linear equation system, default is a
	   sparse linear equation system.

	   Use the Courant-Friedrichs-Lewy criteria for	time step calculation

	   Allow output	files to overwrite existing files

	   Print usage summary

	   Verbose module output

	   Quiet module	output

	   Force launching GUI dialog

       c=nameA [required]
	   The initial concentration in	[kg/m^3]

       phead=nameA [required]
	   The piezometric head	in [m]

       hc_x=nameA [required]
	   The x-part of the hydraulic conductivity tensor in [m/s]

       hc_y=nameA [required]
	   The y-part of the hydraulic conductivity tensor in [m/s]

       status=nameA [required]
	   The status for each cell, = 0 - inactive cell, 1 - active cell, 2 -
	   dirichlet- and 3 - transfer boundary	condition

       diff_x=nameA [required]
	   The x-part of the diffusion tensor in [m^2/s]

       diff_y=nameA [required]
	   The y-part of the diffusion tensor in [m^2/s]

	   Groundwater sources and sinks in [m^3/s]

	   Concentration  sources  and sinks bounded to	a water	source or sink
	   in [kg/s]

       cs=nameA	[required]
	   Concentration of inner sources and inner sinks in  [kg/s]  (i.e.  a
	   chemical reaction)

       rd=nameA	[required]
	   Retardation factor [-]

       nf=nameA	[required]
	   Effective porosity [-]

       top=nameA [required]
	   Top surface of the aquifer in [m]

       bottom=nameA [required]
	   Bottom surface of the aquifer in [m]

       output=nameA [required]
	   The	resulting concentration	of the numerical solute	transport cal-
	   culation will be written to this map. [kg/m^3]

	   Calculate and store the groundwater filter velocity vector part  in
	   x direction [m/s]

	   Calculate  and store	the groundwater	filter velocity	vector part in
	   y direction [m/s]

       dtime=floatA [required]
	   The calculation time	in seconds
	   Default: 86400

	   Maximum number of iteration used to solve the linear	equation  sys-
	   Default: 10000

	   Error break criteria	for iterative solver
	   Default: 0.000001

	   The type of solver which should solve the linear equation system
	   Options: gauss, lu, jacobi, sor, bicgstab
	   Default: bicgstab

	   The	relaxation  parameter  used  by	 the jacobi and	sor solver for
	   speedup or stabilizing
	   Default: 1

	   The longditudinal dispersivity length. [m]
	   Default: 0.0

	   The transversal dispersivity	length.	[m]
	   Default: 0.0

	   Use this number of time loops if the	CFL flag is off. The  timestep
	   will	become dt/loops.
	   Default: 1

	   Set the flow	stabilizing scheme (full or exponential	upwinding).
	   Options: full, exp
	   Default: full

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

       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)

	   o   c -- the	concentration [kg/m^3]

	   o   u -- vector of mean groundwater flow velocity

	   o   dt -- the time step for transient calculation in	seconds	[s]

	   o   R -- the	linear retardation coefficient [-]

	   o   D -- the	diffusion and dispersion tensor	[m^2/s]

	   o   cs -- inner concentration sources/sinks [kg/m^3]

	   o   c_in -- the solute concentration	of influent water [kg/m^3]

	   o   q -- inner well sources/sinks [m^3/s]

	   o   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:

	   o   0  == inactive -	the cell with status 0 will not	be calculated,
	       active cells will have a	no flow	boundary to an inactive	cell

	   o   1 == active - this cell is used for sloute  transport  calcula-
	       tion, inner sources can be defined for those cells

	   o   2  ==  Dirichlet	- cells	of this	type will have a fixed concen-
	       tration value which do not change over time

	   o   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.

       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
       #!/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)


       SA<paragraph>ren	Gebbert

       This work is based on the Diploma Thesis	 of  SA<paragraph>ren  Gebbert
       available here at Technical University Berlin in	Germany.

       Available at: r.solute.transport	source code (history)

       Main  index  | Raster index | Topics index | Keywords index | Graphical
       index | Full index

       A(C) 2003-2020 GRASS Development	Team, GRASS GIS	7.8.4 Reference	Manual

GRASS 7.8.4						 r.solute.transport(1)


Want to link to this manual page? Use this URL:

home | help