- Email: [email protected]
On the discretization and cost-effectiveness of a finite element solution for hillslope subsurface flow
Journal of Hydrology, 110 (1989) 165-179
165
Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands
[3]
ON THE DISCRETIZATION AND COST-EFFECTIVENESS OF A FINITE ELEMENT SOLUTION FOR HILLSLOPE SUBSURFACE FLOW
A. CALVER and W.L. WOOD
Institute of Hydrology, Wallingford, Oxon OXIO 8BB (U.K.) Department of Mathematics, University of Reading, P.O. Box 220, Reading RG6 2AX (U.K.) (Received June 3, 1988; accepted after revision October 13, 1988)
ABSTRACT Calver, A. and Wood, W.L., 1989. On the discretization and cost-effectiveness of a finite element solution for hillslope subsurface flow. J. Hydrol., 110: 165-179. An investigation is carried out into the effects of the size of time and space steps on the discharge and hydraulic potential predictions of a finite element rainfall-runoff model, the Institute of Hydrology Distributed Model version 4. Accuracy of hydrological predictions is viewed in the light of computer processing times and recommendations are given for suitable discretiza. tions for this type of model.
INTRODUCTION
The model discussed in this paper is the Institute of Hydrology Distributed Model version 4 (IHDM4), a rainfall-runoff model based on the physical processes of runoff generation within a catchment. The development of this model is described by Beven et al. (1987) and an application to an Institute of Hydrology instrumented catchment is described by Calver (in press). The Richards equation for flow in porous media (Richards, 1931) is solved by numerical approximation, since practical conditions make it too complex for an analytical solution. A finite element method, the Galerkin method of weighted residuals, is used for the two spatial dimensions and an implicit finite difference scheme for the time dimension. Error in this type of modelling can arise from a number of sources, in particular from uncertainty in the precise value of some of the physical parameters, from the degree of appropriateness-of the mathematical model itself to the physical processes of the real world, and from the numerical approx" nation of the governing equations and boundary conditions. An important aspect of the error from the numerical approximation, which is the subject of the investigation in this paper, is the effect of the size of time and space steps used in a model run. Such an investigation also allows comment on the cost-effectiveness of different time and space discretizations. Because of the 0022-1694/89/$03.50
© 1989 Elsevier Science Publishers B.V.
166
complexity of the problem this investigation can only be carried out empirically. MODELLING METHODOLOGY
In the IHDM4, the subsurface solution is carried out for two-dimensional vertical sections of the catchment hillslopes which can if necessary incorporate variation in the third spatial dimension, slope width. Sufficient such vertical planes are used, and their flows appropriately combined, to represent the three-dimensior.al topography of the catchment. Saturated and unsaturated flows are dealt with together by the Richards equation and, at the downslope end of the hillslopes, outflow can take place from the soil/permeable material into a stream channel. (As the model was developed to cover a wide range of physiographic conditions it also offers a kinematic wave solution for any infiltration-excess or saturation-excess overland flow, though subsurface flow alone is the subject of this investigation.) The Richards equation may be expressed in the form:
o.i in a region V of (x, z) space, where Ox is the horizontal axis, Oz is the vertical axis, 0 is the soil moisture content by volume, ~ = ~ + z is the total hydraulic potential where ~ is the pressure potential, Kz, Kz are the hydraulic conductivity components in the Ox, Oz directions respectively, B is the variable slope width, Q, is the source or sink term and t is the time. The IHDM4 uses a @based formulation of the problem which assmnes that there is a relationship between 0 and ~ which is locally differentiable so that it is possible to put C(~b) = d0/d~, the specific moisture capacity of the soil. Equation (1) can then be rewritten as an equ~,tion for ~(x, z, t):
Ot
Ox
-~-~ B K , ( ¢ )
= Q.
(2)
Initial and boundary conditions are given, i.e. ~(x, z, 0), and values on the boundary, S, of the region V of either the hydraulic potential ~ or the velocity normal to the boundary: IKz~,K~]'n
ffi K V ~ - n
(3)
expressed in shorthand form, where n is the unit vector in the direction of the outward normal to S. If eqn. (2) is written as F (~b) = 0 then:
~v~ W(x, z) F(~) d V = 0
(4)
167
is also true, where W(x, z) is some convenient weighting or test function. Integrating eqn. (4) by parts gives:
Ot
~
BK, -~ + ~
BK~
dV
-IWKV,.nd8--;J[Q,]WdV
(5)
= {~,
Equation (5) is a weak form of the original differential eqn. (2); it is weakened by only being satisfied in a weighted integral or averaged sense. The region V is now divided into quadrilateral finite elements with nodes aligned in vertical columns such that, in this case, the elements are trapezia (Fig. 1). A finite element space discretization is made using piecewise bilinear basis functions N~ (x, z) such that N/(xs, zs) = ~ij (the Kronecker 6) and Ns (x, z) - 0 outside the "neighbourhood" of node j, i.e. the finite elements of which node j is a vertex. Equation (5) is approximated by substituting: m
q' =
Z
z)
jr1
where the j -- 1 , 2 , . . . , m refer to the m nodal points of the finite element grid, and W(x, z) = N/(x, z) in turn to produce m equations: dO
P - ~ + GO
F(t)
=
(6)
where @ is the vector of nodal values of ~, F(t) is the vector of boundard conditions and any source or sink terms. The matrices P arm G are the mass and stiffness matrix respectively given by their components: i I
pR.F,.CI p i ~I.ATION
I i
fiu~
;
i~i.'rm'ior,
i
!
DRAINAGE. i-DIVIDE,
Liq'-I °'~" _~
.~m
J-Y
I
bound,:~rms IMPERMF.A~I.E
,STR.EAM / i CHANNEL. ~ t C [ m_..~
BE,DROCK. 500 m
, pc~strions of nodes common 1"ooil simulations
Fig. I. Diagrammatic representation of hillslope flow domain.
168
Pi~
- ~J BCN~Nj
go =
dxdz
(7)
BK~ Ox Ox + BKz Oz Oz
The time derivative in eqn. (6) is approximated by the backward difference formula as:
where P is an averaged value of F(t) and ~ an average of ~ over the time interval. Equation (9) is then written as: AOt+At
=
(10)
H
The terms in the matrix A and the vector H depend on the value of the unknown ~t+at, hence the nonlinear eqns. (10) are solved by iteration. In the model the time step can change, so that if t~+l = tn + Atn and ~tn = ~ ' then iteration on all time stepv v~cept the first is started by the extrapolation: At~ ~n+l
--
Cn "}" 0.5 ~
(~n -- C n - 1 )
(11)
The iteration is then accelerated by an adaption of the Cooley method (Cooley, 1983) which works very efficiently for this problem. The equations are arranged so that they are solved for the change in the value of the potential ~n+l from one iteration to the next; i.e. for: 'k ~(Pn+l
k
"--
¢n+1 - -
j~k-I
~,,+i
(12)
The Cooley algorithm depends on the assumption that these changes will oscillate in sign from one iteration to the next so that the convergence of the iteration will be accelerated by damping this oscillation, i.e. by underrelaxation. If Akis the maximum change in nodal potentials at the kth iteration then the program uses: (a)s
= 1, f o r k
= 0
s = A~+l/(Co,~),fork > 0
(13)
where oJk is the Cooley relaxation parameter from the previous iteration; (b)& = (3 + s)/(3 + Isl),whens > otherwise & = o.5/Isl
-1 (14)
(c) if &lA~+1[ is greater than e,~, the specified maximum change in the potential allowed in an iteration, then: (Dk+l = emax/l~+ll otherwise
cok+1
--"
(~
(15)
169
The iteration is stopped when [Ak+~[is less than a prescribed tolerance. The flow across the boundary in the neighbourhood of the Dirichlet type nodes (i.e.where the potential is known) is calculated by the method of Lynch (1984). The equations for which there is a Dirichlet boundary condition, and which are consequently not used in the solution for the new potential vector O.+x, are held in reserve and the imbalance on each of these, when the new solution is inserted, gives the flow outwards across the corresponding part of the boundary by the following formula: At .=
i=i 5~i~i
(16)
+ fj
where node~ is a node on the Dirichlet boundary. There is an error in fbrn~ula (16) due to the approximatio~ of the time derivative. Each part of the first term in formula (16) is equal to:
_~tpji[~n+l 1
4~7] = +
1
0
,,+~
2 0 t 2 (qbi)n+l + . . At d~ ] 2 0f (~i)"+1 + higher order terms (17)
Hence the error in the estimate of the flow, for reasonably small values of the time step At, is proport~onai to the second time derivative of the potential solution ~. For the IHDM4 simulations reported here a single realistic but physiographically simple hillslope is used, with ty~ical, non-extreme values of geomorphological and hydrological properties. The hillslope is of ~ spatial scale appropriate to the small basin studies for which the IHDM4 is primarily designed. The hillslope is of 500m horizontal distance between divide and channel, and of a constant 0.1 slope with straight contours (Fig. 1). A permeable soil of 1.5 m depth overlies an impermeable bedrock. (It is frequently the case in the use of the I H D M that the overall domain is, because of the zone producing runoff, much greater in its horizontal than its vertical dimension.) The soil material considered here is spatially uniform with a saturated hydraulic conductivity, vertically and horizontally, of 0.1 m h -I and saturated water content of 40% by volume. The rehtionship of unsaturated .moisture content to moisture potential and the relationship of unsaturated hydraulic conductivity to moisture potential are represented by power curve functions ~tppropriate to soil material of medium or loam texture (Clapp and Hornberger, L978). The initial soil moisture potential is set at - 0.5 m water over the entire ~!ope. The tolerance level of hydrau !ic potential prediction between iterations is se~,at 0.001m water. Rainfall inr~ut to the hillslope is also simplified. The simulation covers a period of 25 h without rainfall followed by 5 h with spatial!y
170
uniform rainfall of 4 mm h-1 followed by a further 25 h without rainfall. Under these conditions of soil, topography and precipitation the hillslope flow response is entirely subsurface in nature. Simulation results are discussed for a range of calculation grids altering, within the fixed flow domain, the absolute size of the elements and the relationship of horizontal dimension, Ax, to vertical dimension, Az, of the elements. The grid remains static during the course of an individual run. In general one expects to place more finite element nodes where most change is expected, but degree of change can be similar throughout the domain for this physical problem. Table i shows the precise meshes used (together with aspects of the results reported below). Eleven of the grids are of regular geometry; the other two have elements of different sizes with finer divisions towards the flux boundaries. There are 24 nodal locations common to all these grids (Fig. 1), namely those at horizontal intervals of 100 m and vertical intervals of 0.5 m. A baseline grid is used as the basis of comparison of results. This is the grid of runs la and lb in which both the horizontal and vertical dimension of the elements is 0.5 m. General finite element experience suggests that this fine mesh will give the most accurate results. Simulations are discussed for fixed solution time steps, At, of 0.5 h and 0.1 h. RESULTS The outcome of the runoff simulation in terms of moisture potential and hydraulic potential is as follows tbr the baseline grid (for both 0.1 and 0.5 h time steps). After the first 15 h of drainage without rainfall, a saturated layer of about 15 cm has developed at the base of the soil profile and the unsaturated zone above this is increasingly drier towards the surface. After 25 h drainage a similar overall pattern prevails but with a slight increase in the vertical extent of the s ~ r ~ . s d zone and slightly drier surface horizons. After 30 h (with rainfall occurring between 25 and 30 h) the basal saturated layer has increased to about 20 cm and the unsaturated zone is wetter than at the start of rainfall, particularly near the ground surface. Two-and-a-half hours later (at about the time of peak discharge from the slope, see below) basal saturation of 27 cm occurs over most of the hillslope. At 1.0 and 0.5 m depth the soil is wetter than at the end of the rainfall but the surface has begun to dry. Fifteen hours into the drainage period following the end of rainfall, the thickness of the saturated layer is about 30 cm. Soil at 1.0 m depth is still becoming wetter, that at 0.5 m showing little change and that at the surface continuing to become drier. At the end of the simulation, 25 h after the end of storm rainfall, little further change has occurred except for the slight downward relocation of water in the soil profile. At the top of the hillslope the base of the soil profile does not show the general saturation described above because of the presence of the no-flow boundary at the divide and the predominance of drainage away from this zone. At the downslope end of the hillslcpe where seepage out of the soil occurs from the saturated zone, the unsaturated-saturated boundary shows a drainage
0.5 100 20 50 25 10 10 5 2 5 2~5 10~), 10' variable 2 0.5 100 20 50 25 I0 10 5 2 5 2.5 100, 10~ variable
la 2a 3a 4a 5a 6a 7a 8~t 9~ lOa lla !2a 13a lb 2b 3b 4b 5b 6b 7b 8b 9b 10b llb ~2b 13b
1 200 200 100 100 100 20 20 20 10 10 200, 20 1 200 200 100 100 100 20 20 20 10 10 200,20
0.5 0.5 O.1 0.5 0.25 0.1 0.5 0.25 0.1 0.5 0.25 0.5
Az/Az
0.5 0.5 0.1 0.5 0.25 0.1 0.5 0.25 0.1 0.5 0.25 0.5
Vertical dimension of elements Az (m) 0.5 0.5 0.5 0.5 05 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.1 0.1 O.1 0.1 0.1 0.1 0.1 0.t 0,1 0.1 0.1 0.1 0.1
Subsurface calculation time step At (h) 4004 24 416 44 147 816 204 707 4016 404 1407 c~O ;~2 40~}4 24 4 i6 t4 147 816 204 707 40!6 404 149~ 60 342
Total no. of nodes 1407 9 258 15 @~ 523 7'7 31,4 2327 145 629 20 2~5 51~];2 34 85~ 57 216 1718 246 1017 8182 533 2160 78 730
C.p.u, time (s)
253 230 252 229 249 253 252 251 25~ 254 252 226 226 853 841 838, 83~ 840 8~S 831 857 867 839 868 830 811
Total no. of iterations
1Lower fifth of hil]slope has the finer mesh. "~'Depths (m) of node rows from surface: 0,00 0.05 0.10 0.15 0.20 0.30 ~.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50. Distances (m) of node columns from di-,ide: 0 100 200 3 ~ 400 4~,) 420 43~ 440 450 460 470 480 490 492 494 496 498 500.
Horizontal dimension of elements Ax (m)
Run
Finite element grid specifications and aspects of simulation results
TABLE 1
Mass balance
-
0.02050 0.02284 0°02671 0~02064 0.02763 0.02516
0.00530 0.00053 0.01107 0.00161 0.01033 0.00945 0.00610 0.00869 0.00969 0.00467 0.00920 0.00354 0.00442 0.02501 0.01552 0.02433 0.02219 (}.02136 0.02454
(m
k=i
172 profile curve, in this case with drawdown extending a horizontal distance of about 20m from the seepage face. In this run (and all others) only the lowest of the nodes on the seepage face becomes saturated and it remains so throughout the course of the run. The overall result of the storm is to promote the saturation of the base of the soil profile, produce an outflow hydrograph (see below) from the seepage boundary at the foot of the hillslope and to increase overall soil moisture content over the hillslope. The sign-~i~cance from the hydrological point of view of the prediction of hydraulic and mois ure potential lies not only in their influence on flow paths and seepage discharge from the hillslope but also in plant water uptake, infiltration capacity and the likelihood of surface runoff, and the behaviour of chemical transfers. An overall measure of the difference in hydraulic potential va~ue, ~, between a given simulation run and the baseline run for a given time is defined as: (18) where the ~bb are the baseline run values of hydraulic potential and Nn is the number of node positions common to all runs. The build-up of this ~ deviation is seen to occur at an approximately constant rate during the course of each simulation. Figure 2 shows that for the 0.5 h time step runs the ~ deviation, determined for the end of the simulation, decreases with total number of nodes in the run and, further, that for ~he regular meshes (all but runs 12a and 13a) the ~ deviation is reduced, given the total number of nodes, for smaller values of Ax/Az. The rate of lowering of ~bdeviation decreases with total node number and this rate, for a given total node number, is less for Ax/Az values of 10 and 20 than for values of 100 and 200. The two non-regular meshes show comparatively h~gb ¢ deviations from the baseline run for their total node numbers. These patterns of deviation of ~b with node number are almost identical to those obtained for the simulations with a solute.on time step of 0.1 h (runs lb-13b): the ~ deviation bet~l,~en the baseline grid with a 0.5 h and a 0.1 h time step is 1.9 x 10-3m of wate~ at the end of the simulation. Inspection of the pressure potentials, ~, at all nodes shows for some runs small oscillations in part of the body of the hillslope and ~lightly larger such effects nearer the seepage face. In the case of run 5a wit ~ Ax[Az of 100, for example, oscillations in ¢ ~f a fraction of a millimetre of water occur between nodal values at a given depth at i00-200m from the channel and oscillations of a few centimetres occur 0-100m from the channel. Degree of oscillation is reduced for finer discretizations and is absent from the baseline simulation (to 10-5m water). This result is closely related to the well-known ~ i d Peclet number effect (see e.g., Huyahorn and Pinder, 1983). Pressure potential oscilh~ions produce,minor changes in magnitude of gradient of hydraulic potential, ~b; they do not change the direction of the gradient. Figure 3 shows the hillslope discharge calculated at 0.5 h intervals during
173 0.20 -
¢ deviation (m) \ e Zo e 1?.a
ei3(J
0.~5
O.tO
005
0 10
I fO0
I fOoo f o t a l number of nodes
., tO000
Fig. :~.Relationship of deviation in hydraulic potential, ~, from baseline run (~ee text for detail of derivation) and number of nodes for different calculation grids (~s defi:~e," i,~ Table 1).
the course of the simulations using the ~ soluLions from the 0.5 h time steps. The hydrographs calculated at 0.5 h intervals from the 0.1 b ~ solutions are virtually identical for a given grid. For the baseline run, la, the storm rainfall promotes a low hydrograph response, peaking at 9.4 x 10-3m3h -1 per metre width of hillslope four hours after the e~d of rainfall and increasing baseflow from 5.0 x 10 -3 m 3h- l m - 1prior to the sterm to 7.4 x 10 -3 m 3h- 1m- 1,25 h after the end of rainfall. The degree of devi~,tion from this situation using other calculation grids is apparent from Fig. 3 with coarser grids predicting a more extreme hydrograph response, though of similar timing, compared with the baseline run. The coarsest grids promote a small oscillatioii in outflow at the beginning of occurrence of precipitation. Such spurious oscillation is related to occurrences of the greatest os :illat~o,s in mag~litude of hydraulic gradient noted above. It is clear from Fig. 3 that there is a "running-~n" period at the beginning of the simul~_t~on be~'ore a ~teady basefiow is achieved, and that this settling down period is longer and shows mcro ~a~:~stion in flow values for grids with large elements near the base of the seepage face. Experiments with a
174 ~ 0.074
~ 0.09Z
0.06
0.05
rainf~I1 0.04 discherge m'~h .4m-I
I I
zQ
0.0~
\
\ zQ
002
0.01
to -tla, t3o
9a
t3a o
0
'-"
,
5
I
10
~
1
,
20
,
;)5 h
-
,
30
,
35
,
40
,~
,
50
,
5~
Fig. 3. Discharges per metre hillslope width at the foot of the hillslope during course of simulations using different grids.
second storm, of the same magnitude and duration as the first, occurring after 55 h from the start of ~imulation, yield similar findings to those of Fig. 3. The results are not identical because of the slightly higher soil moisture contents prevailing on ~he hillslope at the beginning of the second rainfall event. A measure ~)f the deviation in discharge, Q, from the baseline run:
where Qb indicates baseline run discharge values a M Nt the number of 0.5 h time steps, is shown in Fig. 4 in relation to the total number of calculation nodes in a run for the 0.5 h time step simulations. A<~ in the case of the deviation, this Q deviation for regular meshes decreases with increasing number of nodes and, for a given number of nodes, is less for lower Ax/Az values. In the case of the non-regular grids the Q deviation is relatively low for a given total nod~ number because the nodal values used for calculating discharge are those at or connected ~:o the seepage face which are withil, the region where the irregular grids are finer and where we expect the solution to be more accurate. Following from the qualitative comment on hydrograph form above, it will be apparent that a situation very similar to that of Fig. 4 applies to the 0.1 h ~ calculation time step simulations. The Q deviation between th,v baselino grid runs with 0.5 h and 0.1 h ~ime step is 9.0 x 10 -5 m3h -~.
175 0.t Q deviation (m3 h -~)
,,~2a
0.0t Ax/Az = ~~~---',~,,~Ax~ 5Oe
~ z : ZOO
e~J
0.00t
® 15a
0.000t t0
I 100
J 1000 total number of nodes
J 1OOO0
Fig. 4. Rdat~.vn~b~? ,~f deviation in discharge Q from baseline run (see text for details)and number of nodes for dii~fere~'~.tgrids.
Further asp~ct:;of the results of the simulations are now considered, namely the required nr~n!)er of iterations for solution and the mass balances achieved. The iteration r~ethod using the Cooley acceleration parameter is shown to be very efficient{Table 1). With the specified tolerance in potential of 0.001 m water, an averag ~ of between 2 and 3 iterations per 0.5 h ;tep is required and an average of 1-2 iterations per 0.1 h step. Iterations per tlme step are greatest in the initial Q~ett '" |ing down period and after the introduction of rainfall. For a given solution tkne step, the total number of iterations differs littlebetween grids. ][terqti~nnumber reqaired for solution is therefore not a basis for making a choice between grid geometries. Table I shows the mass balance error as calculated by the Lynch method for the various grids tested. It is interesting to see from this that the Lynch result is of no use for indicating which grid is bes~ in a time-dependent problem. It gives an apparently very good result for run 2a which, together with runs 2b and 2c (another run with the same grid ancl time step At = 0.05 h), snows a particularly high peak in the flow prediction in the run-in period before the rain starts. The error in the calculation of the mass balance is proportional to ~2~/~t2 as shown in eqn. (17) because of the approximation of the time derivative. ~Tow the fixed head node at the base of the seepage face remains at
176
the same fixed level throughout a run. Hence the change in outflow at this fixed head boundary node corresponds to the change in the values of ~ in the neighbourhood of this node. The spurious peak in the outflow in runs 2a, 2b and 2c cerresponds to a spurious peak in the ~ values for this grid. The error in runs 2a, 2b and 2c is proportional to the second time derivative of the ~ solution as At tends to zero but the grid remains fixed. Hence it is possible to deduce that the cumulative mass balance error give~ by the Lynch formula is the result of adding terms with changing sign of ~2~ / ~ t 2 and this makes it speciously lower than the results obtained with other runs using other grids where the second time derivative be,naves correctly. It is conjectured that the greater mass balance error associated with the 0.1 h time step is a function of a cumulative effect of performh~g the non-linear solution over the greater number of the smaller time steps. Consideration is now given to the central processor unit (c.p.u.) time involved in achieving the above results. Though not of primary hydrological significance per se, this may influencv the choice of simulation details when practical aspects of achieving a prediction are taken into account. For mainframe computing the processing time and therefore cost may be an important consideration in running simulations. For microcomputers (of sufficient storage capacity) run time may be primarily a matter of convenience. Such considerations are of particular relevance given that in many applications it is svund practice to calibrate the model for some of the physical parameter values involved rather than conduct a single run based on a priori best estimates. Figure 5 shows the computer processing time used in running' the simulations on an IBM mainframe (4381 model group 3 dual processor). Run ~ime increases approximately linearly with total node number for the 0.5 h time step simulations and for the 0.1 h runs, the latter group using three to four times the processing time of the fondler tb~ ~he same grid. There is not a large c.p.u. "overhead" independent of the number of calculation nodes. For equal time step and number of nodes (or approximately so, given the data set here) a lower value of Ax/Az reduces run time (cf. runs la and 9a, lb and 9b, 10a and 3a, 10b and 3b). Since smaller tlme steps involve a greater number of loops through the model but may be associated with fewer iterations per time step to reach a solution, it is not apparent until runs are undertaken whether there is some optimum value of Ai between 0.1 and 0.5 h for which a m i n i m u m run th-~e is achieved. This however is not the case and runs using, for example, At = 0.25h yield run times between those of the At -- 0.1 and At = 0.5 runs. For grid 10 and At = 0.25, run time is 276s and for grid 13 it is 359s. The model is primarily designed for use over the period of a storm hydrograph and it is inappropriate in that context to consider calculation time steps m u c h in excess of 0.5 h because of ti~e expected rate of change o~ ~oil moisture potentials and discharges. The results with a time step of < 0.5 h differ very littlefrom those of the 0.5 h steps, and the former use considerably greater processing times than the latter,suggesting that the 0.5 h step is an appropriate and efficientchoice for this type of problem. [It is a subject of further research
177 2500
~,b i 8187. s
tb ~ s~s2
9a • ,1tb
20oo C..p.U. time
(s)
*6b
1500
la,
,8b
lOOO
* t3b ,110 *1Oh
OUU
5b*
?.b za~
0
to
4bo 4~,
®7b
,3a *t3a * lOa
e60
,8o
,JZb 50® ®Ta ®l?.a , I t00 t000 total nu.mber of nodes
I
10000
Fig. 5o Relationship of computer processing time to number of nodes in grid.
to investigatv the use of such a model for baseflow, in which case lon~er time steps could be considered.] In everyday use of the program the subsurface calculation time step can be flexible depending on the number of iterations necessary for solution: an initial and 8 rn,~ximum subsurface calculation time step are specified. With practical cons.~.~lerationsin mind, two measures of discharge predictive success are plotted against processing time for the 0.5 h time step runs in Fig. 6. The discharge errors relate, again because of prac~:i~alconsiderations, to the period after mo&~l run-in effects can be expected to have dispersed. To exclude the run-in effects of all the runs in this comparative work, discharge predictions are considered after 25 h of simulation, that is,from the beginning of the rainstorm. The two measures plotted are the percentage of half-hourly
178 t.0 baseline run
ul (i) =3
erie A
t0a~e :t3a e
[3
Z0
8a
e70
. °0.~ t=
A
:ga
o
U
6a
(..
&
o
o.61_%
g,o
e
3a
® within ~ 10~ of |~!~oseof baseline run
&
A
,,
+ - 5~
.
.
.
.
.
.
.
.
.
.
e
Q,
0.41 ~a 0
I
500
I
I
t000 t500 c.p.u.fime (s)
I
Z000
I
ZS00
Fig. 6. Relationship of processing time to accuracy of discharge predictions.
v~.lues of discharge within _+10% and _+5% of the values of the baseline run. It should perhaps be noted at this stage that in practical use the mode] is likely to be calibrated against discharge. F i b r e 6 suggests that, in terms of cost-effectiveness, simulations 7a, 8a, 10a, 12a and 13a are the best at the + 10% level, offering reasonably accurate predictions without excessive run times. These runs are those with regular meshes of Ax/Az not greater than 20, or of non-regular mesh as used here, and with total node numbers less than about 750 (which averages about one node per square metre of the hillslope section in th~se cases). The runs ~uich perform best at the _+5% prediction level are runs 8a, 10a and 13a. They have Ax/Az values of 10, 20 and variable bat, of the five runs above, have the grep~test numbers of nodes, that is, the generally smaller elvments. CONCLUSIONS
Given the type of rainfall-runoff model methodology and the nature of the subsurface storm~ow response considered here, certain guidvlines car be given Zor goo~l hydrological prediction and cost-effective modelling. These points are as follows. (1) With regard to the overall prediction of hydraulic potential values over the hillslope and the reduction of oscillation in the ~olution, generally smaller elements are favoured with a ratio of horizontal to vertical element dimension of ~ 20 for regular grids, ideally ~<10. (2) With xegard to hillslope discharge prediction over the course of the hydrograph, it is again wise to have generally smaller elements ~nd Ax/Az ~ 20 for regular grids. It is also acceptable to have an i~egular grid (see details above) with fining of the mesh towards the seepage face and the infiltration boundary.
179
(3) In addition to accuracy of prediction, finer grids have the advantage of reducing the necessary run-in time at the beginning of ~imuiation. Lower Ax/Az values for a given total number of nodes reduce run times. (4) With this methodology neither number of iterations required nor cumulative mass balance serve as a distinguishing factor for choice between different calculation meshes. (5) If At = 0.5 h is the maximum time step considered for the solution time step for the response of soil conditions to storm rainfall, this value is seen to serve efficiently, giving little if any loss of accuracy over shorter time steps and using shorter c.p.u, times. (6) There is a conflict between the better discharge predictions of smaller elements and the associated longer c.p.u, times. These experiments suggest that with the meshes recomme~ded above (a regular mesh of Ax/Az ~ 20 or the irregular meshes as used here) a higher number of nodes (here an average of one node per 1-2 m ~ of domain area) should be used if the emphasis is on more accurate prediction and a lower number of nodes (here one node per 2-10 m 2) if shorter execution time is the priority. It is appreciated that in particular applications physical inhomogeneities may modify these mesh considerations. The above guidelines, used within the context of a particular case, should ensure that inaccuracies of numerical approximation as influenced by time and space steps are reduced as far as practically possible. ACKNOWLEDGEMENT
This work has been funded by the U.K. Ministry of Agriculture, Fisheries and Food. REFERI~INCES
Beven, K., Calver, A. and Morris, E.M., 1987. The Institute of Hydrology Distributed Model. ?ep. No. 98, Institute of Hydrology, Wallingford. Calver, A., in press. Calibration, sensitivity and validation of a physically-based rainfall-nmoff model. J. Hydrol., 103: 103-115. Clapp, R.B. and Hornberger, G.M., 1978. Empirical equations for some soil hydraulic pror~rties. Water Re,our. Rcs., 14: 601--604. Coole,v, R.L., 1983. Some new procedure~ for numerical ~olution of variably saturated flow problems. Water Resour. Res., 19: 1271-1285. Huyakorn, P.S. and Pinder, G.F., 1983. Computational Methods in Subsurface Flow. Academic Press, New York, N.Y., 473 pp. Lynch, D.R., 1984. Mass conservation in finite element groundwat~.r models. Adv. Water Resour., 7: 67-75. Richards, L.A., 1931. Capillary conduction of liquids through ~)orous mediums. Physics, 1: 318-333.
165
Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands
[3]
ON THE DISCRETIZATION AND COST-EFFECTIVENESS OF A FINITE ELEMENT SOLUTION FOR HILLSLOPE SUBSURFACE FLOW
A. CALVER and W.L. WOOD
Institute of Hydrology, Wallingford, Oxon OXIO 8BB (U.K.) Department of Mathematics, University of Reading, P.O. Box 220, Reading RG6 2AX (U.K.) (Received June 3, 1988; accepted after revision October 13, 1988)
ABSTRACT Calver, A. and Wood, W.L., 1989. On the discretization and cost-effectiveness of a finite element solution for hillslope subsurface flow. J. Hydrol., 110: 165-179. An investigation is carried out into the effects of the size of time and space steps on the discharge and hydraulic potential predictions of a finite element rainfall-runoff model, the Institute of Hydrology Distributed Model version 4. Accuracy of hydrological predictions is viewed in the light of computer processing times and recommendations are given for suitable discretiza. tions for this type of model.
INTRODUCTION
The model discussed in this paper is the Institute of Hydrology Distributed Model version 4 (IHDM4), a rainfall-runoff model based on the physical processes of runoff generation within a catchment. The development of this model is described by Beven et al. (1987) and an application to an Institute of Hydrology instrumented catchment is described by Calver (in press). The Richards equation for flow in porous media (Richards, 1931) is solved by numerical approximation, since practical conditions make it too complex for an analytical solution. A finite element method, the Galerkin method of weighted residuals, is used for the two spatial dimensions and an implicit finite difference scheme for the time dimension. Error in this type of modelling can arise from a number of sources, in particular from uncertainty in the precise value of some of the physical parameters, from the degree of appropriateness-of the mathematical model itself to the physical processes of the real world, and from the numerical approx" nation of the governing equations and boundary conditions. An important aspect of the error from the numerical approximation, which is the subject of the investigation in this paper, is the effect of the size of time and space steps used in a model run. Such an investigation also allows comment on the cost-effectiveness of different time and space discretizations. Because of the 0022-1694/89/$03.50
© 1989 Elsevier Science Publishers B.V.
166
complexity of the problem this investigation can only be carried out empirically. MODELLING METHODOLOGY
In the IHDM4, the subsurface solution is carried out for two-dimensional vertical sections of the catchment hillslopes which can if necessary incorporate variation in the third spatial dimension, slope width. Sufficient such vertical planes are used, and their flows appropriately combined, to represent the three-dimensior.al topography of the catchment. Saturated and unsaturated flows are dealt with together by the Richards equation and, at the downslope end of the hillslopes, outflow can take place from the soil/permeable material into a stream channel. (As the model was developed to cover a wide range of physiographic conditions it also offers a kinematic wave solution for any infiltration-excess or saturation-excess overland flow, though subsurface flow alone is the subject of this investigation.) The Richards equation may be expressed in the form:
o.i in a region V of (x, z) space, where Ox is the horizontal axis, Oz is the vertical axis, 0 is the soil moisture content by volume, ~ = ~ + z is the total hydraulic potential where ~ is the pressure potential, Kz, Kz are the hydraulic conductivity components in the Ox, Oz directions respectively, B is the variable slope width, Q, is the source or sink term and t is the time. The IHDM4 uses a @based formulation of the problem which assmnes that there is a relationship between 0 and ~ which is locally differentiable so that it is possible to put C(~b) = d0/d~, the specific moisture capacity of the soil. Equation (1) can then be rewritten as an equ~,tion for ~(x, z, t):
Ot
Ox
-~-~ B K , ( ¢ )
= Q.
(2)
Initial and boundary conditions are given, i.e. ~(x, z, 0), and values on the boundary, S, of the region V of either the hydraulic potential ~ or the velocity normal to the boundary: IKz~,K~]'n
ffi K V ~ - n
(3)
expressed in shorthand form, where n is the unit vector in the direction of the outward normal to S. If eqn. (2) is written as F (~b) = 0 then:
~v~ W(x, z) F(~) d V = 0
(4)
167
is also true, where W(x, z) is some convenient weighting or test function. Integrating eqn. (4) by parts gives:
Ot
~
BK, -~ + ~
BK~
dV
-IWKV,.nd8--;J[Q,]WdV
(5)
= {~,
Equation (5) is a weak form of the original differential eqn. (2); it is weakened by only being satisfied in a weighted integral or averaged sense. The region V is now divided into quadrilateral finite elements with nodes aligned in vertical columns such that, in this case, the elements are trapezia (Fig. 1). A finite element space discretization is made using piecewise bilinear basis functions N~ (x, z) such that N/(xs, zs) = ~ij (the Kronecker 6) and Ns (x, z) - 0 outside the "neighbourhood" of node j, i.e. the finite elements of which node j is a vertex. Equation (5) is approximated by substituting: m
q' =
Z
z)
jr1
where the j -- 1 , 2 , . . . , m refer to the m nodal points of the finite element grid, and W(x, z) = N/(x, z) in turn to produce m equations: dO
P - ~ + GO
F(t)
=
(6)
where @ is the vector of nodal values of ~, F(t) is the vector of boundard conditions and any source or sink terms. The matrices P arm G are the mass and stiffness matrix respectively given by their components: i I
pR.F,.CI p i ~I.ATION
I i
fiu~
;
i~i.'rm'ior,
i
!
DRAINAGE. i-DIVIDE,
Liq'-I °'~" _~
.~m
J-Y
I
bound,:~rms IMPERMF.A~I.E
,STR.EAM / i CHANNEL. ~ t C [ m_..~
BE,DROCK. 500 m
, pc~strions of nodes common 1"ooil simulations
Fig. I. Diagrammatic representation of hillslope flow domain.
168
Pi~
- ~J BCN~Nj
go =
dxdz
(7)
BK~ Ox Ox + BKz Oz Oz
The time derivative in eqn. (6) is approximated by the backward difference formula as:
where P is an averaged value of F(t) and ~ an average of ~ over the time interval. Equation (9) is then written as: AOt+At
=
(10)
H
The terms in the matrix A and the vector H depend on the value of the unknown ~t+at, hence the nonlinear eqns. (10) are solved by iteration. In the model the time step can change, so that if t~+l = tn + Atn and ~tn = ~ ' then iteration on all time stepv v~cept the first is started by the extrapolation: At~ ~n+l
--
Cn "}" 0.5 ~
(~n -- C n - 1 )
(11)
The iteration is then accelerated by an adaption of the Cooley method (Cooley, 1983) which works very efficiently for this problem. The equations are arranged so that they are solved for the change in the value of the potential ~n+l from one iteration to the next; i.e. for: 'k ~(Pn+l
k
"--
¢n+1 - -
j~k-I
~,,+i
(12)
The Cooley algorithm depends on the assumption that these changes will oscillate in sign from one iteration to the next so that the convergence of the iteration will be accelerated by damping this oscillation, i.e. by underrelaxation. If Akis the maximum change in nodal potentials at the kth iteration then the program uses: (a)s
= 1, f o r k
= 0
s = A~+l/(Co,~),fork > 0
(13)
where oJk is the Cooley relaxation parameter from the previous iteration; (b)& = (3 + s)/(3 + Isl),whens > otherwise & = o.5/Isl
-1 (14)
(c) if &lA~+1[ is greater than e,~, the specified maximum change in the potential allowed in an iteration, then: (Dk+l = emax/l~+ll otherwise
cok+1
--"
(~
(15)
169
The iteration is stopped when [Ak+~[is less than a prescribed tolerance. The flow across the boundary in the neighbourhood of the Dirichlet type nodes (i.e.where the potential is known) is calculated by the method of Lynch (1984). The equations for which there is a Dirichlet boundary condition, and which are consequently not used in the solution for the new potential vector O.+x, are held in reserve and the imbalance on each of these, when the new solution is inserted, gives the flow outwards across the corresponding part of the boundary by the following formula: At .=
i=i 5~i~i
(16)
+ fj
where node~ is a node on the Dirichlet boundary. There is an error in fbrn~ula (16) due to the approximatio~ of the time derivative. Each part of the first term in formula (16) is equal to:
_~tpji[~n+l 1
4~7] = +
1
0
,,+~
2 0 t 2 (qbi)n+l + . . At d~ ] 2 0f (~i)"+1 + higher order terms (17)
Hence the error in the estimate of the flow, for reasonably small values of the time step At, is proport~onai to the second time derivative of the potential solution ~. For the IHDM4 simulations reported here a single realistic but physiographically simple hillslope is used, with ty~ical, non-extreme values of geomorphological and hydrological properties. The hillslope is of ~ spatial scale appropriate to the small basin studies for which the IHDM4 is primarily designed. The hillslope is of 500m horizontal distance between divide and channel, and of a constant 0.1 slope with straight contours (Fig. 1). A permeable soil of 1.5 m depth overlies an impermeable bedrock. (It is frequently the case in the use of the I H D M that the overall domain is, because of the zone producing runoff, much greater in its horizontal than its vertical dimension.) The soil material considered here is spatially uniform with a saturated hydraulic conductivity, vertically and horizontally, of 0.1 m h -I and saturated water content of 40% by volume. The rehtionship of unsaturated .moisture content to moisture potential and the relationship of unsaturated hydraulic conductivity to moisture potential are represented by power curve functions ~tppropriate to soil material of medium or loam texture (Clapp and Hornberger, L978). The initial soil moisture potential is set at - 0.5 m water over the entire ~!ope. The tolerance level of hydrau !ic potential prediction between iterations is se~,at 0.001m water. Rainfall inr~ut to the hillslope is also simplified. The simulation covers a period of 25 h without rainfall followed by 5 h with spatial!y
170
uniform rainfall of 4 mm h-1 followed by a further 25 h without rainfall. Under these conditions of soil, topography and precipitation the hillslope flow response is entirely subsurface in nature. Simulation results are discussed for a range of calculation grids altering, within the fixed flow domain, the absolute size of the elements and the relationship of horizontal dimension, Ax, to vertical dimension, Az, of the elements. The grid remains static during the course of an individual run. In general one expects to place more finite element nodes where most change is expected, but degree of change can be similar throughout the domain for this physical problem. Table i shows the precise meshes used (together with aspects of the results reported below). Eleven of the grids are of regular geometry; the other two have elements of different sizes with finer divisions towards the flux boundaries. There are 24 nodal locations common to all these grids (Fig. 1), namely those at horizontal intervals of 100 m and vertical intervals of 0.5 m. A baseline grid is used as the basis of comparison of results. This is the grid of runs la and lb in which both the horizontal and vertical dimension of the elements is 0.5 m. General finite element experience suggests that this fine mesh will give the most accurate results. Simulations are discussed for fixed solution time steps, At, of 0.5 h and 0.1 h. RESULTS The outcome of the runoff simulation in terms of moisture potential and hydraulic potential is as follows tbr the baseline grid (for both 0.1 and 0.5 h time steps). After the first 15 h of drainage without rainfall, a saturated layer of about 15 cm has developed at the base of the soil profile and the unsaturated zone above this is increasingly drier towards the surface. After 25 h drainage a similar overall pattern prevails but with a slight increase in the vertical extent of the s ~ r ~ . s d zone and slightly drier surface horizons. After 30 h (with rainfall occurring between 25 and 30 h) the basal saturated layer has increased to about 20 cm and the unsaturated zone is wetter than at the start of rainfall, particularly near the ground surface. Two-and-a-half hours later (at about the time of peak discharge from the slope, see below) basal saturation of 27 cm occurs over most of the hillslope. At 1.0 and 0.5 m depth the soil is wetter than at the end of the rainfall but the surface has begun to dry. Fifteen hours into the drainage period following the end of rainfall, the thickness of the saturated layer is about 30 cm. Soil at 1.0 m depth is still becoming wetter, that at 0.5 m showing little change and that at the surface continuing to become drier. At the end of the simulation, 25 h after the end of storm rainfall, little further change has occurred except for the slight downward relocation of water in the soil profile. At the top of the hillslope the base of the soil profile does not show the general saturation described above because of the presence of the no-flow boundary at the divide and the predominance of drainage away from this zone. At the downslope end of the hillslcpe where seepage out of the soil occurs from the saturated zone, the unsaturated-saturated boundary shows a drainage
0.5 100 20 50 25 10 10 5 2 5 2~5 10~), 10' variable 2 0.5 100 20 50 25 I0 10 5 2 5 2.5 100, 10~ variable
la 2a 3a 4a 5a 6a 7a 8~t 9~ lOa lla !2a 13a lb 2b 3b 4b 5b 6b 7b 8b 9b 10b llb ~2b 13b
1 200 200 100 100 100 20 20 20 10 10 200, 20 1 200 200 100 100 100 20 20 20 10 10 200,20
0.5 0.5 O.1 0.5 0.25 0.1 0.5 0.25 0.1 0.5 0.25 0.5
Az/Az
0.5 0.5 0.1 0.5 0.25 0.1 0.5 0.25 0.1 0.5 0.25 0.5
Vertical dimension of elements Az (m) 0.5 0.5 0.5 0.5 05 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.1 0.1 O.1 0.1 0.1 0.1 0.1 0.t 0,1 0.1 0.1 0.1 0.1
Subsurface calculation time step At (h) 4004 24 416 44 147 816 204 707 4016 404 1407 c~O ;~2 40~}4 24 4 i6 t4 147 816 204 707 40!6 404 149~ 60 342
Total no. of nodes 1407 9 258 15 @~ 523 7'7 31,4 2327 145 629 20 2~5 51~];2 34 85~ 57 216 1718 246 1017 8182 533 2160 78 730
C.p.u, time (s)
253 230 252 229 249 253 252 251 25~ 254 252 226 226 853 841 838, 83~ 840 8~S 831 857 867 839 868 830 811
Total no. of iterations
1Lower fifth of hil]slope has the finer mesh. "~'Depths (m) of node rows from surface: 0,00 0.05 0.10 0.15 0.20 0.30 ~.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50. Distances (m) of node columns from di-,ide: 0 100 200 3 ~ 400 4~,) 420 43~ 440 450 460 470 480 490 492 494 496 498 500.
Horizontal dimension of elements Ax (m)
Run
Finite element grid specifications and aspects of simulation results
TABLE 1
Mass balance
-
0.02050 0.02284 0°02671 0~02064 0.02763 0.02516
0.00530 0.00053 0.01107 0.00161 0.01033 0.00945 0.00610 0.00869 0.00969 0.00467 0.00920 0.00354 0.00442 0.02501 0.01552 0.02433 0.02219 (}.02136 0.02454
(m
k=i
172 profile curve, in this case with drawdown extending a horizontal distance of about 20m from the seepage face. In this run (and all others) only the lowest of the nodes on the seepage face becomes saturated and it remains so throughout the course of the run. The overall result of the storm is to promote the saturation of the base of the soil profile, produce an outflow hydrograph (see below) from the seepage boundary at the foot of the hillslope and to increase overall soil moisture content over the hillslope. The sign-~i~cance from the hydrological point of view of the prediction of hydraulic and mois ure potential lies not only in their influence on flow paths and seepage discharge from the hillslope but also in plant water uptake, infiltration capacity and the likelihood of surface runoff, and the behaviour of chemical transfers. An overall measure of the difference in hydraulic potential va~ue, ~, between a given simulation run and the baseline run for a given time is defined as: (18) where the ~bb are the baseline run values of hydraulic potential and Nn is the number of node positions common to all runs. The build-up of this ~ deviation is seen to occur at an approximately constant rate during the course of each simulation. Figure 2 shows that for the 0.5 h time step runs the ~ deviation, determined for the end of the simulation, decreases with total number of nodes in the run and, further, that for ~he regular meshes (all but runs 12a and 13a) the ~ deviation is reduced, given the total number of nodes, for smaller values of Ax/Az. The rate of lowering of ~bdeviation decreases with total node number and this rate, for a given total node number, is less for Ax/Az values of 10 and 20 than for values of 100 and 200. The two non-regular meshes show comparatively h~gb ¢ deviations from the baseline run for their total node numbers. These patterns of deviation of ~b with node number are almost identical to those obtained for the simulations with a solute.on time step of 0.1 h (runs lb-13b): the ~ deviation bet~l,~en the baseline grid with a 0.5 h and a 0.1 h time step is 1.9 x 10-3m of wate~ at the end of the simulation. Inspection of the pressure potentials, ~, at all nodes shows for some runs small oscillations in part of the body of the hillslope and ~lightly larger such effects nearer the seepage face. In the case of run 5a wit ~ Ax[Az of 100, for example, oscillations in ¢ ~f a fraction of a millimetre of water occur between nodal values at a given depth at i00-200m from the channel and oscillations of a few centimetres occur 0-100m from the channel. Degree of oscillation is reduced for finer discretizations and is absent from the baseline simulation (to 10-5m water). This result is closely related to the well-known ~ i d Peclet number effect (see e.g., Huyahorn and Pinder, 1983). Pressure potential oscilh~ions produce,minor changes in magnitude of gradient of hydraulic potential, ~b; they do not change the direction of the gradient. Figure 3 shows the hillslope discharge calculated at 0.5 h intervals during
173 0.20 -
¢ deviation (m) \ e Zo e 1?.a
ei3(J
0.~5
O.tO
005
0 10
I fO0
I fOoo f o t a l number of nodes
., tO000
Fig. :~.Relationship of deviation in hydraulic potential, ~, from baseline run (~ee text for detail of derivation) and number of nodes for different calculation grids (~s defi:~e," i,~ Table 1).
the course of the simulations using the ~ soluLions from the 0.5 h time steps. The hydrographs calculated at 0.5 h intervals from the 0.1 b ~ solutions are virtually identical for a given grid. For the baseline run, la, the storm rainfall promotes a low hydrograph response, peaking at 9.4 x 10-3m3h -1 per metre width of hillslope four hours after the e~d of rainfall and increasing baseflow from 5.0 x 10 -3 m 3h- l m - 1prior to the sterm to 7.4 x 10 -3 m 3h- 1m- 1,25 h after the end of rainfall. The degree of devi~,tion from this situation using other calculation grids is apparent from Fig. 3 with coarser grids predicting a more extreme hydrograph response, though of similar timing, compared with the baseline run. The coarsest grids promote a small oscillatioii in outflow at the beginning of occurrence of precipitation. Such spurious oscillation is related to occurrences of the greatest os :illat~o,s in mag~litude of hydraulic gradient noted above. It is clear from Fig. 3 that there is a "running-~n" period at the beginning of the simul~_t~on be~'ore a ~teady basefiow is achieved, and that this settling down period is longer and shows mcro ~a~:~stion in flow values for grids with large elements near the base of the seepage face. Experiments with a
174 ~ 0.074
~ 0.09Z
0.06
0.05
rainf~I1 0.04 discherge m'~h .4m-I
I I
zQ
0.0~
\
\ zQ
002
0.01
to -tla, t3o
9a
t3a o
0
'-"
,
5
I
10
~
1
,
20
,
;)5 h
-
,
30
,
35
,
40
,~
,
50
,
5~
Fig. 3. Discharges per metre hillslope width at the foot of the hillslope during course of simulations using different grids.
second storm, of the same magnitude and duration as the first, occurring after 55 h from the start of ~imulation, yield similar findings to those of Fig. 3. The results are not identical because of the slightly higher soil moisture contents prevailing on ~he hillslope at the beginning of the second rainfall event. A measure ~)f the deviation in discharge, Q, from the baseline run:
where Qb indicates baseline run discharge values a M Nt the number of 0.5 h time steps, is shown in Fig. 4 in relation to the total number of calculation nodes in a run for the 0.5 h time step simulations. A<~ in the case of the deviation, this Q deviation for regular meshes decreases with increasing number of nodes and, for a given number of nodes, is less for lower Ax/Az values. In the case of the non-regular grids the Q deviation is relatively low for a given total nod~ number because the nodal values used for calculating discharge are those at or connected ~:o the seepage face which are withil, the region where the irregular grids are finer and where we expect the solution to be more accurate. Following from the qualitative comment on hydrograph form above, it will be apparent that a situation very similar to that of Fig. 4 applies to the 0.1 h ~ calculation time step simulations. The Q deviation between th,v baselino grid runs with 0.5 h and 0.1 h ~ime step is 9.0 x 10 -5 m3h -~.
175 0.t Q deviation (m3 h -~)
,,~2a
0.0t Ax/Az = ~~~---',~,,~Ax~ 5Oe
~ z : ZOO
e~J
0.00t
® 15a
0.000t t0
I 100
J 1000 total number of nodes
J 1OOO0
Fig. 4. Rdat~.vn~b~? ,~f deviation in discharge Q from baseline run (see text for details)and number of nodes for dii~fere~'~.tgrids.
Further asp~ct:;of the results of the simulations are now considered, namely the required nr~n!)er of iterations for solution and the mass balances achieved. The iteration r~ethod using the Cooley acceleration parameter is shown to be very efficient{Table 1). With the specified tolerance in potential of 0.001 m water, an averag ~ of between 2 and 3 iterations per 0.5 h ;tep is required and an average of 1-2 iterations per 0.1 h step. Iterations per tlme step are greatest in the initial Q~ett '" |ing down period and after the introduction of rainfall. For a given solution tkne step, the total number of iterations differs littlebetween grids. ][terqti~nnumber reqaired for solution is therefore not a basis for making a choice between grid geometries. Table I shows the mass balance error as calculated by the Lynch method for the various grids tested. It is interesting to see from this that the Lynch result is of no use for indicating which grid is bes~ in a time-dependent problem. It gives an apparently very good result for run 2a which, together with runs 2b and 2c (another run with the same grid ancl time step At = 0.05 h), snows a particularly high peak in the flow prediction in the run-in period before the rain starts. The error in the calculation of the mass balance is proportional to ~2~/~t2 as shown in eqn. (17) because of the approximation of the time derivative. ~Tow the fixed head node at the base of the seepage face remains at
176
the same fixed level throughout a run. Hence the change in outflow at this fixed head boundary node corresponds to the change in the values of ~ in the neighbourhood of this node. The spurious peak in the outflow in runs 2a, 2b and 2c cerresponds to a spurious peak in the ~ values for this grid. The error in runs 2a, 2b and 2c is proportional to the second time derivative of the ~ solution as At tends to zero but the grid remains fixed. Hence it is possible to deduce that the cumulative mass balance error give~ by the Lynch formula is the result of adding terms with changing sign of ~2~ / ~ t 2 and this makes it speciously lower than the results obtained with other runs using other grids where the second time derivative be,naves correctly. It is conjectured that the greater mass balance error associated with the 0.1 h time step is a function of a cumulative effect of performh~g the non-linear solution over the greater number of the smaller time steps. Consideration is now given to the central processor unit (c.p.u.) time involved in achieving the above results. Though not of primary hydrological significance per se, this may influencv the choice of simulation details when practical aspects of achieving a prediction are taken into account. For mainframe computing the processing time and therefore cost may be an important consideration in running simulations. For microcomputers (of sufficient storage capacity) run time may be primarily a matter of convenience. Such considerations are of particular relevance given that in many applications it is svund practice to calibrate the model for some of the physical parameter values involved rather than conduct a single run based on a priori best estimates. Figure 5 shows the computer processing time used in running' the simulations on an IBM mainframe (4381 model group 3 dual processor). Run ~ime increases approximately linearly with total node number for the 0.5 h time step simulations and for the 0.1 h runs, the latter group using three to four times the processing time of the fondler tb~ ~he same grid. There is not a large c.p.u. "overhead" independent of the number of calculation nodes. For equal time step and number of nodes (or approximately so, given the data set here) a lower value of Ax/Az reduces run time (cf. runs la and 9a, lb and 9b, 10a and 3a, 10b and 3b). Since smaller tlme steps involve a greater number of loops through the model but may be associated with fewer iterations per time step to reach a solution, it is not apparent until runs are undertaken whether there is some optimum value of Ai between 0.1 and 0.5 h for which a m i n i m u m run th-~e is achieved. This however is not the case and runs using, for example, At = 0.25h yield run times between those of the At -- 0.1 and At = 0.5 runs. For grid 10 and At = 0.25, run time is 276s and for grid 13 it is 359s. The model is primarily designed for use over the period of a storm hydrograph and it is inappropriate in that context to consider calculation time steps m u c h in excess of 0.5 h because of ti~e expected rate of change o~ ~oil moisture potentials and discharges. The results with a time step of < 0.5 h differ very littlefrom those of the 0.5 h steps, and the former use considerably greater processing times than the latter,suggesting that the 0.5 h step is an appropriate and efficientchoice for this type of problem. [It is a subject of further research
177 2500
~,b i 8187. s
tb ~ s~s2
9a • ,1tb
20oo C..p.U. time
(s)
*6b
1500
la,
,8b
lOOO
* t3b ,110 *1Oh
OUU
5b*
?.b za~
0
to
4bo 4~,
®7b
,3a *t3a * lOa
e60
,8o
,JZb 50® ®Ta ®l?.a , I t00 t000 total nu.mber of nodes
I
10000
Fig. 5o Relationship of computer processing time to number of nodes in grid.
to investigatv the use of such a model for baseflow, in which case lon~er time steps could be considered.] In everyday use of the program the subsurface calculation time step can be flexible depending on the number of iterations necessary for solution: an initial and 8 rn,~ximum subsurface calculation time step are specified. With practical cons.~.~lerationsin mind, two measures of discharge predictive success are plotted against processing time for the 0.5 h time step runs in Fig. 6. The discharge errors relate, again because of prac~:i~alconsiderations, to the period after mo&~l run-in effects can be expected to have dispersed. To exclude the run-in effects of all the runs in this comparative work, discharge predictions are considered after 25 h of simulation, that is,from the beginning of the rainstorm. The two measures plotted are the percentage of half-hourly
178 t.0 baseline run
ul (i) =3
erie A
t0a~e :t3a e
[3
Z0
8a
e70
. °0.~ t=
A
:ga
o
U
6a
(..
&
o
o.61_%
g,o
e
3a
® within ~ 10~ of |~!~oseof baseline run
&
A
,,
+ - 5~
.
.
.
.
.
.
.
.
.
.
e
Q,
0.41 ~a 0
I
500
I
I
t000 t500 c.p.u.fime (s)
I
Z000
I
ZS00
Fig. 6. Relationship of processing time to accuracy of discharge predictions.
v~.lues of discharge within _+10% and _+5% of the values of the baseline run. It should perhaps be noted at this stage that in practical use the mode] is likely to be calibrated against discharge. F i b r e 6 suggests that, in terms of cost-effectiveness, simulations 7a, 8a, 10a, 12a and 13a are the best at the + 10% level, offering reasonably accurate predictions without excessive run times. These runs are those with regular meshes of Ax/Az not greater than 20, or of non-regular mesh as used here, and with total node numbers less than about 750 (which averages about one node per square metre of the hillslope section in th~se cases). The runs ~uich perform best at the _+5% prediction level are runs 8a, 10a and 13a. They have Ax/Az values of 10, 20 and variable bat, of the five runs above, have the grep~test numbers of nodes, that is, the generally smaller elvments. CONCLUSIONS
Given the type of rainfall-runoff model methodology and the nature of the subsurface storm~ow response considered here, certain guidvlines car be given Zor goo~l hydrological prediction and cost-effective modelling. These points are as follows. (1) With regard to the overall prediction of hydraulic potential values over the hillslope and the reduction of oscillation in the ~olution, generally smaller elements are favoured with a ratio of horizontal to vertical element dimension of ~ 20 for regular grids, ideally ~<10. (2) With xegard to hillslope discharge prediction over the course of the hydrograph, it is again wise to have generally smaller elements ~nd Ax/Az ~ 20 for regular grids. It is also acceptable to have an i~egular grid (see details above) with fining of the mesh towards the seepage face and the infiltration boundary.
179
(3) In addition to accuracy of prediction, finer grids have the advantage of reducing the necessary run-in time at the beginning of ~imuiation. Lower Ax/Az values for a given total number of nodes reduce run times. (4) With this methodology neither number of iterations required nor cumulative mass balance serve as a distinguishing factor for choice between different calculation meshes. (5) If At = 0.5 h is the maximum time step considered for the solution time step for the response of soil conditions to storm rainfall, this value is seen to serve efficiently, giving little if any loss of accuracy over shorter time steps and using shorter c.p.u, times. (6) There is a conflict between the better discharge predictions of smaller elements and the associated longer c.p.u, times. These experiments suggest that with the meshes recomme~ded above (a regular mesh of Ax/Az ~ 20 or the irregular meshes as used here) a higher number of nodes (here an average of one node per 1-2 m ~ of domain area) should be used if the emphasis is on more accurate prediction and a lower number of nodes (here one node per 2-10 m 2) if shorter execution time is the priority. It is appreciated that in particular applications physical inhomogeneities may modify these mesh considerations. The above guidelines, used within the context of a particular case, should ensure that inaccuracies of numerical approximation as influenced by time and space steps are reduced as far as practically possible. ACKNOWLEDGEMENT
This work has been funded by the U.K. Ministry of Agriculture, Fisheries and Food. REFERI~INCES
Beven, K., Calver, A. and Morris, E.M., 1987. The Institute of Hydrology Distributed Model. ?ep. No. 98, Institute of Hydrology, Wallingford. Calver, A., in press. Calibration, sensitivity and validation of a physically-based rainfall-nmoff model. J. Hydrol., 103: 103-115. Clapp, R.B. and Hornberger, G.M., 1978. Empirical equations for some soil hydraulic pror~rties. Water Re,our. Rcs., 14: 601--604. Coole,v, R.L., 1983. Some new procedure~ for numerical ~olution of variably saturated flow problems. Water Resour. Res., 19: 1271-1285. Huyakorn, P.S. and Pinder, G.F., 1983. Computational Methods in Subsurface Flow. Academic Press, New York, N.Y., 473 pp. Lynch, D.R., 1984. Mass conservation in finite element groundwat~.r models. Adv. Water Resour., 7: 67-75. Richards, L.A., 1931. Capillary conduction of liquids through ~)orous mediums. Physics, 1: 318-333.