Provided by: grass-doc_7.6.0-1_all bug


       r.solute.transport  - Numerical calculation program for transient, confined 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  status=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]    [relax=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=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]

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

           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  calculation  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=float [required]
           The calculation time in seconds
           Default: 86400

           Maximum number of iteration used to solve the linear equation system
           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
           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 region 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
       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 seconds).
       To reduce the numerical dispersion, which is a consequence of the convection 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  differential  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 status 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  calculation,  inner  sources
               can be defined for those cells

           ·   2 == Dirichlet - cells of this type will have a fixed concentration 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. Additionally the linear
       equation system will be allocated, so the memory consumption 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 quadratic 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 available.  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 projection.
       #!/usr/bin/env python
       # 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)




       Sören Gebbert

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

       Last changed: $Date: 2017-02-17 22:14:15 +0100 (Fri, 17 Feb 2017) $


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

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

       © 2003-2019 GRASS Development Team, GRASS GIS 7.6.0 Reference Manual