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

**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 calculation 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 system
- 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*

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

(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

#!/usr/bin/env python3 # This is an example script how groundwater flow and solute transport are # computed within GRASS GIS import sys import os import grass.script as gs # Overwrite existing maps gs.run_command("g.gisenv", set="OVERWRITE=1") gs.message(_("Set the region")) # The area is 200m x 100m with a cell size of 1m x 1m gs.run_command("g.region", res=1, res3=1, t=10, b=0, n=100, s=0, w=0, e=200) gs.run_command("r.mapcalc", expression="phead = if(col() == 1 , 50, 40)") gs.run_command("r.mapcalc", expression="phead = if(col() ==200 , 45 + row()/40, phead)") gs.run_command("r.mapcalc", expression="status = if(col() == 1 || col() == 200 , 2, 1)") gs.run_command("r.mapcalc", expression="well = if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.001, 0)") gs.run_command("r.mapcalc", expression="hydcond = 0.00005") gs.run_command("r.mapcalc", expression="recharge = 0") gs.run_command("r.mapcalc", expression="top_conf = 20") gs.run_command("r.mapcalc", expression="bottom = 0") gs.run_command("r.mapcalc", expression="poros = 0.17") gs.run_command("r.mapcalc", expression="syield = 0.0001") gs.run_command("r.mapcalc", expression="null = 0.0") gs.message(_("Compute a steady state groundwater flow")) gs.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") gs.message(_("generate the transport data")) gs.run_command("r.mapcalc", expression="c = if(col() == 15 && row() == 75 , 500.0, 0.0)") gs.run_command("r.mapcalc", expression="cs = if(col() == 15 && row() == 75 , 0.0, 0.0)") gs.run_command("r.mapcalc", expression="tstatus = if(col() == 1 || col() == 200 , 3, 1)") gs.run_command("r.mapcalc", expression="diff = 0.0000001") gs.run_command("r.mapcalc", expression="R = 1.0") # Compute the initial state gs.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): gs.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)

This work is based on the Diploma Thesis of Sö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

© 2003-2020 GRASS Development Team, GRASS GIS 7.9.dev Reference Manual