[GRASS-SVN] r70765 - grass/trunk/raster/r.topmodel
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Mar 18 07:48:37 PDT 2017
Author: hcho
Date: 2017-03-18 07:48:37 -0700 (Sat, 18 Mar 2017)
New Revision: 70765
Modified:
grass/trunk/raster/r.topmodel/infiltration.c
grass/trunk/raster/r.topmodel/r.topmodel.html
Log:
r.topmodel: Rename variables for infiltration calculation based on the original paper; Add more comments to the infiltration routine; Add more references
Modified: grass/trunk/raster/r.topmodel/infiltration.c
===================================================================
--- grass/trunk/raster/r.topmodel/infiltration.c 2017-03-17 21:29:45 UTC (rev 70764)
+++ grass/trunk/raster/r.topmodel/infiltration.c 2017-03-18 14:48:37 UTC (rev 70765)
@@ -3,64 +3,125 @@
#include <grass/glocale.h>
#include "global.h"
-/* The Green-and-Ampt Model */
+/* Calculates the infiltrate rate (m/hr) for the given timestep. For variable
+ * names and equation numbers in comments, refer to Beven (1984).
+ *
+ * Beven, K. J., 1984. Infiltration into a class of vertically non-uniform
+ * soils. Hydrological Sciences Journal 29 (4), 425-434.
+ *
+ * Beven, K. J., Kirkby, M. J., 1979. A physically based, variable contributing
+ * area model of basin hydrology. Hydrological Sciences Bulletin 24 (1), 43-69.
+ *
+ * Morel-Seytoux, H. J., Khanji, J., 1974. Derivation of an equation of
+ * infiltration. Water Resources Research 10 (4), 795-800.
+ */
double calculate_infiltration(int timestep, double R)
{
- static double cumf = 0.0, f = 0.0;
+ /* params.K0 Surface hydraulic conductivity (m/h)
+ * params.psi Wetting front suction (m)
+ * params.dtheta Water content change across the wetting front
+ * dtheta = saturated moisture content
+ * - initial moisture content
+ * params.m Parameter controlling the decline rate of
+ * transmissivity (m)
+ *
+ * Beven and Kirkby (1979) introduced the scaling
+ * parameter m.
+ *
+ * K(z) = K0 * exp(-f * z)
+ *
+ * where K(z) is hydraulic conductivity at depth z,
+ * z is the soil depth
+ * f is the parameter controlling the decline rate
+ * of transmissivity (1/m); can be defined by m as
+ * f = dtheta / m
+ *
+ * Now, m = dtheta / f.
+ *
+ * R Rainfall intensity (m/h)
+ * r Infiltration rate (m/h)
+ * cumI Cumulative infiltration at the start of timestep (m)
+ * I Cumulative infiltration at the end of timestep (m)
+ * dIdt Infiltration rate for the current timestep (m/hr)
+ * C Storage-suction factor (m) (Morel-Seytoux and Khanji,
+ * 1974); C = psi * dtheta
+ * IC I + C
+ * lambda lambda in Eq. (8); Note that this lambda is different
+ * from params.lambda
+ * t Current time (hr)
+ * tp Time to ponding (hr)
+ */
+ static double cumI = 0.0, I = 0.0;
static char ponding = 0;
- double t, df, f1, f2, fc, R2, cnst, pt, psi_dtheta, sum;
- int factorial;
+ double r, C, IC, lambda, dIdt, t, tp;
+ double f1, f2, df, sum;
+ int fact;
int i, j;
/* reset if there is no rainfall */
if (R <= 0.0) {
- cumf = 0.0;
- f = 0.0;
+ cumI = 0.0;
+ I = 0.0;
ponding = 0;
return 0.0;
}
t = timestep * input.dt;
- f1 = cnst = pt = 0.0;
- psi_dtheta = params.psi * params.dtheta;
+ lambda = tp = f1 = 0.0;
+ C = params.psi * params.dtheta;
+
+ /* if no ponding occurred yet */
if (!ponding) {
- if (cumf) {
- f1 = cumf;
- R2 = -params.K0 / params.m *
- (psi_dtheta + f1) / (1 - exp(f1 / params.m));
+ if (cumI > 0) {
+ f1 = cumI;
+ /* infiltration rate: Eq. (6)
+ * Note that his Ks = K0 * exp(f * z) in Eq. (1a) instead of
+ * Ks = K0 * exp(-f * z) used in his TOPMODEL code, TMOD9502.F.
+ * Substitute f=-dtheta/m for f in Eq. (1a), hence -K0 and
+ * exp(f1/m), slightly different from the original Eq. (6).
+ */
+ r = -params.K0 / params.m * (C + f1) / (1 - exp(f1 / params.m));
/* rainfall intensity is greater than infiltration
* rate, so ponding starts */
- if (R2 < R) {
- f = cumf;
- pt = t - input.dt;
+ if (r < R) {
+ I = cumI;
+ tp = t - input.dt;
ponding = 1;
goto cont1;
}
}
- f2 = cumf + R * input.dt;
- R2 = -params.K0 / params.m *
- (psi_dtheta + f2) / (1 - exp(f2 / params.m));
- /* rainfall intensity is less than infiltration rate. all
- * rainfall will be infiltrated. */
- if (f2 == 0.0 || R2 > R)
+ /* try full infiltration */
+ f2 = cumI + R * input.dt;
+ /* infiltration rate */
+ r = -params.K0 / params.m * (C + f2) / (1 - exp(f2 / params.m));
+ /* if infiltration rate is greater than rainfall intensity, all
+ * rainfall will be infiltrated */
+ if (f2 == 0.0 || r > R)
goto cont2;
- /* rainfall intensity is greater than infiltration rate. */
- f = cumf + R2 * input.dt;
+ /* infiltration rate is less than rainfall intensity */
+ /* Newton-Raphson iteration to solve Eq. (6) for f */
+ /* guess new cumulative infiltration */
+ I = cumI + r * input.dt;
for (i = 0; i < MAXITER; i++) {
- R2 = -params.K0 / params.m *
- (psi_dtheta + f) / (1 - exp(f / params.m));
- if (R2 > R) {
- f1 = f;
- f = (f + f2) / 2.0;
- df = f - f1;
+ /* new infiltration rate */
+ r = -params.K0 / params.m * (C + I) / (1 - exp(I / params.m));
+ /* if new infiltration rate is greater than rainfall intensity,
+ * increase cumulative infiltration
+ */
+ if (r > R) {
+ f1 = I;
+ I = (I + f2) / 2.0;
+ df = I - f1;
}
+ /* otherwise, decrease cumulative infiltration */
else {
- f2 = f;
- f = (f + f1) / 2.0;
- df = f - f2;
+ f2 = I;
+ I = (I + f1) / 2.0;
+ df = I - f2;
}
+ /* stop if cumulative infiltration converged */
if (fabs(df) <= TOLERANCE)
break;
}
@@ -69,38 +130,58 @@
_("Maximum number of iterations exceeded at timestep %d"),
timestep);
- pt = t - input.dt + (f - cumf) / R;
- if (pt > t)
+ /* ponding time */
+ tp = t - input.dt + (I - cumI) / R;
+ /* if ponding time is greater than the current time,
+ * tp = t - dt + (I - cumI) / R > t
+ * (I - cumI) / R > dt
+ * I - cumI > R * dt
+ * means that additional infiltration (I - cumI) is greater than the
+ * total rainfall (R * dt), which cannot happen when there is no
+ * ponding, so infiltrate all rainfall
+ */
+ if (tp > t)
goto cont2;
cont1:
- cnst = 0.0;
- factorial = 1;
- fc = f + psi_dtheta;
+ /* if additional infiltration is less than the total rainfall, ponding
+ * starts
+ */
+ lambda = 0.0;
+ fact = 1;
+ IC = I + C;
for (j = 1; j <= NTERMS; j++) {
- factorial *= j;
- cnst += pow(fc / params.m, (double)j) / (double)(j * factorial);
+ fact *= j;
+ lambda += pow(IC / params.m, (double)j) / (double)(j * fact);
}
- cnst = log(fc) - (log(fc) + cnst) / exp(psi_dtheta / params.m);
- f += R * (t - pt) / 2.0;
+ /* lambda in Eq. (8) */
+ lambda = log(IC) - (log(IC) + lambda) / exp(C / params.m);
+ I += R * (t - tp) / 2.0;
ponding = 1;
}
- /* Newton-Raphson */
+ /* Newton-Raphson iteration to solve Eq. (8) for f */
for (i = 0; i < MAXITER; i++) {
- fc = f + psi_dtheta;
+ IC = I + C;
sum = 0.0;
- factorial = 1;
+ fact = 1;
for (j = 1; j <= NTERMS; j++) {
- factorial *= j;
- sum += pow(fc / params.m, (double)j) / (double)(j * factorial);
+ fact *= j;
+ sum += pow(IC / params.m, (double)j) / (double)(j * fact);
}
- f1 = -(log(fc) - (log(fc) + sum) /
- exp(psi_dtheta / params.m) - cnst) /
- (params.K0 / params.m) - (t - pt);
- f2 = (exp(f / params.m) - 1.0) / (fc * params.K0 / params.m);
+ /* Eq. (8) - (t - tp) in hr: should converge to 0
+ * Note that sum is outside 1/exp(C/m) in Eq. (8), but inside in his
+ * TMOD9502.F. Based on lambda and his code, it looks like a typo in
+ * Eq. (8).
+ */
+ f1 = -(log(IC) - (log(IC) + sum) / exp(C / params.m) - lambda) /
+ (params.K0 / params.m) - (t - tp);
+ /* inverse of Eq. (7) in hr/m */
+ f2 = (exp(I / params.m) - 1.0) / (IC * params.K0 / params.m);
+ /* -(Eq. (8) - (t-tp)) * Eq. (7): cumulative infiltration in a short
+ * time period */
df = -f1 / f2;
- f += df;
+ I += df;
if (fabs(df) <= TOLERANCE)
break;
}
@@ -108,18 +189,23 @@
G_warning(_("Maximum number of iterations exceeded at timestep %d"),
timestep);
- if (f < cumf + R * input.dt) {
- df = (f - cumf) / input.dt;
- cumf = f;
+ /* if new cumulative infiltration is less than the previous cumulative
+ * infiltration plus the total rainfall, update the current cumulative
+ * infiltration and guess cumulative infiltration for the next time step
+ */
+ if (I < cumI + R * input.dt) {
+ dIdt = (I - cumI) / input.dt;
+ cumI = I;
/* initial guess for next time step */
- f += df * input.dt;
- return df;
+ I += dIdt * input.dt;
+ return dIdt;
}
cont2:
- df = R;
- cumf += df * input.dt;
+ /* infiltrate all rainfall */
+ dIdt = R;
+ cumI += dIdt * input.dt;
ponding = 0;
- return df;
+ return dIdt;
}
Modified: grass/trunk/raster/r.topmodel/r.topmodel.html
===================================================================
--- grass/trunk/raster/r.topmodel/r.topmodel.html 2017-03-17 21:29:45 UTC (rev 70764)
+++ grass/trunk/raster/r.topmodel/r.topmodel.html 2017-03-18 14:48:37 UTC (rev 70765)
@@ -28,7 +28,8 @@
# lnTe [ln(m^2/h)]: Areal average of the soil surface transmissivity
4.
-# m [m]: Scaling parameter
+# m [m]: Parameter controlling the decline rate of transmissivity
+# See Beven and Kirkby (1979)
0.0125
# Sr0 [m]: Initial root zone storage deficit
@@ -132,16 +133,27 @@
<h2>REFERENCES</h2>
<ul>
- <li>Cho, H., 2000. GIS Hydrological Modeling System by Using Programming Interface
- of GRASS. Master's Thesis, Department of Civil Engineering, Kyungpook National
- University, Korea.
- <li>Beven K., R. Lamb, P. Quinn, R. Romanowicz, and J. Freer, 1995. TOPMODEL, in
- V.P. Singh (Ed.). Computer Models of Watershed Hydrology. Water Resources
- Publications.
- <li>
- Liaw, S.C., 1988. Streamflow Simulation Using a Physically Based Hydrologic
- Model in Humid Forested Watersheds. Dissertation, Colorado State University,
- CO. p163.
+ <li>Beven, K. J., 1984. Infiltration into a class of vertically non-uniform
+ soils. Hydrological Sciences Journal 29 (4), 425-434.
+
+ <li>Beven, K. J., Kirkby, M. J., 1979. A physically based, variable
+ contributing area model of basin hydrology. Hydrological Sciences Bulletin 24
+ (1), 43-69.
+
+ <li>Beven K. J., R. Lamb, P. Quinn, R. Romanowicz, and J. Freer, 1995.
+ TOPMODEL, in V.P. Singh (Ed.). Computer Models of Watershed Hydrology. Water
+ Resources Publications.
+
+ <li>Cho, H., 2000. GIS Hydrological Modeling System by Using Programming
+ Interface of GRASS. Master's Thesis, Department of Civil Engineering,
+ Kyungpook National University, Korea.
+
+ <li>Liaw, S. C., 1988. Streamflow Simulation Using a Physically Based
+ Hydrologic Model in Humid Forested Watersheds. Dissertation, Colorado State
+ University, CO. p163.
+
+ <li>Morel-Seytoux, H. J., Khanji, J., 1974. Derivation of an equation of
+ infiltration. Water Resources Research 10 (4), 795-800.
</ul>
@@ -157,7 +169,7 @@
<h2>AUTHORS</h2>
-<a href="mailto:grass4u at gmail com">Huidae Cho</a>,
+<a href="mailto:grass4u at gmail com">Huidae Cho</a>,
Hydro Laboratory, Kyungpook National University, South Korea
<p>
Based on TMOD9502.FOR by Keith Beven.
More information about the grass-commit
mailing list