Daniel P. Ames, Nigel W. T. Quinn, Andrea E. Rizzoli (Eds.) http://www.iemss.org/society/index.php/iemss-2014-proceedings
Feedback versus uncertainty
Ronald R. P. van Nooijena, Markus Hrachowitzaand Alla G. Kolechkinab aDelft University of Technology, Delft, the Netherlands (r.r.p.vannooyen@tudelft.nl)
bAronwis, Den Hoorn Z.H., the Netherlands
Abstract: Even without uncertainty about the model structure or parameters, the output of a
hydrolog-ical model run still contains several sources of uncertainty. These are: measurement errors affecting the input, the transition from continuous time and space to discrete time and space, which causes loss of information about the input, discretization of the model equations resulting in errors due to the discretization scheme and the use of finite precision calculations in model evaluation. Interval analy-sis can provide upper bounds on the output error due to all of these sources. This paper focuses on tracking uncertainty about input values and the effects of finite precision calculation. A conventional hydrological model was recoded in interval arithmetic and run with small, but non-zero uncertainty bounds on the input time series values to see whether the internal feedback loops would keep the output uncertainty bounded.
Keywords: interval analysis; conceptual model; uncertainty; feedback.
1 INTRODUCTION
Conceptual models of environmental systems are attractive tools in water resources management because their relative simplicity seems a good match for the relative lack of data still common in this field. When using models there are several obstacles to overcome. Some of these are general problems shared by all fields of scientific research. Measurements of quantities in a natural system invariably entail a certain amount of uncertainty, see for example Beven et al. [2011] or Refsgaard et al. [2007], but for this there is a long tradition of error estimation in science. Incomplete knowledge of the variation of a quantity or parameter in time also has a long history. However, dealing with the incomplete knowledge of parameters or variables in space is still a field of active study, see for example Torquato [2002]. In conceptual models large areas are often modeled by a essentially zero dimensional sub-system, which is one way of dealing with incomplete knowledge of parameters and variables in space. Usually these models are continuous time models. When numerical results are desired there are two options: exact or approximate analytical expressions for the solution or a discretization of the continuous time model.
Evidently the analysis of measurement error, the derivation of bounds on variability in time and space and the development of effective ways of dealing with variability in space, all depend strongly on hy-drological and water resource management knowledge. However, the interaction between uncertainty about model input, the model characteristics and the method of approximation of the model solution are of a more mathematical nature. Given the relatively large uncertainties in input that add uncertainty with each time step, there must be a mechanism that causes a decrease in uncertainty or else hydro-logical models would be useless for long term studies. One such mechanism is feedback. The most obvious feedback loop in hydrological models is caused by the proportionality of the outflow to the level in the reservoir. This tends to move the separate reservoirs in the system towards zero storage. In the absence of inflow the level of a linear reservoir would move towards zero as exp (−ct) for some c > 0.
This paper examines the interaction of input uncertainty and feedback using interval analysis, a tech-nique discussed in Alefeld and Herzberger [1983], Birdie and Surana [1998], Corliss [1989], Hayes [2003], Jaulin et al. [2001] and Moore et al. [2009]. Interval analysis was used in earlier paper van Nooijen and Kolechkina [2012] to study the effect of loss of information due to averaging in time [Kirch-ner et al., 2004] for a simple model. In contrast this paper uses interval analysis to study the effect of uncertainty in the measurement value on the state and the output of the discretized model. The effect of input errors on the model output can be studied in different ways as discussed in Beven and West-erberg [2011], McMillan et al. [2012], Beven [2006], Beven and Binley [1992], G ¨otzinger and B ´ardossy [2008] and Hrachowitz et al. [2013a]. Interval analysis has both a theoretical and a practical side. It is performed on a computer, but provides hard bounds on the error for a given range of input errors. Now if it were necessary to rewrite the model code to perform this analysis or if it took inordinate amounts of computer time then this method would not be interesting. The reason we examine the method is that in a modern computer language we may not need to write a completely new program and we may not need different code for production runs in normal arithmetic and exploratory runs in interval arithmetic. Moreover, a run of a program in interval arithmetic may not be that much slower than a run of production code. Finally, even if it were 1000 times slower, one run can provide lower and upper bounds for all possible combinations of errors up to 40% in a 10 year daily precipitation and potential evaporation time series. Especially for non linear models, where fully exploring the space of possible inputs is likely to be important, this can be very useful.
This paper reports of the application of interval analysis to the discretization of a non-linear model, taken from Hrachowitz et al. [2013b]. It discusses problems encountered while writing the software used and reports whether or not the results are unduly pessimistic. This paper uses explicit (forward) Euler as its discretization method. The code contains a check on the appropriateness of the time step size. More information on the role of discretization schemes can be found in Clark and Kavetski [2010] and Kavetski and Clark [2010].
2 NOTATION
In conceptual hydrological models there are usually special conditions and “if ... then ... else ...” constructs. To explicitly include these in the written equations is not always easy. It can be helpful to be able to translate a statement that is either false or true to 0 (for false) and 1 (for true). For instance, to put the rule that evaporation eevap(t)from a reservoir with inflow p (t) and current stored volume s (t) occurs at the potential evaporation rate epot(t), but there is no evaporation from an empty reservoir, so when s (t) = 0, we would need to write
ds (t) dt = f (s, t) (1) where f (s, t) = ( 0 s (t) ≤ 0and p (t) − eevap(t) < 0 p (t) − eevap(t) s (t) > 0or p (t) − eevap(t) ≥ 0 (2) For such cases Kenneth E. Iverson introduced the concept of the Iverson bracket, which Donald E. Knuth extended in Knuth [1992]. It is a notational convenience. Suppose we have a condition C that can be true or false, like x < 1 or 0 ≤ y ≤ 1, in that case we define
JC K = (
0 : C is false
1 : C is true (3)
which allows us to put all information contained in Eq. 1 and Eq. 2 into one equation: ds (t)
dt =Js (t) > 0KJp (t) − eevap(t) ≥ 0K (p (t) − eevap(t)) (4) Please note that while the Iverson bracket is well established, the use of this type of brackets,J. . .K , is not a universal standard, there are several competing notations. We will also use a variation on strong zeros as discussed in Knuth [1992], when we write
we will interpret this as
JC (x)Kf (x) = (
0 C (x) is false
f (x) C (x) is true (6)
which allows us to ignore any problems that f might have when C (x) is false. The notation allows us to force a function to be zero even where it is not defined (hence the term strong zero). This notation also highlights points where we may expect problems. Two other nice abbreviations that will be used are f+for the positive part of a function
f+(t) = f (t)Jf (t) ≥ 0K = max (0, f (t)) (7) and f−for the absolute value of the negative part of a function
f−(t) = −f (t)Jf (t) ≤ 0K = − min (0, f (t)) (8) which allows us to write Eq. 4 as
ds (t) dt = (p (t) − eevap(t)) + −Js (t) > 0K (p (t) − eevap(t)) − (9) In the remainder of the paper we will often omit the function argument of functions that are a function of time only and use ds/dt in stead of ds (t) /dt.
3 INTERVAL ANALYSIS
The basic principle of interval analysis is that we replace all numbers by closed intervals. All function and operations on numbers are likewise replaced by functions and operations that work on intervals. When executing this replacement the following conditions must be satisfied.
For functions: if we replace a function f on real numbers by a function ˜f then
{f (a) : a ∈ A} ⊆ ˜f (A) (10)
must hold.
For operations: if we replace addition (+) by then
{a + b : a ∈ A, b ∈ B} ⊆ A B (11)
and if we replace multiplication (·) by then
{a · b : a ∈ A, b ∈ B} ⊆ A B (12)
and so on.
3.1 A basic problem with standard interval analysis
Suppose we have a simple linear reservoir with linear outflow, ds dt = p − cs (13) s (t0) = s0 (14) If we define ¯ pk= ´tk+1 τ =tkp (τ ) dτ tk+1− tk (15)
then explicit Euler approximation would be
s (tk+1) = s (tk) + (tk+1− tk) pk− (tk+1− tk) cs (tk) (16) or in a mathematically equivalent form
s (tk+1) = (1 − (tk+1− tk) c) s (tk) + (tk+1− tk) pk (17) In standard interval analysis the following problem occurs for explicit Euler. If we program Eq. 16, then for normal interval arithmetic the computer has no way of knowing that the first and the second occurrence of s (tk)refer to the same variable. Once s (tk)has a non-zero width this will result in quickly expanding intervals for s (tk)in the iteration. If however we use the Eq. 17 then the same condition that is needed for the stability of the explicit Euler iteration, that is |1 − (tk+1− tk) c| < 1, guarantees that our intervals will not expand. Evidently this difference in implementation conflicts with the natural desire to program the algorithm without having to think about whether or not we are using interval arithmetic. For this simple case a variation on interval arithmetic called affine interval arithmetic may solve the problem. The question is “what to do with more complex right hand sides”, further study is needed to answer this question. For this paper the code was adapted by hand to avoid the cancellation problem of interval arithmetic, that is the fact that the general rule for subtraction, [a, b] − [c, d] = [a − d, b − c], results in [a, b] − [a, b] = [a − b, b − a] even when the two intervals represent the same number.
4 ACONCEPTUAL MODEL
In Hrachowitz et al. [2013b] three conceptual models are given, we will use the “Allt Coire nan Con” model without the tracer part. The relevant equations are reproduced here.
cR Runoff generation coefficient [-] qSF Runoff from fast reservoir [L T-1]
eP Potential evaporation/transpiration [L T-1] qOF Overland flow [L T-1]
kF Storage coefficient of fast reservoir [ T-1] qSS Runoff from slow reservoir [L T-1]
kS Storage coefficient of slow reservoir [ T-1] sF Storage in fast reservoir [L]
lP Transpiration threshold [-] sF,max Storage capacity of fast reservoir [L]
pmax Percolation capacity [L T-1] sS Storage in slow reservoir [L]
p Rainfall [L T-1] s
U Storage in unsaturated reservoir [L]
β Shape parameter sU,max Capacity of unsaturated reservoir [L]
4.1 Equations for the model
The selected model uses only the unsaturated zone, the slow and the fast reservoir. It does not use preferential recharge. Moreover, for this model the calibration showed that the convolution is not needed at a daily time scale. If the notation introduced earlier is used then we can write the equations as follows dsU(t) dt = f + u s + U, min s + F, sF,max , t −JsU> 0Kf − u sU, min s + F, sF,max , t (18) dsF(t) dt =JsF < sF,maxKf + f s + U, s + F, t −JsF> 0Kf − f s + U, min (sF, sF,max) , t (19) dsS(t) dt = f + s s + U, s + S, t −JsU> 0Kf − s sU, s + S, t (20) where fu(sU, sF, t) = (1 − cR(sU)) pE−JsU> 0Kpmax sU sU,max − min 1, sU sU,maxlP JsU> 0K s U sU+ s + F eP (21) ff(sU, sF, t) = cR(sU) pE−JsF > 0K sF s+U + sF eP−JsF> 0KkFsF (22) fs(sU, sS, t) =JsU> 0K pmax sU,max sU−JsS> 0KkSsS (23)
with cR(y) = 1 1 + exp1β12− y sU,max (24)
and the drainage flows are given by
qSF= kFs + F (25) qOF=JsF≥ sF,maxKf + f s + U, sF,max, t (26) qSS= kSs + S (27) 4.2 Feedback
In this model the feedback terms are the term containing pmaxin Eq. 21, the term containing kF in Eq.
22 and the term containing kSin Eq. 23. These terms will tend to move the levels in the continuous
model back to zero. As long as our discretization is well behaved the same will hold for the discretized equations.
It is not surprising that hydrological models can deal with small errors in input and initial conditions, after all these are never known exactly, so any model that is intended to extend system behavior into the future should still be able to do so after a small perturbation. What is somewhat unexpected is that even the non-linear models can deal with quite large errors. It is also not surprising that the discretized models approximate the continuous solution, what is surprising is that they seem to do so even when the time steps are relatively large. The fact that the good behavior is not entirely obvious suggests that it would be nice to have methods to detect potential problems with time step size or non-linearity.
5 NUMERICAL EXPERIMENTS
The model parameters from Hrachowitz et al. [2013b] were used. The aim of this paper was to exam-ine whether or not the feedback in the model would keep propagation of known errors under control. This should either work or not work independent of any relation between input data and calibration data. So, as the data used in Hrachowitz et al. [2013b] is not yet publicly available, other precipitation and evaporation series were used to generate input. Hourly precipitation and daily Makkink evapora-tion, made available by the “Koninklijk Nederlands Meteorologisch Instituut” KNMI (Royal Netherlands Meteorological Institute), were used to create input files. The hourly precipitation was first converted to daily precipitation.
5.1 Interval experiments
We perform an interval experiment with inputs [0.6, 1.4] · p and [0.6, 1.4] · eP(so we imposed 40% error
margins). Runs were done for the period from 1982 to 1992. Precipitation (Fig. 1), evaporation (Fig. 2), bounds on unsaturated zone storage (Fig. 3) and bounds on the sum of slow, fast and overland flow (total outflow) (Fig. 4) for 1987 are shown here.
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 0 5 10 15 20 25 30 35 40 daily precipitation [mm] Figure 1. Precipitation 1987
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
0 1 2 3 4 5 6
daily potential evapotranspiration [mm]
Figure 2. Evaporation 1987
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan
80 100 120 140 160 180 200 220
unsaturated zone storage [mm]
Figure 3. Unsaturated zone storage 1987
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan
-10 0 10 20 30 40 50 60 70 80 total outflow
Figure 4. Total outflow 1987
The model under study does not seem to be very sensitive to input error for the given range of input data. The full 10 year output is not shown, but it has the same yearly pattern (more uncertainty in summer than in winter) for the state as found for 1987 (Fig. 3). A run where the evaporation was given 0% error showed a more even distribution of uncertainty over the year. Recall that all runs were done with fixed parameters not related to the input, no calibration on the input data was done as the goal was to study how the model dealt with input uncertainty for known model parameters and structure.
5.2 Check on degree of pessimism of interval arithmetic
To check the results obtained using interval arithmetic the input data was modified as follows. Run “a” was done with 1.4 · p and 0.6 · eP and run “b” with 0.6 · p and 1.4 · eP. For the interval run the
following inputs were used: [0.6, 1.4] · p and [0.6, 1.4] · eP. We examined the differences between the
lower bounds of the intervals and the time step wise minimum of runs “a” and “b” and of the upper bounds of the intervals and and the time step wise maximum of runs “a” and “b” and found them to be small (less than 0.1 percent of the level) for the unsaturated zone and the slow reservoir, so the interval method is not unduly pessimistic there. However, for the fast reservoir there were some spikes in the differences because in the version of the code used for these runs there is no attempt to avoid cancellation problems for the fast reservoir. Moreover, the fast reservoir has a “all or nothing” threshold for overland flow that complicates matters for interval arithmetic.
6 DISCUSSION
From the results it is clear that the feedback provided by the relation between reservoir level and reservoir outflow does indeed limit the uncertainty in the output error. Comparison of runs with and without evaporation shows that uncertainty in evaporation is a major contributor to over all uncertainty. This is clearly seen in the multi-year simulation where the additional uncertainty built up during the relatively wet winters is removed each year.
The interval method allows calculation of error bounds on the computer model output from error bounds on the model input in one model run, both for linear and highly non-linear models. This is a quick way to examine whether or not the model is sensitive to errors in the input and if so, for which regimes. The same method can also be used to study sensitivity to errors in the model parameters for given parameter ranges. The feedback of the continuous system and the stability and consistency of the numerical scheme keep the error propagation under control. Whether or not the method can be applied without some degree of rewriting the model code is a matter for future investigation.
REFERENCES
Alefeld, G. and Herzberger, J. (1983). Introduction to Interval Computation. Academic Press. Beven, K. (2006). A manifesto for the equifinality thesis. Journal of Hydrology, 320:18–36.
Beven, K. and Binley, A. (1992). The future of distributed models: model calibration and uncertainty prediction. Hydrological processes, 6(3):279–298.
Beven, K., Smith, P. J., and Wood, A. (2011). On the colour and spin of epistemic error (and what we might do about it). Hydrol. Earth Syst. Sci., 15(10):3123–3133.
Beven, K. and Westerberg, I. (2011). On red herrings and real herrings: disinformation and information in hydrological inference. Hydrological Processes, 25(10):1676–1680.
Birdie, T. R. and Surana, K. S. (1998). The use of interval analysis in hydrologic systems. Reliable Computing, 4:269–281.
Clark, M. and Kavetski, D. (2010). Ancient numerical daemons of conceptual hydrological modeling: 1. fidelity and efficiency of time stepping schemes. Water Resources Research, 46:W10510 pp. 1–23. Corliss, G. F. (1989). Survey of interval algorithms for ordinary differential equations. Applied
Mathe-matics and Computation, 31:112 – 120.
G ¨otzinger, J. and B ´ardossy, A. (2008). Generic error model for calibration and uncertainty estimation of hydrological models. Water Resources Research, 44:W00B07 pp.1–18.
Hayes, B. (2003). A lucid interval. American Scientist, 91(6):484–488.
Hrachowitz, M., Savenije, H., Bl ¨oschl, G., McDonnell, J., Sivapalan, M., Pomeroy, J., Arheimer, B., Blume, T., Clark, M., Ehret, U., Fenicia, F., Freer, J., Gelfan, A., Gupta, H., Hughes, D., Hut, R., Montanari, A., Pande, S., Tetzlaff, D., Troch, P., Uhlenbrook, S., Wagener, T., Winsemius, H., Woods, R., Zehe, E., and Cudennec, C. (2013a). A decade of predictions in ungauged basins (PUB)–a review. Hydrological Sciences Journal, 58(6):1198–1255.
Hrachowitz, M., Savenije, H., Bogaard, T. A., Tetzlaff, D., and Soulsby, C. (2013b). What can flux tracking teach us about water age distribution patterns and their temporal dynamics? Hydrology and Earth System Sciences, 17(2):533–564.
Jaulin, L., Kieffer, M., Didrit, O., and Walter, ´E. (2001). Applied interval analysis. Springer-Verlag London Ltd., London.
Kavetski, D. and Clark, M. P. (2010). Ancient numerical daemons of conceptual hydrological modeling: 2. impact of time stepping schemes on model analysis and prediction. Water Resources Research, 46:W10511 pp. 1–27.
Kirchner, J. W., Feng, X., Neal, C., and Robson, A. J. (2004). The fine structure of water-quality dynamics: the (high-frequency) wave of the future. Hydrological Processes, 18(7):1353–1359. Knuth, D. E. (1992). Two notes on notation. American Mathematical Monthly, 99(5):403–422.
McMillan, H., Krueger, T., and Freer, J. (2012). Benchmarking observational uncertainties for hydrol-ogy: rainfall, river discharge and water quality. Hydrological Processes, 26(26):4078–4111.
Moore, R. E., Kearfott, R. B., and Cloud, M. J. (2009). Introduction to interval analysis. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA.
Refsgaard, J. C., van der Sluijs, J. P., Højberg, A. L., and Vanrolleghem, P. A. (2007). Uncertainty in the environmental modelling process – a framework and guidance. Environmental Modelling & Software, 22(11):1543–1556.
Torquato, S. (2002). Random heterogeneous materials: Microstructure and macroscopic properties. Springer-Verlag, New York.
van Nooijen, R. and Kolechkina, A. (2012). A problem in hydrological model calibration in the case of averaged flux input and flux output. Environmental Modelling & Software, 37(0):167 – 178.
