Fixed wrong and misleading documentation about the groundwater flow concept.

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@48665 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/gpde/N_gwflow.c

@@ -307,7 +307,7 @@ N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data * geom, int col,
 
     /*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 @@ N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data * geom, int col,
     /*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 @@ N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
     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,6 +495,14 @@ N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
     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,
@@ -551,15 +559,6 @@ N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
     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 @@ N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
     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 @@ N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
     }
 
     /* 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 @@ N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
     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);

raster/r.gwflow/main.c

@@ -78,7 +78,7 @@ void set_params(void)
 
     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 @@ void set_params(void)
     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";

raster/r.gwflow/r.gwflow.html

@@ -28,7 +28,7 @@ with <em>r.out.vtk</em>.
 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 @@ In detail for 2 dimensions:
 <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="top_unconf=70.0"
 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 @@ r.mapcalc --o expression="drain_leak=if(col() == 5 , 0.01, null())"
 #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 \

raster/r.gwflow/valid_calc_7x7.py

@@ -25,19 +25,19 @@ grass.run_command("r.mapcalc", expression="hydcond=0.0005")
 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
 

raster/r.gwflow/valid_calc_excavation.py

@@ -29,10 +29,10 @@ grass.run_command("r.mapcalc", expression="hydcond=0.001")
 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")