1 #include <math.h>
2 #include <grass/gis.h>
3 #include <grass/glocale.h>
4 #include "global.h"
5
6 #define TOLERANCE 0.00001
7 #define MAX_ITER 20
8 #define NUM_TERMS 10
9 #define PONDING_NO 0
10 #define PONDING_STARTED 1
11 #define PONDING_YES 2
12
13 /* Calculates the infiltrate rate (m/hr) for the given time step. For variable
14 * names and equation numbers in comments, refer to Beven (1984).
15 *
16 * Beven, K. J., 1984. Infiltration into a class of vertically non-uniform
17 * soils. Hydrological Sciences Journal 29 (4), 425-434.
18 *
19 * Beven, K. J., Kirkby, M. J., 1979. A physically based, variable contributing
20 * area model of basin hydrology. Hydrological Sciences Bulletin 24 (1), 43-69.
21 *
22 * Morel-Seytoux, H. J., Khanji, J., 1974. Derivation of an equation of
23 * infiltration. Water Resources Research 10 (4), 795-800.
24 */
calculate_infiltration(int timestep,double R)25 double calculate_infiltration(int timestep, double R)
26 {
27 /* params.K0 Surface hydraulic conductivity (m/h)
28 * params.psi Wetting front suction (m)
29 * params.dtheta Water content change across the wetting front
30 * dtheta = saturated moisture content
31 * - initial moisture content
32 * params.m Parameter controlling the decline rate of
33 * transmissivity (m)
34 *
35 * Beven and Kirkby (1979) introduced the scaling
36 * parameter m.
37 *
38 * K(z) = K0 * exp(-f * z)
39 *
40 * where K(z) is hydraulic conductivity at depth z,
41 * z is the soil depth, and
42 * f is the parameter controlling the decline rate
43 * of transmissivity (1/m); can be defined by m as
44 * f = dtheta / m
45 *
46 * Now, m = dtheta / f.
47 *
48 * R Rainfall intensity (m/h)
49 * r Infiltration rate (m/h)
50 * cumI Cumulative infiltration at the start of time step (m)
51 * I Cumulative infiltration at the end of time step (m)
52 * dIdt Infiltration rate for the current time step (m/hr)
53 * C Storage-suction factor (m) (Morel-Seytoux and Khanji,
54 * 1974); C = psi * dtheta
55 * IC I + C
56 * lambda lambda in Eq. (8); Note that this lambda is different
57 * from params.lambda
58 * t Current time (hr)
59 * tp Time to ponding (hr)
60 * ponding Ponding indicator
61 * PONDING_NO: no ponding
62 * PONDING_STARTED: ponding started in this time step
63 * PONDING_YES: ponding has started in a previous time step
64 */
65 static double cumI = 0.0, I = 0.0, lambda = 0.0, tp = 0.0;
66 static int ponding = PONDING_NO;
67 double r, C, IC, dIdt, t;
68 double f1, f2, df, sum;
69 int fact;
70 int i, j;
71
72 /* reset if there is no rainfall */
73 if (R <= 0.0) {
74 cumI = I = lambda = tp = 0.0;
75 ponding = PONDING_NO;
76 return 0.0;
77 }
78
79 t = timestep * input.dt;
80 C = params.psi * params.dtheta;
81 f1 = 0.0;
82
83 /* if ponding hasn't started and cumulative infiltration is greater than 0
84 */
85 if (ponding == PONDING_NO && cumI > 0) {
86 f1 = cumI;
87 /* infiltration rate: Eq. (6)
88 * Note that his Ks=K0*exp(f*z) in Eq. (1a) instead of Ks=K0*exp(-f*z)
89 * used in his TOPMODEL code, TMOD9502.F. Substitute f=-dtheta/m for f
90 * in Eq. (1a), hence -K0 and exp(f1/m), slightly different from the
91 * original Eq. (6).
92 */
93 r = -params.K0 / params.m * (C + f1) / (1 - exp(f1 / params.m));
94 /* if infiltration rate is less than rainfall intensity, ponding starts
95 */
96 if (r < R) {
97 I = cumI;
98 /* ponding time: tp will remain the same until next ponding occurs
99 */
100 tp = t - input.dt;
101 ponding = PONDING_STARTED;
102 }
103 }
104
105 /* if ponding hasn't started yet */
106 if (ponding == PONDING_NO) {
107 /* try full infiltration */
108 f2 = cumI + R * input.dt;
109 /* infiltration rate */
110 r = -params.K0 / params.m * (C + f2) / (1 - exp(f2 / params.m));
111 /* if potential cumulative infiltration is 0 or infiltration rate is
112 * greater than rainfall intensity, all rainfall will be infiltrated
113 */
114 if (f2 == 0.0 || r > R) {
115 /* full infiltration */
116 dIdt = R;
117 cumI += dIdt * input.dt;
118 return dIdt;
119 }
120
121 /* infiltration rate is less than rainfall intensity */
122 /* Newton-Raphson iteration to solve Eq. (6) for I */
123 /* guess new cumulative infiltration */
124 I = cumI + r * input.dt;
125 for (i = 0; i < MAX_ITER; i++) {
126 /* new infiltration rate */
127 r = -params.K0 / params.m * (C + I) / (1 - exp(I / params.m));
128 /* if new infiltration rate is greater than rainfall intensity,
129 * increase cumulative infiltration
130 */
131 if (r > R) {
132 f1 = I;
133 I = (I + f2) / 2.0;
134 df = I - f1;
135 }
136 /* otherwise, decrease cumulative infiltration */
137 else {
138 f2 = I;
139 I = (I + f1) / 2.0;
140 df = I - f2;
141 }
142 /* stop if cumulative infiltration converged */
143 if (fabs(df) <= TOLERANCE)
144 break;
145 }
146 if (i == MAX_ITER)
147 G_warning(
148 _("Maximum number of iterations exceeded at time step %d"),
149 timestep);
150
151 /* ponding time: tp will remain the same until next ponding occurs */
152 tp = t - input.dt + (I - cumI) / R;
153 /* if ponding time is greater than the current time,
154 * tp = t - dt + (I - cumI) / R > t
155 * (I - cumI) / R > dt
156 * I - cumI > R * dt
157 * means that additional infiltration (I - cumI) is greater than the
158 * total rainfall (R * dt), which cannot happen when there is no
159 * ponding, so infiltrate all rainfall
160 */
161 if (tp > t) {
162 /* full infiltration */
163 dIdt = R;
164 cumI += dIdt * input.dt;
165 return dIdt;
166 }
167
168 /* ponding starts if additional infiltration is less than the total
169 * rainfall because not all rainfall can be infiltrated in this time
170 * step
171 */
172 ponding = PONDING_STARTED;
173 }
174
175 /* if ponding just started */
176 if (ponding == PONDING_STARTED) {
177 /* lambda will remain the same until next ponding occurs */
178 lambda = 0.0;
179 fact = 1;
180 IC = I + C;
181 for (j = 1; j <= NUM_TERMS; j++) {
182 fact *= j;
183 lambda += pow(IC / params.m, (double)j) / (double)(j * fact);
184 }
185 /* lambda in Eq. (8) */
186 lambda = log(IC) - (log(IC) + lambda) / exp(C / params.m);
187 I += R * (t - tp) / 2.0;
188 }
189
190 /* Newton-Raphson iteration to solve Eq. (8) for I */
191 for (i = 0; i < MAX_ITER; i++) {
192 IC = I + C;
193 sum = 0.0;
194 fact = 1;
195 for (j = 1; j <= NUM_TERMS; j++) {
196 fact *= j;
197 sum += pow(IC / params.m, (double)j) / (double)(j * fact);
198 }
199 /* Eq. (8) - (t - tp) in hr: should converge to 0
200 * Note that sum is outside 1/exp(C/m) in Eq. (8), but inside in his
201 * TMOD9502.F. Based on lambda and his code, it looks like a typo in
202 * Eq. (8).
203 */
204 f1 = -(log(IC) - (log(IC) + sum) / exp(C / params.m) - lambda) /
205 (params.K0 / params.m) - (t - tp);
206 /* inverse of Eq. (7) in hr/m */
207 f2 = (exp(I / params.m) - 1.0) / (IC * params.K0 / params.m);
208 /* -(Eq. (8) - (t-tp)) * Eq. (7): cumulative infiltration in a short
209 * time period
210 */
211 df = -f1 / f2;
212 I += df;
213 if (fabs(df) <= TOLERANCE)
214 break;
215 }
216 if (i == MAX_ITER)
217 G_warning(_("Maximum number of iterations exceeded at time step %d"),
218 timestep);
219
220 /* if new cumulative infiltration is less than the previous cumulative
221 * infiltration plus the total rainfall, update the current cumulative
222 * infiltration and guess cumulative infiltration for the next time step
223 */
224 if (I < cumI + R * input.dt) {
225 /* less than full infiltration */
226 dIdt = (I - cumI) / input.dt;
227 cumI = I;
228 /* initial guess for next time step */
229 I += dIdt * input.dt;
230 ponding = PONDING_YES;
231 } else {
232 /* full infiltration */
233 dIdt = R;
234 cumI += dIdt * input.dt;
235 ponding = PONDING_NO;
236 }
237
238 return dIdt;
239 }
240