[GRASS-SVN] r48665 - in grass/trunk: lib/gpde raster/r.gwflow
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Oct 6 17:13:16 EDT 2011
Author: huhabla
Date: 2011-10-06 14:13:16 -0700 (Thu, 06 Oct 2011)
New Revision: 48665
Modified:
grass/trunk/lib/gpde/N_gwflow.c
grass/trunk/raster/r.gwflow/main.c
grass/trunk/raster/r.gwflow/r.gwflow.html
grass/trunk/raster/r.gwflow/valid_calc_7x7.py
grass/trunk/raster/r.gwflow/valid_calc_excavation.py
Log:
Fixed wrong and misleading documentation about the groundwater flow concept.
Modified: grass/trunk/lib/gpde/N_gwflow.c
===================================================================
--- grass/trunk/lib/gpde/N_gwflow.c 2011-10-06 12:56:02 UTC (rev 48664)
+++ grass/trunk/lib/gpde/N_gwflow.c 2011-10-06 21:13:16 UTC (rev 48665)
@@ -307,7 +307,7 @@
/*inner sources */
q = N_get_array_3d_d_value(data->q, col, row, depth);
- /*specific yield */
+ /*storativity */
Ss = N_get_array_3d_d_value(data->s, col, row, depth);
/*porosity */
nf = N_get_array_3d_d_value(data->nf, col, row, depth);
@@ -325,7 +325,7 @@
/*mass balance center cell to bottom cell */
B = -1 * Az * hc_b / dz;
- /*specific yield */
+ /*storativity */
Ss = Az * dz * Ss;
/*the diagonal entry of the matrix */
@@ -474,7 +474,7 @@
double hc_xe, hc_ys;
double z_xe, z_ys;
double hc, hc_start;
- double Ss, r, nf, q;
+ double Ss, r, q;
double C, W, E, N, S, V;
N_gwflow_data2d *data;
N_data_star *mat_pos;
@@ -495,7 +495,15 @@
hc = N_get_array_2d_d_value(data->phead, col, row);
top = N_get_array_2d_d_value(data->top, col, row);
+ /* Inner sources */
+ q = N_get_array_2d_d_value(data->q, col, row);
+ /* storativity or porosity of current cell face [-]*/
+ Ss = N_get_array_2d_d_value(data->s, col, row);
+ /* recharge */
+ r = N_get_array_2d_d_value(data->r, col, row) * Az;
+
+
if (hc > top) { /*If the aquifer is confined */
z = N_get_array_2d_d_value(data->top, col,
row) -
@@ -551,15 +559,6 @@
else
z_s = z;
- /* Inner sources */
- q = N_get_array_2d_d_value(data->q, col, row);
- nf = N_get_array_2d_d_value(data->nf, col, row);
-
- /* specific yield of current cell face */
- Ss = N_get_array_2d_d_value(data->s, col, row) * Az;
- /* recharge */
- r = N_get_array_2d_d_value(data->r, col, row) * Az;
-
/*get the surrounding permeabilities */
hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
@@ -575,16 +574,17 @@
T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
/* Compute the river leakage, this is an explicit method
- * Rivers are only enabled, if the river bed is lower or equal to the surface
+ * Influent and effluent flow is computed.
*/
if (data->river_leak &&
(N_get_array_2d_d_value(data->river_leak, col, row) != 0) &&
N_get_array_2d_d_value(data->river_bed, col, row) <= top) {
+ /* Groundwater surface is above the river bed*/
if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
N_get_array_2d_d_value(data->river_leak, col, row);
river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
- }
+ } /* Groundwater surface is below the river bed */
else if (hc < N_get_array_2d_d_value(data->river_bed, col, row)) {
river_vect = (N_get_array_2d_d_value(data->river_head, col, row) -
N_get_array_2d_d_value(data->river_bed, col, row))
@@ -594,7 +594,7 @@
}
/* compute the drainage, this is an explicit method
- * Drainage is only enabled, if the drain bed is lower or equal to the surface
+ * Drainage is only enabled, if the drain bed is lower the groundwater surface
*/
if (data->drain_leak &&
(N_get_array_2d_d_value(data->drain_leak, col, row) != 0) &&
@@ -620,11 +620,11 @@
S = -1 * T_s * dx / dy;
/*the diagonal entry of the matrix */
- C = -1 * (W + E + N + S - Ss / data->dt - river_mat * Az -
+ C = -1 * (W + E + N + S - Az *Ss / data->dt - river_mat * Az -
drain_mat * Az);
/*the entry in the right side b of Ax = b */
- V = (q + hc_start * Ss / data->dt) + r + river_vect * Az +
+ V = (q + hc_start * Az * Ss / data->dt) + r + river_vect * Az +
drain_vect * Az;
G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
Modified: grass/trunk/raster/r.gwflow/main.c
===================================================================
--- grass/trunk/raster/r.gwflow/main.c 2011-10-06 12:56:02 UTC (rev 48664)
+++ grass/trunk/raster/r.gwflow/main.c 2011-10-06 21:13:16 UTC (rev 48665)
@@ -78,7 +78,7 @@
param.s = G_define_standard_option(G_OPT_R_INPUT);
param.s->key = "s";
- param.s->description = _("Input raster map with specific yield in [1/m]");
+ param.s->description = _("Input raster map with storativity for confined or effective porosity for unconfined groundwater flow booth in [-] ");
param.r = G_define_standard_option(G_OPT_R_INPUT);
param.r->key = "recharge";
@@ -103,19 +103,19 @@
param.vector_x->key = "vx";
param.vector_x->required = NO;
param.vector_x->description =
- _("Output raster map storing the groundwater filter velocity vector part in x direction [m/s]\n");
+ _("Output raster map to store the groundwater filter velocity vector part in x direction [m/s]\n");
param.vector_y = G_define_standard_option(G_OPT_R_OUTPUT);
param.vector_y->key = "vy";
param.vector_y->required = NO;
param.vector_y->description =
- _("Output raster map storing the groundwater filter velocity vector part in y direction [m/s]\n");
+ _("Output raster map to store the groundwater filter velocity vector part in y direction [m/s]\n");
param.budget = G_define_standard_option(G_OPT_R_OUTPUT);
param.budget->key = "budget";
param.budget->required = NO;
param.budget->description =
- _("Output raster map storing the groundwater budget for each cell [m^3/s]\n");
+ _("Output raster map to store the groundwater budget for each cell [m^3/s]\n");
param.type = G_define_option();
param.type->key = "type";
Modified: grass/trunk/raster/r.gwflow/r.gwflow.html
===================================================================
--- grass/trunk/raster/r.gwflow/r.gwflow.html 2011-10-06 12:56:02 UTC (rev 48664)
+++ grass/trunk/raster/r.gwflow/r.gwflow.html 2011-10-06 21:13:16 UTC (rev 48665)
@@ -28,7 +28,7 @@
The groundwater flow will always be calculated transient.
For stady state computation set the timestep
to a large number (billions of seconds) or set the
-specific yield/ effective porosity raster map to zero.
+storativity/ effective porosity raster map to zero.
<br>
<br>
The water budget is calculated for each non inactive cell. The
@@ -52,14 +52,19 @@
<ul>
<li>h -- the piezometric head im [m]</li>
<li>dt -- the time step for transient calculation in [s]</li>
-<li>S -- the specific yield [1/m]</li>
+<li>S -- the specific storage [1/m]</li>
<li>Kxx -- the hydraulic conductivity tensor part in x direction in [m/s]</li>
<li>Kyy -- the hydraulic conductivity tensor part in y direction in [m/s]</li>
<li>q - inner source/sink in meter per second [1/s]</li>
</ul>
-<br>
-<br>
+<p>
+Confined and unconfined groundwater flow is supported. Be aware that the storativity input parameter
+is handled differently in case of unconfined flow. Instead of the storativity, the effective porosity is expected.
+<p>
+To compute unconfined groundwater flow, a simple Picard based linearization scheme is used to
+solve the resulting non-linear equation system.
+<p>
Two different boundary conditions are implemented,
the Dirichlet and Neumann conditions. By default the calculation area is surrounded by homogeneous Neumann boundary conditions.
The calculation and boundary status of single cells must be set with a status map,
@@ -105,7 +110,7 @@
r.mapcalc --o expression="bottom=0.0"
r.mapcalc --o expression="null=0.0"
r.mapcalc --o expression="poros=0.15"
-r.mapcalc --o expression="syield=0.0001"
+r.mapcalc --o expression="s=0.0001"
# The maps of the river
r.mapcalc --o expression="river_bed=if(col() == 35 , 48, null())"
@@ -119,10 +124,11 @@
#confined groundwater flow with cg solver and sparse matrix, river and drain
#do not work with this confined aquifer (top == 20m)
r.gwflow --o 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 \
+hc_x=hydcond hc_y=hydcond q=well s=s recharge=recharge output=gwresult_conf \
dt=8640000 type=confined vx=gwresult_conf_velocity_x vy=gwresult_conf_velocity_y budget=budget_conf
#unconfined groundwater flow with cg solver and sparse matrix, river and drain are enabled
+# We use the effective porosity as storativity parameter
r.gwflow --o solver=cg top=top_unconf bottom=bottom phead=phead \
status=status hc_x=hydcond hc_y=hydcond q=well s=poros recharge=recharge \
river_bed=river_bed river_head=river_head river_leak=river_leak \
Modified: grass/trunk/raster/r.gwflow/valid_calc_7x7.py
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_7x7.py 2011-10-06 12:56:02 UTC (rev 48664)
+++ grass/trunk/raster/r.gwflow/valid_calc_7x7.py 2011-10-06 21:13:16 UTC (rev 48665)
@@ -25,19 +25,19 @@
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="syield=0.0001")
+grass.run_command("r.mapcalc", expression="s=0.0001")
grass.run_command("r.mapcalc", expression="null=0.0")
-#First compute the initial groundwater flow
+#First compute the groundwater flow
grass.run_command("r.gwflow", "f", solver="cholesky", top="top_conf", bottom="bottom", phead="phead",\
- status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
+ status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="s",\
recharge="recharge", output="gwresult_conf", dt=500, type="confined", budget="water_budget")
count=500
# loop over the timesteps
for i in range(20):
grass.run_command("r.gwflow", "f", solver="cholesky", top="top_conf", bottom="bottom", phead="gwresult_conf",\
- status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
+ status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="s",\
recharge="recharge", output="gwresult_conf", dt=500, type="confined", budget="water_budget")
count += 500
Modified: grass/trunk/raster/r.gwflow/valid_calc_excavation.py
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_excavation.py 2011-10-06 12:56:02 UTC (rev 48664)
+++ grass/trunk/raster/r.gwflow/valid_calc_excavation.py 2011-10-06 21:13:16 UTC (rev 48665)
@@ -29,10 +29,10 @@
grass.run_command("r.mapcalc", expression="recharge=0.000000006")
grass.run_command("r.mapcalc", expression="top=20")
grass.run_command("r.mapcalc", expression="bottom=0")
-grass.run_command("r.mapcalc", expression="syield=0.001")
+grass.run_command("r.mapcalc", expression="poros=0.1")
grass.run_command("r.mapcalc", expression="null=0.0")
#compute a steady state groundwater flow
grass.run_command("r.gwflow", "f", solver="cholesky", top="top", bottom="bottom", phead="phead", \
- status="status", hc_x="hydcond", hc_y="hydcond", s="syield", \
+ status="status", hc_x="hydcond", hc_y="hydcond", s="poros", \
recharge="recharge", output="gwresult", dt=864000000000, type="unconfined", budget="water_budget")
More information about the grass-commit
mailing list