- Email: [email protected]
A combined power-law and exponential model for streamflow recessions
Journal of Hydrology (2008) 352, 157– 167
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jhydrol
A combined power-law and exponential model for streamflow recessions Qiuming Cheng
*
Department of Earth and Space Science and Engineering, Department of Geography, York University, 4700 Keele Street, Toronto, Canada M3J 1P3 State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan and Beijing, China Received 6 November 2007; accepted 9 January 2008
KEYWORDS Non-linear process; Fractals; River flow; Power-law; Exponential model; Combined model
Summary This paper proposes a non-linear model to describe streamflow recessions. A model of combining power-law and exponential functions is derived from a non-linear differential equation associating outflow change rate and flow itself with a time-dependent decreasing decay rate. The decay rate is approximated by the first three terms in a Taylor series expansion. The solution of the resulting differential equation includes three components of the expansion determining two exponential functions and a power-law function, respectively. It is demonstrated analytically and empirically that the introduction of a non-linear differential equation comprising a mixture of power-law and exponential forms is superior to using an ordinary differential equation with exponential solution, or a nonlinear equation with power-law solution, for modeling the decay patterns of flow in river systems. Analytical properties of the mixing model are discussed and the model is subsequently validated by characterizing the relationships between decay of peak flow as a function of time after the initial occurrence of flow peak in 10 datasets of flow events chosen from five gauging stations in the Oak Ridges Moraine (ORM) area, Canada. ª 2008 Elsevier B.V. All rights reserved.
Introduction Quantitative description of the flow recession process in a river system and its relationship with precipitation is one * Address: Department of Earth and Space Science and Engineering, Department of Geography, York University, 4700 Keele Street, Toronto, Canada M3J 1P3. E-mail address: [email protected]
of the most important tasks of contemporary hydrology. As an important surface water system, river flow has significant influence on many other environmental systems. Understanding river flow recession processes is not only scientifically important but also essential for water resource management and geo-hazard (flooding, landslide, etc.) prevention. Therefore, river gauging and flow monitoring have being conducted worldwide for many river systems. These programs have generated a large amount of data which
0022-1694/$ - see front matter ª 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2008.01.017
158
Q. Cheng
can be used for empirical or analytical modeling of the properties of river systems. Detailed study of river flow recession processes has been increasingly attracting attention from researchers in academic institutions, industries and governments at various levels. Prediction of runoff in river systems is also important for the insurance industry, transportation and other relevant sectors. To develop a practical system for prediction of river flow requires several new developments including sophisticated models, calibration–validation techniques, model generalization and an efficient on-line system to disseminate the results publicly. A number of mathematical models have been proposed in the literature for quantifying the recession process of river flow. Some of these models, for example, rainfall–runoff models (RRMs), have become standard tools that are routinely used for hydrological modeling in engineering and environmental science (Littlewood et al., 1997; Sefton and Howarth, 1998; Wagener et al., 2004; Littlewood, 2002; Kokkonen et al., 2003; Yuan and Cheng, 2005). However, an optimum model for prediction purposes in all types of catchment environments at various scales does not yet exist. While most authors use a constant-coefficient exponential model derived from a linear differential equation, nonlinear models have also been investigated and applied. It is generally accepted that linear models are suitable for simple catchment environments but non-linear models might be the most appropriate tool for more complex catchment systems. There are two kinds of recession models in which non-linearity has been introduced: the best known form of non-linearity is related to the basic non-linear differential equation governing unsteady flow from a large unconfined aquifer in a stream channel presented by Boussinesq in 1877 (Hall, 1968; Tallaksen, 1995; Brutsaert and Nieber, 1977). In this context, recession rate of flow Q(t) can be associated with storage S(t) as SðtÞ ¼ sQ n ðtÞ
ð1Þ
where s is the storage coefficient and n is the non-linearity parameter. Main cause of the non-linearity of the relationship between rainfall and runoff is the antecedent moisture condition (Beven, 2000; Tallaksen, 1995; Wittenberg, 1999). If n = 1, the model becomes the most commonly used linear form, which has the advantage of computational efficiency. Another mathematical description of flow systems is that changes in storage S(t) of a hydrological system in time depend on inflow U(t) and outflow Q(t) that are related by dSðtÞ ¼ UðtÞ Q ðtÞ dt
ð2Þ
In dealing with river flow mainly caused by rainfall the U(t) can be considered as effective rainfall, being the contribution of rainfall to the river flow. If there is no inflow, U(t) = 0, or if the amount of U(t) is so small that it can be neglected, then Eqs. (1) and (2) can be combined as dQ ðtÞ 1 ¼ Q 2n ðtÞ dt sn
ð3Þ
This approach results in different types of non-linearity depending on the value of n: (1) When n 5 1, the non-linear differential equation can be solved and its solution is
1
Q ðtÞ ¼ Q ðt0 Þ½1 þ cðt t0 Þn1
ð4Þ
Q 1n ðt0 Þ is a constant and Q(t0) is the value of where c ¼ n1 ns Q(t) at initial time t0. This result is similar to that given by Brutsaert and Nieber (1977), Coutagne (1948) and Moore (1997). The relation between log Q(t) and log [1 + c(t t0)] can be plotted as a straight-line on log–log paper. When n = 1/2, the model (4) generates a hyperbola form as has been discussed by many authors (Boussinesq, 1904; Hall, 1968; Brutsaert and Nieber, 1977). (2)
When n = 1, Eq. (3) has the exponential solution 1 Q ðtÞ ¼ Q ðt0 Þ exp ðt t0 Þ ð5Þ s where 1/s is a time independent parameter determining the decay rate of runoff Q(t) without effective rainfall (Maillet, 1905; Mitchell, 1972). Q(t0) represents initial river flow at time t0, and Q(t) is river flow at time t. Since the parameter s is the time it takes for the value Q(t0) to recess to exp (1) of its original value, Q(t0) can be considered as 1 for unit hydrograph. The function (5) can be plotted on semi-logarithmic paper as straight-line with slope equal to 1/s. This is the most commonly used model for describing hydrograph recessions and river peak flow recession. However, a number of studies have shown that the function (5) is often not sufficient to represent the relationship because the plot shows curvature and this has led to the investigation of other, non-linear, models. The other type of generalization of non-linearity is based on modification of the coefficient s in model (3). Croke (2006) has recently pointed out that the assumption that the unit hydrograph has a constant-coefficient exponential form is not valid from a conservation of mass point of view. He suggested a power-law with monotonically decreasing coefficient as an alternative form for the unit hydrograph. Cheng et al. (2006) noticed this suggesting the use of a differential equation with time-decreasing coefficient yielding a power-law form for describing decay patterns of flow. From these points of view the constantcoefficient exponential model only represents a special solution for catchment channel modeling with limited applicability. In Croke (2006), a power-law model was proposed to describe river flow without effective rainfall and it was implemented using a sum of multiple exponential terms for approximation. That author also pointed out that the exponential terms used for representing the power-law do not necessarily have any physical significance. A similar model involving a sum of multiple exponential terms had already been used by others including Werner and Sundquist (1951). Other non-linearity treatments include Napiorkowski and Strupczewski’s (1979) cascade of non-linear reservoirs described by a Volterra series. Ishihara and Takagi (1965) modeled recession flow as a sum of two components: outflow from a confined aquifer subject to a simple exponential decay (5), and outflow from an unconfined aquifer expressed by means of a hyperbola form (4). It can be seen from the above review that although exponential and power-law models are two most commonly used models, application of either of these is often not sufficient for describing the behavior of river flow recession. Some efforts have been made to investigate
A combined power-law and exponential model for streamflow recessions combined multiple parameter models, for example using a sum of multiple exponential terms, but these types of combinations of the same types of functions are neither physically based nor straightforward for implementation. The current paper attempts to introduce a new mixing model combining exponential and power-law forms derived from a non-linear differential equation with a variable coefficient function that is approximated by the first three terms of a Taylor series expansion. This new combined model is demonstrated to be appropriate for characterizing streamflow recessions in various types of catchment. Comparisons between the power-law, exponential and combined models were conducted both analytically and empirically using river flow datasets for five gauging stations selected from the Greater Toronto Area, Canada.
Non-linear model for describing river peak flow recession In order to discuss the two suggested non-linear extensions to the flow model, we will start with Eq. (3) with the coefficient s as a function of time t: dQ ðtÞ 1 ¼ Q ðtÞ dt sðtÞ
ð6Þ
In order to solve Eq. (6), one needs to know the functional form of s(t). It is usually not possible to know the exact expression of this function since it may have different forms depending upon the dominant system of river flow and the characteristics of the catchments channels. For example, in some simple cases if the change of flow rate dQdtðtÞ is mainly due to magnitude of Q(t) or if the flow is dominated by direct surface flow in small catchments with simple regular channel geometries, the coefficient should be independent of the time duration after the peak flow. In other catchments, however, decreasing recession rate may not only be proportional to the flow itself but may also decrease through time due to channel geometry changes; for example, when flow in the channel reduces, wetted channel cross-sectional area and saturation of the catchment basin reduce as well. In this case, the coefficient might be expressed as 1/[s(t)] = 1/[st] which reduces monotonically as t increases. The solution of the flow differential equation assumes the power-law form Q ðtÞ ¼ Q ðt0 Þ
1=s t t0
ð7Þ
More complex situations might be applicable to actual river channel systems because wetted channel cross-sectional area and saturation of the catchment basin reduce as Q(t) decreases through time t. This situation can be characterized by 1/[s(t)] = 1/[stm] where m is a positive constant number (m > 1) which yields the solution 1 ðt1m t1m Q ðtÞ ¼ Q ðt0 Þ exp Þ ð8Þ 0 sð1 mÞ In order to understand the properties of this solution for different non-linearity (m), we will determine the ratio of flow for two different time periods.
159
(1) The case of m = 0. When m = 0, the solution satisfies the simple exponential form (5) with time interval t t0. At time interval t t0 = s, Q(t) recesses to Q(t0)exp (1) with the following properties: (a) When t ! 1, then Q(t) ! 0. (b) The ratio of Q(t) at two different points of time, say t1 and t2, is Q ðt1 Þ 1 ¼ exp ðt1 t2 Þ ð9Þ Q ðt2 Þ s This property indicates that the ratio of Q(t1) and Q(t2) depends only on the time interval (t1 t2) and is independent of t1 and t2 (t2 < t1). For a given time interval, for example, t1 t2 = k where k is a constant, the ratio of Q(t1) and Q(t2) will remain constant no matter how much time has expired after the initial time t0 (e.g. if t1 and t2 t0). The parameter 1/s = [log Q(t1) log Q(t2)]/[t1 t2], represents the ratio of log difference between Q(t1) and Q(t2) and the difference between t1 and t2. (2) The case of m = 1. When m = 1 in (6), the decay rate of Q(t) depends not only on Q(t) itself but also on 1/[st] implying the decay rate decreases with time t. The decay of Q(t) follows a power-law relationship with time t. For time ratio t/ t0 = ks, Q(t) is reduced by the factor k to Q(t0)/k with the following properties: (a) When t/t0 ! 1, then Q(t) ! 0. (b) The ratio of Q(t) at two different points of time, say t1 and t2, is 1s Q ðt1 Þ t1 ¼ ð10Þ Q ðt2 Þ t2 This property indicates that the ratio Q(t1)/Q(t2) is proportional to the ratio (t1/t2)1/s, and the parameter 1/s = [log Q(t1) log Q(t2)]/[log t1 log t2], represents the ratio of log difference between Q(t1) and Q(t2) and the log difference between t1 and t2. For comparison purposes we can check the ratio of Q(t2 + k)/Q(t2) at the two different times (t1 and t2) with constant time interval t1 t2 = k. This illustrates that, for any given constant k, when t1 = t2+k ! 1, Q(t1) tends to become the same as Q(t2). (c) According to (10), if t1/t2 = 1/k (constant), then Q(t1)/ Q(t2) = k1/s. For example, if k = 2, then Q(t)/Q(2t) = 21/s, implying that in order to maintain the same reduced rate from Q(t) to Q(2t) one needs to keep constant ratio of t1/t2 = 1/2. Then t2 = 2 days if t1 = 1 day after the initial time t0; and t2 = 10 days if t1 = 5 days after the initial time t0. This indicates that the decay becomes slower moving away from the initial time. (3) The case of m > 1. When m > 1, the solution (8) shows that the decay of Q(t) follows an exponential relationship with time interval t1m t1m . For time interval t1m t1m ¼ sð1 mÞ, Q(t) 0 0 is reduced to simpler forms with the following properties:
160
Q. Cheng
h i t1m 0 (a) When t ! 1, then Q ðtÞ ! Q ðt0 Þ exp sðm1Þ . (b) The ratio of Q(t) at two different points of time, say t1 and t2, is 1m Q ðt1 Þ 1 t1 t1m ð11Þ ¼ exp 2 Q ðt2 Þ sð1 mÞ This property indicates that the ratio of the Q(t) at different times is proportional to the difference t1m t1m , and 1 2 the parameter 1=s ¼ ð1 mÞ½log Q ðt1 Þ log Q ðt2 Þ=½t1m 1 t1m , represents the ratio of log difference between 2 Q(t1) and Q(t2) and (1 m) times the difference between t1m and t1m . The form (11) is the same as the non-linear 1 2 relationship proposed by Horton (1933). It has been shown that the model (6) involves a parameter m determining the non-linearity of the system. The ordinary system with constant-coefficient s becomes the special case of m = 0. Given that the different values of m indicate different decay rates which might be dominated by different types of flow and drainage basins as proposed by many authors, we can explore mixed flow systems by combining the functions with m = 0, 1 and 2 in the preceding expression. Taylor series expansion can be applied to 1/[s(1/t*)] where t*=1/t. In this paper the coefficient function 1/s(t) is approximated by means of Taylor series expansion of 1/[s(1/t*)] resulting in 2 m 1 1 1 1 1 1 1 þ c2 þ þ o ¼ c0 þ c1 sðtÞ t t0 t t0 t t0 ð12Þ where t0 must be positive and t usually is greater than 1 (t0 > 0 and t > 1) so that 1/t 1/t0 is usually negative and much smaller than 1 which guarantees convergence of the summation. The term o( ) is the remainder term and states that there exists a small error for the approximation. If we take the first three terms for approximation, then the form (12) can be rearranged into 1 1 1 1 ¼ a0 þ a1 þ a2 2 þ o 2 ð13Þ sðtÞ t t t where the components a0, a1, and a2 correspond to the coefficients of the first three orders of 1/t. The remainder term states that the error for the approximation approaches zero when t increases infinitely and the order of the error is smaller than the order of (1/t2). Using this approximation one can express the differential model as dQ ðtÞ 1 1 ¼ a0 þ a1 þ a2 2 Q ðtÞ ð14Þ dt t t It can be seen that the Eq. (14) is a combination of Eq. (6) with m = 0, 1 and 2, respectively. Therefore, the model (14) represents the combination of three different types of forms: normal exponential, power-law and exponential of inverse time i a1 h a0 ðtt0 Þa2 1t t1 t 1 0 Q ðtÞ ¼ Q ðt0 Þ ¼ cta1 ea0 tþa2 t e t0 ð15Þ The model (15) gives a new four-parameter model with the parameter c as constant and a0, a1, and a2 representing the coefficients of the first three increment terms,
respectively. In the model (15) both the first term, exp (a0t), and the third term, exp (a2/t), are monotonically decreasing functions of time t representing the decay regulations of flow Q(t) as time increases. An advantage of model (15) in comparison with other multiple parameter non-linear models is that it can be converted into a simple form by applying logarithmic transformation to both sides to give log Q ðtÞ ¼ c a0 t a1 logðtÞ þ a2 t1
ð16Þ
Using this expression, it is straightforward to validate the model graphically and to estimate the parameters using multiple regression. The parameters involved in the combined model can be easily estimated by means of regression applied to log Q(t) against t, log t and 1/t. The multiple correlation coefficient (R2) then can be used as a measure of overall likelihood fitness. This model will be implemented and validated using daily flow data from five gauging stations from the Oak Ridge Moraine area, Canada.
Validation of the new combined power-law and exponential model for characterizing flow recessions The study area on the Oak Ridges Moraine (ORM) is located in southern Ontario, Canada. This is one of the most developed areas in Canada. The Oak Ridges Moraine is one of the most significant landforms in southern Ontario. The moraine is characterized by its rolling hills and river valleys extending 160 km from the Niagara Escarpment to Rice Lake and was formed 12,000 years ago by advancing and retreating glaciers. The moraine contains the headwaters of 65 river systems and has a wide diversity of streams, woodlands, wetlands, kettle lakes, kettle bogs and significant flora and fauna. It is one of the last remaining continuous green corridors in southern Ontario containing 30% forested area. The moraine’s sands and gravel deposits absorb rain and snow melt. The underground water is then stored in layers of sand and gravel (aquifers), filtered and slowly released as cool fresh water to the 65 rivers and streams flowing north into Lakes Simcoe and Scugog and south into Lake Ontario. The greatest threat to the function of the moraine is inappropriate land-use at the surface of the moraine, particularly in headwater areas. Situated on the north shore of Lake Ontario, the Oak Ridges Moraine area experiences a wide range of weather conditions through an average year. Residents of the area have come to expect bitter cold winter and oppressive summer heat. Precipitation includes deep winter snow, heavy spring and summer thunderstorms, and, sometimes, lengthy summer dry spells. The temperature is moderate most of the time and, annually, there is approximately 600– 1000 mm of precipitation including snow in winter season and rainfall in other seasons across the ORM area. Periods of wet or dry weather are not uncommon, often making a year significantly different from normal. Snowfall mainly in December to March accounts for about 15% of the total annual precipitation. In the summer, widespread daily-long rainfalls are rare. Instead, showery precipitation is the gen-
A combined power-law and exponential model for streamflow recessions eral rule. This frequently results in precipitation varying widely from place to place. In addition, this kind of precipitation may come as short but intense downpours that quickly run into drains and streams. A thunderstorm that results in 20 or 30 mm of rain in just one hour does little to relieve the moisture deficit that comes with a long period of dry weather. On occasion, this area receives the effects of the remnants of a tropical storm or hurricane moving through the eastern part of the United States. Occurring from midsummer to autumn, these systems may result in an allday rainfall that can be quite heavy. Precipitation amounts from these systems can also be somewhat variable over relatively short distances. When Hurricane Hazel hit the area in October 1954, 181.6 mm of rain was recorded at Snelgrove station: however, only 38.1 mm at Alton station just 22 km northwest of Snelgrove at the same time. Land-use in the ORM area is largely rural, with forest and agricultural practices dominating the landscape although it includes the Greater Toronto Area that is a highly industrialized area. Increasing population and residential development in this area has led to residential development expanding into the rural areas. Land-use changes, primarily the building of residential subdivisions, the construction of roads and the paving of parking lots,
02HB001
161
increase the imperviousness of the ground surface. Consequently, this surface runoff results in dramatic increases in wet weather flows of the headwater streams on the Moraine causing erosion and degradation of these fragile systems. Various modeling exercises conducted on the study area have demonstrated that runoff volumes observed at the river gauging stations in the drainage basin networks are highly correlated with the baseflow caused by ground water discharge, snow melting, direct flow caused by rainfall as well as the physical properties of drainage basins (Cheng et al., 2006). It has been found that the delay response of surface runoff (river peak flow) to heavy precipitation events varies across the study area depending on the physical properties of drainage networks and stream systems such as the complexity of drainage basin boundaries, glacier sediment types, land-use types, and surface inclination (Ko and Cheng, 2004). It also has been demonstrated that long-term persistency of river flow is related to the characteristics of drainage basins (Ko and Cheng, 2004). The interaction between groundwater and surface runoff also has been investigated so that the influence of groundwater discharge on river flow could be modeled (Cheng et al., 2001, 2006; Cheng, 2004). The conventional SCS method (Soil Conservation Service) of the US Department of Agriculture, Soil Conservation Ser-
DCB
HRB 02HC003
02HC009
CRB
DCB – Duffins Creek Basin 02HC003
02HB002
HRB - Humber River Basin CRB - Credit River Basin
Figure 1 Shaded relief of the digital elevation model of the ORM (DEM data from Kenny, 1997). Lines represent the major basins of the area. White dots represent the locations of five river gauging stations chosen for the study. Labels in the boxes are the IDs of gauging stations.
162 Table 1 Date
Q. Cheng Peak river flow and precipitation observed from five gauging stations in the ORM Precipitation P (mm)
River flow (Q) (m3/s)
t (Day)
Log (t)
Log (Q)
0 1 2 3 4 5 6 7 8 9
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451 0.9031 0.9542
1.3909 0.6590 0.5977 0.4564 0.4564 0.4298 0.3560 0.3096 0.2455
0 1 2 3 4 5
0.0000 0.3010 0.4771 0.6021 0.6990
0.9420 0.6314 0.4564 0.3598 0.2648
Pickering Station on Duffins Creek 02HC006 (201 km2) Dataset 3 for figure (E) and (F) 1954-10-15 78.7 35.7 1954-10-16 0.8 51.8 1954-10-17 1.0 11.5 1954-10-18 0.0 6.17 1954-10-19 0.0 4.84 1954-10-20 0.0 4.28 1954-10-21 0.0 3.96
0 1 2 3 4 5 6
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782
1.7143 1.0607 0.7903 0.6848 0.6314 0.5977
Dataset 4 for figure (G) and (H) 1960-5-12 10.7 1960-5-13 0.0 1960-5-14 0.0 1960-5-15 0.0 1960-5-16 0.0 1960-5-17 0.0
11.4 14.4 7.28 5.64 4.16 3.37
0 1 2 3 4 5
0.0000 0.3010 0.4771 0.6021 0.6990
1.1584 0.8621 0.7513 0.6191 0.5276
Weston Station on East Humber River 02HC003 (820 km2) DataSet 5 for figure (I) and (J) 1954-10-15 125.3 3.79 1954-10-16 1.3 838 1954-10-17 1.3 498 1954-10-18 0 155 1954-10-19 0 87.8 1954-10-20 0 63.7 1954-10-21 0 44.2
0 1 2 3 4 5 6
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782
2.9232 2.4972 2.1903 1.9435 1.8041 1.6454
Erindale Station on Credit River 02HB002 (806 km2) Dataset 6 for figure (M) and (N) 1954-10-15 129.6 1954-10-16 2.3 1954-10-17 1.4 1954-10-18 0 1954-10-19 0 1954-10-20 0 1954-10-21 0 1954-10-22 0
0 1 2 3 4 5 6 7
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451
2.5011 1.8021 1.5092 1.2923 1.0792 1.0000 0.9117
2
Pine Grove Station on East Humber River 02HC009 (260 km ) Dataset 1 for figure (A) and (B) 1956-5-11 23.1 1.95 1956-5-12 3.3 24.6 1956-5-13 3.7 4.56 1956-5-14 0 3.96 1956-5-15 6.8 2.86 1956-5-16 0.9 2.86 1956-5-17 3.1 2.69 1956-5-18 1.9 2.27 1956-5-19 0 2.04 1956-5-20 0 1.76 Dataset 2 for figure (C) and (D) 1960-5-23 11.9 1960-5-24 1.3 1960-5-25 2.5 1960-5-26 0.0 1960-5-27 0.0 1960-5-28 0.0 1960-5-29 0.0
3.23 8.75 4.28 2.86 2.29 1.84 3.48
7.08 317 63.4 32.3 19.6 12 10 8.16
A combined power-law and exponential model for streamflow recessions
163
Table 1 (continued) Date
River flow (Q) (m3/s)
Precipitation P (mm)
t (Day)
Log (t)
Log (Q)
0 1 2 3 4 5 6 7
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451
1.3674 0.9009 0.7284 0.6160 0.6010 0.5378 0.4771
2
Cataract Station on Credit River 02HB001 (174 km ) Dataset 7 for figure (O) and (P) 1954-10-15 69.6 10.2 1954-10-16 1.4 22.2 1954-10-17 3 23.3 1954-10-18 0 7.96 1954-10-19 0 5.35 1954-10-20 0 4.13 1954-10-21 0 3.99 1954-10-22 0 3.45 1954-10-23 0 3 Dataset 8 for figure (Q) and (R) 1974-5-16 29.7 1974-5-17 0.0 1974-5-18 0.0 1974-5-19 0.0 1974-5-20 0.0 1974-5-21 0.0 1974-5-22 1.0 1974-5-23 0.0 1974-5-24 0.0 1974-5-25 0.0 1974-5-26 0.0
4.53 15.3 9.12 5.27 3.74 3.17 2.83 2.89 2.55 2.24 2
0 1 2 3 4 5 6 7 8 9 10
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451 0.9031 0.9542 1.0000
1.1847 0.9600 0.7218 0.5729 0.5011 0.4518 0.4609 0.4065 0.3502 0.3010
Dataset 9 for figure (S) and (T) 1956-5-10 0.0 1956-5-11 30.5 1956-5-12 2.5 1956-5-13 7.1 1956-5-14 0.0 1956-5-15 6.1 1956-5-16 0.0 1956-5-17 0.0
4.62 17.6 8.83 7.48 5.44 4.93 4.56 4.13
0 1 2 3 4 5 6 7
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451
1.2455 0.9460 0.8739 0.7356 0.6928 0.6590 0.6160
Dataset 10 for figure (U) and (V) 1960-5-8 44.5 1960-5-9 19.1 1960-5-10 9.1 1960-5-11 5.6 1960-5-12 6.4 1960-5-13 2.0 1960-5-14 0.0
4.7 14 9.03 8.81 7.96 6.68 5.89
0 1 2 3 4 5 6
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782
1.1461 0.9557 0.9450 0.9009 0.8248 0.7701
vice (1972) and IHACRES (Identification of Unit Hydrographs and Component Flows from Rainfall, Evaporation, Stream flow data) developed by Jakeman et al. (1990), Jakeman and Hornberger (1993) and Jakeman et al. (1994) have been applied to predict annual runoff volume in the river system including ungauged basins from observed precipitation records (Cheng et al., 2004). The baseflow component can be separated from total river flow using integrated spatial-frequency analysis. This is essential for establishing a spatial-temporal stochastic on-line model for predicting river flow as multi-scale time series events in the study area (Cheng et al., 2004).
In order to validate the model (16) for fitting the shape of recessions starting at peak flows, we have selected five river gauging stations from the Oak Ridges Moraine (ORM) area (data from HYDAT CD-ROM User’s Manual, 1996). These stations have records of mean daily flow (m3/s) and daily rainfall data since 1950. The locations of the five stream flow gauging stations are shown in Fig. 1 superimposed on a digital elevation model (DEM) showing the landscape of the ORM area (Kenny, 1997). Peak flow events with a recession period longer than 4 days were selected for analysis. Altogether 10 events were selected from five stations and the dates, river flow (Q) and
164
Q. Cheng 100
A
B
100 Q(t) = 11.58e-0.198t, R2 = 0.839
Log Q(t)
Log Q(t)
Q(t) = 14.8t-0.883, R2 = 0.978
10
10
y = 11. 584e
- 0. 1 979 x
2
R = 0. 839
1
1 1
10
1
100
C
D
100 Log Q(t) = 1.531 + 0.06t – 1.77log(t) 2 – 0.40/t, R = 0.988
Log Q(t)
Q(t) = 2.013e2.258/t, R2 = 0.984
Log Q(t)
10
Log t (day)
Log t (day)
10
1
10
1 1
10
Log t (day)
1
10
Log t (day)
Figure 2 Results obtained from dataset 8 showing relationships between flow Q(t) and time t in log–log scale. Dots represent the observed data of log (Q) against log t and lines represent model fits using power-law (A), normal exponential (B), inverse time exponential (C), and the combined model (D), respectively. Log is 10-base logarithmic transformation. The parameters estimated from these models are shown in Table 2.
precipitation (P) are shown in Table 1. The data show that most river flow events reach flow peak after 1 day from the day that heavy rainfall occurs. After reaching the peaks there were no further significant rainfalls; therefore we can use model (16) assuming U(t) is close to zero. Taking the day with heavy rainfall as the time of reference, in most cases the peak flows occur on the first day which can be set as initial time t0. Three additional columns were added to Table 1 to represent the values of t, log (t) and log (Q). Four graphs were used to fit the data (Figs. 2–4) representing log (Q) vs. t, log (t), 1/t, and the combination of a0t + a1log (t) a2(1/t), respectively. The four graphs are used to compare the four models: power-law model (A), normal exponential model (B), exponential model with inverse time (C), and the combined model (D), respectively. The values of the parameters c, a0, a1, and a2 in the model (16) were estimated using multiple regression analysis with the data from the 10 flow events. The calculated results including the regression coefficients and multiple correlation coefficient R2, the correlation coefficient between log (Q) and the calculated value of log (Q) as c a0t a1 log (t) + a2t1, are shown in Table 2. The results illustrate that likelihood fitness for the combined model (16) is significantly better than those obtained with the exponential
model and power-law model alone. It can be concluded that the combined power-law and exponential model (16) is superior to the exponential model and power-law models on their own in characterizing the decay patterns of the peak flows at all chosen stations. Comparing the results in Table 2, one can see that: (1) in general, the parameter (a0) is relatively small (less than 0.3) implying that the contribution of the normal exponential is not very significant; (2) some datasets (datasets 2, 4, 6, and 8) show values of parameters, which are dominated by a1, implying that the model is power-law dominated. Fig. 2 shows the results obtained for dataset 8; (3) some datasets (datasets 3, 5, and 7) show values of parameters that are dominated by a2 implying that the model is dominated by the exponential form with inverse time. This can be seen from the results obtained from dataset 3 (Fig. 3); and (4) some datasets (1, 9, and 10) show values of parameters that are dominated by a1 and a2 implying that the model is dominated by both power-law and inverse time exponential forms. This can be seen from the results obtained from dataset 1 (Fig. 4). These variations might be related to the characteristics of the catchment as well as channel geometry. Further investigation is needed to explain the variations in more detail. The results estimated for those parameters involved in the
A combined power-law and exponential model for streamflow recessions 100
165
100 2
2
Q(t) = 27.52e-0.318t, R = 0.710
Log Q(t)
Log Q(t)
Q(t) = 34.9t-1.259, R = 0.918
10
10
1
1 1
1
10
10
Log t (day)
Log t (day) 100
100 Q(t) = 2.265e3.14/t, R2 = 0.997
Log Q(t) = 0.356 – 0.03log(t) + 1.37/t,
Log Q(t)
Log Q(t)
R2 = 0 998 10
1
10
1 1
10
Log t (day)
1
10
Log t (day)
Figure 3 Results obtained from dataset 3 showing relationships between flow Q(t) and time t in log–log scale. Dots represent the observed data of log (Q) against log t and lines represent model fits using power-law (A), normal exponential (B), inverse time exponential (C), and the combined model (D), respectively. Log is 10-base logarithmic transformation. The parameters estimated from these models are shown in Table 2.
new model using the ordinary linear multivariate regression gives negative values which may be prevented using special type of regression so that only positive regression coefficients can be guaranteed. This will be further investigated.
Concluding remarks It has been discussed in this paper that there are two ways in which non-linearity can be introduced into the flow differential equation. In sections 2 and 3 the implication and implementation of these two types of non-linearity extensions were discussed separately. These two types of non-linearity can be combined into a single equation to derive a solution of combined non-linearity. The differential equation can be expressed as follows: dQ ðtÞ 1 ¼ Q 2n ðtÞ dt nsðtÞ
ð17Þ
where n is the non-linearity index. Using the same approximation with the first three terms of Taylor series expansion applied to the coefficient function s(t) we can get 1 1 n1 Q ðtÞ ¼ c a0 t a1 logðtÞ þ a2 ð18Þ t
In comparison with model (16), model (18) includes an additional parameter (n). This model and its properties will be further discussed in future papers. Two aspects of the river flow recession model can be generalized to incorporate non-linearity: the most commonly used non-linear model (1) involves non-linearity (n) as power transformation applied to the flow, i.e. Qn(t). The other generalization proposed in the paper involves a non-constant monotonically decreasing coefficient function. For the non-linear differential equation with monotonically decreasing coefficient function the solution is a combination of normal exponential, powerlaw and exponential with inverse time forms, which was derived by approximating the coefficient function by the first three terms of a Taylor series expansion. In comparison with other existing multiple parameter flow recession models, the new model involves the combination of components (normal exponential, power-law and exponential with inverse time) which previously had been used as single models for describing flow recession processes. The physical interpretations of the parameters in these earlier models remain approximately valid when the new model is applied. The non-linear model can be expressed as a linear equation by logarithmic transformation and this makes it straightforward to implement the model and validate it graphically for observed river gauging data. The
166
Q. Cheng
A
100
B
100
2
2
Q(t) = 11.38e-0.236t, R = 0.668
Log Q(t)
Log Q(t)
Q(t) = 15.0t-1.027, R = 0.874
10
1
10
1 1
10
1
10
Log t (day)
C
Log t (day)
100
D
100
Q(t) = 1.479e2.736/t, R2 = 0.975
Log Q(t) = -1.198 +0.17t – 2.8log(t) + 2.75/t,
Log Q(t)
Log Q(t)
R2 = 0.995
10
1
10
1 1
10
1
10
Log t (day)
Log t (day)
Figure 4 Results obtained from Dataset 1 showing relationship between flow Q(t) and time t in log–log sacle. Dots represent the observed data of log (Q) against log t and lines represent model fits using power-law (A), normal exponential (B), inverse time exponential (C), and the combined model (D), respectively. Log is 10-base logarithmic transformation. The parameters estimated from these models are shown in Table 2.
Table 2
Results obtained by means of multiple regression
Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset
1 2 3 4 5 6 7 8 9 10
Constant C
Exponential term a0
Power-law term a1
Exponential term with reversetime a2
Multiple correlation coefficient R2
1.198 0.797 0.356 1.409 2.520 2.210 0.387 1.531 0.915 0.072
0.17 0.00 0.00 0.11 0.20 0.06 0.03 0.06 0.01 0.24
2.80 0.80 0.03 1.86 0.15 2.04 0.21 1.77 0.35 2.62
2.75 0.14 1.37 0.36 0.60 0.23 1.01 0.40 0.34 1.46
0.995 1.000 0.998 0.989 0.995 0.999 0.998 0.988 0.992 0.997
new four-parameter model has been validated with river peak flow data from the ORM area, Canada. The results obtained from 10 datasets from the ORM have clearly demonstrated that the combined model provides better likelihood fitness to the observed flow data.
It should be remembered that the paper deals only with flow recessions in the absence of rainfall. More general results related to continuous flow simulation from rainfall and model validation with simulation will be further investigated.
A combined power-law and exponential model for streamflow recessions
Acknowledgements Thanks are due to Dr. F. P. Agterberg at the Geological Survey of Ottawa, Canada for his constructive comments and improvement of the representation of the manuscript. Four anonymous reviewers are thanked for their constructive comments on reviewing the earlier version of the manuscript. This research has been jointly supported by NSERC Discovery Research ‘‘development of GIS technology for mineral exploration and environmental assessments’’ (ERC-OGP0183993) and Chinese National Foundation of Science Projects ‘‘Distinguished Young Researcher Grant’’ (40525009) and Strategic Research project (40638041).
References Beven, K.J., 2000. Rainfall–runoff Modeling: The Primer. John Wiley and Sons, Chichester, UK, p. 360. Boussinesq, J., 1904. Recherches theoretique sur l’ecoulement des nappes d’eau infiltrees dans le sol et sur le debit des sources. J. Math. Pure Appl. 10 (5th series), 5–78. Brutsaert, W., Nieber, J.L., 1977. Regionalized drought flow hydrographs from a mature glaciated plateau. Water Resour. Res. 13 (3), 637–643. Cheng, Q., 2004. Weights of evidence modeling of flowing wells in the Greater Toronto Area, Canada. Nat. Resour. Res. 13 (2), 77– 86. Cheng, Q., Russell, H., Sharpe, D., Kenny, F., Pin, Q., 2001. GISbased statistical and fractal/multifractal analysis of surface stream patterns in the Oak Ridges Moraine. Comput. Geosci. 27 (5), 513–526. Cheng, Q., Zhang, G., Lu, C., Ko, C., 2004. GIS spatial-temporal modeling of water systems in the Greater Toronto Area, Canada, Earth Sciences. J. China Univ. Geosci. 15 (3), 275–282. Cheng, Q., Ko, C., Yuan, Y., Ge, Y., Zhang, S., 2006. GIS Modeling for predicting river runoff volume in ungauged drainages in the Greater Toronto Area, Canada. Comput. Geosci. 32 (8), 1108– 1119. Coutagne, A., 1948. Etude generale des variations de debits en function des facteurs qui les conditionnent, 2eme partie, les variations de debit en periode non influence par precipitations’, La Houille Blanche, Sept. Oct., pp. 416–436. Croke, B.F.W., 2006. A technique for deriving an average event unit hydrograph from streamflow – only data for ephemeral quickflow-dominant catchments. Adv. Water Resour. 29, 493–502. Hall, F.R., 1968. Base flow recession – a review. Water Resour. Res. 4 (5), 973–983. Horton, R.E., 1933. The role of infiltration in the hydrologic cycle. Trans. Am. Geophys. Union 14, 446–460. HYDAT CD-ROM User’s Manual, 1996. Surface Water and Sediment Data, Atmospheric Environment Program, Version 96 – 1.04 User’s Manual, Environment Canada, 95p. Ishihara, T., Takagi, F., 1965. A study of the variation of low flow. Disaster Prev. Res. Inst. Kyoto Univ. Bull. 15 (2), 75–98. Jakeman, A.J., Hornberger, G.M., 1993. How much complexity is warranted in a rainfall–runoff model? Water Resour. Res. 29 (8), 2637–2649.
167
Jakeman, A.J., Littlewood, I.G., Whitehead, P.G., 1990. Computation of the instantaneous unit hydrograph and identifiable component flows with application to two small upland catchments. J. Hydrol. 117, 275–300. Jakeman, A.J., Post, D.A., Beck, M.B., 1994. From data and theory to environmental model: the case of rainfall–runoff. Environmetrics 5, 297–314. Kenny, F.M., 1997. A chromostereo enhanced digital elevation model of the Oak Ridges Moraine Area, Southern Ontario and Lake Ontario Bathymetry, Geological Survey of Canada, Open file 3423, Scale 1:200,000. Ko, C., Cheng, Q., 2004. GIS spatial modeling of river flows and precipitation in the ORM, Ontario. Comput. Geosci. 30 (4), 379– 389. Kokkonen, T.S., Jakeman, A.J., Young, P.C., Koivusalo, H.J., 2003. Predicting daily flows in ungauged catchments: model regionalization from catchment descriptors at the Coweeta Hydrologic Laboratory, North Carolina. Hydrol. Process. 17, 2219–2238. Littlewood, I.G., 2002. Improved unit hydrograph characterisation of the daily flow regime (including low flows) for the River Teifi, Wales: towards better rainfall–streamflow models for regionalisation. Hydrol. Earth Syst. Sci. 6 (5), 899–911. Littlewood, I.G., Down, K., Parker, J.R., Post, D.A., 1997. IHACRES v. 1.0 User Guide. Center for Ecology and Hydrology, Wallingford, UK & Integrated Catchment Assessment and Management Center, Australian National University, Canberra, 94pp. Maillet, E., 1905. Essai d’hydraulique souterraine et fluviale, Libraire Sci., A, Herman, Paris (Cited by Hall, 1968). Mitchell, W.D., 1972. Model hydrographs. US Geol. Surv. Water Supply Pap. 2005. Moore, R.D., 1997. Storage-outflow modeling of streamflow recession, with application to a shallow-soil forested catchment. J. Hydrol. 198, 260–270. Napiorkowski, J.J., Strupczewski, W.G., 1979. The analytical determination of the Kernels of the Volterra Series describing the cascade of nonlinear reservoirs. J. Hydrol. Sci. 6 (3-4), 121– 142. Sefton, C.E.M., Howarth, S.M., 1998. Relationships between dynamic response characteristics and physical descriptors of catchments in England and Wales. J. Hydrol. 211, 1–16. Tallaksen, L.M., 1995. A review of baseflow recession analysis. J. Hydrol. 165, 349–370. US Department of Agriculture, Soil Conservation Service, 1972. National Engineering Handbook, Section 4, Hydrology, US Government of Printing Office, Washington, DC, 544p. Wagener, T., Wheater, H., Gupta, H.V., 2004. Rainfall–runoff Modeling in Gauged and Ungauged Catchments. Imperial College Press, London, 306p. Werner, P.W., Sundquist, K.J., 1951. On the groundwater recession curve for large watersheds, IAHS General Assembly, Brussels, IAHS Publ., 33, pp. 202–212. Wittenberg, H., 1999. Baseflow recession and recharge as nonlinear storage processes. Hydrol. Process. 13, 715–726. Yuan, Y., Cheng, Q., 2005. Predicting streamflow in ungauged basins with GIS and rainfall–runoff models: a case study of model regionalization in the Oak Ridges Moraine area, Southern Ontario, Canada. In: Qiuming, Cheng, Graeme, Bonham-Carter (Eds.), Proc. of the Int. Ass. for Math. Geol. Annual Conference, 1, pp. 428–433.
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jhydrol
A combined power-law and exponential model for streamflow recessions Qiuming Cheng
*
Department of Earth and Space Science and Engineering, Department of Geography, York University, 4700 Keele Street, Toronto, Canada M3J 1P3 State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan and Beijing, China Received 6 November 2007; accepted 9 January 2008
KEYWORDS Non-linear process; Fractals; River flow; Power-law; Exponential model; Combined model
Summary This paper proposes a non-linear model to describe streamflow recessions. A model of combining power-law and exponential functions is derived from a non-linear differential equation associating outflow change rate and flow itself with a time-dependent decreasing decay rate. The decay rate is approximated by the first three terms in a Taylor series expansion. The solution of the resulting differential equation includes three components of the expansion determining two exponential functions and a power-law function, respectively. It is demonstrated analytically and empirically that the introduction of a non-linear differential equation comprising a mixture of power-law and exponential forms is superior to using an ordinary differential equation with exponential solution, or a nonlinear equation with power-law solution, for modeling the decay patterns of flow in river systems. Analytical properties of the mixing model are discussed and the model is subsequently validated by characterizing the relationships between decay of peak flow as a function of time after the initial occurrence of flow peak in 10 datasets of flow events chosen from five gauging stations in the Oak Ridges Moraine (ORM) area, Canada. ª 2008 Elsevier B.V. All rights reserved.
Introduction Quantitative description of the flow recession process in a river system and its relationship with precipitation is one * Address: Department of Earth and Space Science and Engineering, Department of Geography, York University, 4700 Keele Street, Toronto, Canada M3J 1P3. E-mail address: [email protected]
of the most important tasks of contemporary hydrology. As an important surface water system, river flow has significant influence on many other environmental systems. Understanding river flow recession processes is not only scientifically important but also essential for water resource management and geo-hazard (flooding, landslide, etc.) prevention. Therefore, river gauging and flow monitoring have being conducted worldwide for many river systems. These programs have generated a large amount of data which
0022-1694/$ - see front matter ª 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2008.01.017
158
Q. Cheng
can be used for empirical or analytical modeling of the properties of river systems. Detailed study of river flow recession processes has been increasingly attracting attention from researchers in academic institutions, industries and governments at various levels. Prediction of runoff in river systems is also important for the insurance industry, transportation and other relevant sectors. To develop a practical system for prediction of river flow requires several new developments including sophisticated models, calibration–validation techniques, model generalization and an efficient on-line system to disseminate the results publicly. A number of mathematical models have been proposed in the literature for quantifying the recession process of river flow. Some of these models, for example, rainfall–runoff models (RRMs), have become standard tools that are routinely used for hydrological modeling in engineering and environmental science (Littlewood et al., 1997; Sefton and Howarth, 1998; Wagener et al., 2004; Littlewood, 2002; Kokkonen et al., 2003; Yuan and Cheng, 2005). However, an optimum model for prediction purposes in all types of catchment environments at various scales does not yet exist. While most authors use a constant-coefficient exponential model derived from a linear differential equation, nonlinear models have also been investigated and applied. It is generally accepted that linear models are suitable for simple catchment environments but non-linear models might be the most appropriate tool for more complex catchment systems. There are two kinds of recession models in which non-linearity has been introduced: the best known form of non-linearity is related to the basic non-linear differential equation governing unsteady flow from a large unconfined aquifer in a stream channel presented by Boussinesq in 1877 (Hall, 1968; Tallaksen, 1995; Brutsaert and Nieber, 1977). In this context, recession rate of flow Q(t) can be associated with storage S(t) as SðtÞ ¼ sQ n ðtÞ
ð1Þ
where s is the storage coefficient and n is the non-linearity parameter. Main cause of the non-linearity of the relationship between rainfall and runoff is the antecedent moisture condition (Beven, 2000; Tallaksen, 1995; Wittenberg, 1999). If n = 1, the model becomes the most commonly used linear form, which has the advantage of computational efficiency. Another mathematical description of flow systems is that changes in storage S(t) of a hydrological system in time depend on inflow U(t) and outflow Q(t) that are related by dSðtÞ ¼ UðtÞ Q ðtÞ dt
ð2Þ
In dealing with river flow mainly caused by rainfall the U(t) can be considered as effective rainfall, being the contribution of rainfall to the river flow. If there is no inflow, U(t) = 0, or if the amount of U(t) is so small that it can be neglected, then Eqs. (1) and (2) can be combined as dQ ðtÞ 1 ¼ Q 2n ðtÞ dt sn
ð3Þ
This approach results in different types of non-linearity depending on the value of n: (1) When n 5 1, the non-linear differential equation can be solved and its solution is
1
Q ðtÞ ¼ Q ðt0 Þ½1 þ cðt t0 Þn1
ð4Þ
Q 1n ðt0 Þ is a constant and Q(t0) is the value of where c ¼ n1 ns Q(t) at initial time t0. This result is similar to that given by Brutsaert and Nieber (1977), Coutagne (1948) and Moore (1997). The relation between log Q(t) and log [1 + c(t t0)] can be plotted as a straight-line on log–log paper. When n = 1/2, the model (4) generates a hyperbola form as has been discussed by many authors (Boussinesq, 1904; Hall, 1968; Brutsaert and Nieber, 1977). (2)
When n = 1, Eq. (3) has the exponential solution 1 Q ðtÞ ¼ Q ðt0 Þ exp ðt t0 Þ ð5Þ s where 1/s is a time independent parameter determining the decay rate of runoff Q(t) without effective rainfall (Maillet, 1905; Mitchell, 1972). Q(t0) represents initial river flow at time t0, and Q(t) is river flow at time t. Since the parameter s is the time it takes for the value Q(t0) to recess to exp (1) of its original value, Q(t0) can be considered as 1 for unit hydrograph. The function (5) can be plotted on semi-logarithmic paper as straight-line with slope equal to 1/s. This is the most commonly used model for describing hydrograph recessions and river peak flow recession. However, a number of studies have shown that the function (5) is often not sufficient to represent the relationship because the plot shows curvature and this has led to the investigation of other, non-linear, models. The other type of generalization of non-linearity is based on modification of the coefficient s in model (3). Croke (2006) has recently pointed out that the assumption that the unit hydrograph has a constant-coefficient exponential form is not valid from a conservation of mass point of view. He suggested a power-law with monotonically decreasing coefficient as an alternative form for the unit hydrograph. Cheng et al. (2006) noticed this suggesting the use of a differential equation with time-decreasing coefficient yielding a power-law form for describing decay patterns of flow. From these points of view the constantcoefficient exponential model only represents a special solution for catchment channel modeling with limited applicability. In Croke (2006), a power-law model was proposed to describe river flow without effective rainfall and it was implemented using a sum of multiple exponential terms for approximation. That author also pointed out that the exponential terms used for representing the power-law do not necessarily have any physical significance. A similar model involving a sum of multiple exponential terms had already been used by others including Werner and Sundquist (1951). Other non-linearity treatments include Napiorkowski and Strupczewski’s (1979) cascade of non-linear reservoirs described by a Volterra series. Ishihara and Takagi (1965) modeled recession flow as a sum of two components: outflow from a confined aquifer subject to a simple exponential decay (5), and outflow from an unconfined aquifer expressed by means of a hyperbola form (4). It can be seen from the above review that although exponential and power-law models are two most commonly used models, application of either of these is often not sufficient for describing the behavior of river flow recession. Some efforts have been made to investigate
A combined power-law and exponential model for streamflow recessions combined multiple parameter models, for example using a sum of multiple exponential terms, but these types of combinations of the same types of functions are neither physically based nor straightforward for implementation. The current paper attempts to introduce a new mixing model combining exponential and power-law forms derived from a non-linear differential equation with a variable coefficient function that is approximated by the first three terms of a Taylor series expansion. This new combined model is demonstrated to be appropriate for characterizing streamflow recessions in various types of catchment. Comparisons between the power-law, exponential and combined models were conducted both analytically and empirically using river flow datasets for five gauging stations selected from the Greater Toronto Area, Canada.
Non-linear model for describing river peak flow recession In order to discuss the two suggested non-linear extensions to the flow model, we will start with Eq. (3) with the coefficient s as a function of time t: dQ ðtÞ 1 ¼ Q ðtÞ dt sðtÞ
ð6Þ
In order to solve Eq. (6), one needs to know the functional form of s(t). It is usually not possible to know the exact expression of this function since it may have different forms depending upon the dominant system of river flow and the characteristics of the catchments channels. For example, in some simple cases if the change of flow rate dQdtðtÞ is mainly due to magnitude of Q(t) or if the flow is dominated by direct surface flow in small catchments with simple regular channel geometries, the coefficient should be independent of the time duration after the peak flow. In other catchments, however, decreasing recession rate may not only be proportional to the flow itself but may also decrease through time due to channel geometry changes; for example, when flow in the channel reduces, wetted channel cross-sectional area and saturation of the catchment basin reduce as well. In this case, the coefficient might be expressed as 1/[s(t)] = 1/[st] which reduces monotonically as t increases. The solution of the flow differential equation assumes the power-law form Q ðtÞ ¼ Q ðt0 Þ
1=s t t0
ð7Þ
More complex situations might be applicable to actual river channel systems because wetted channel cross-sectional area and saturation of the catchment basin reduce as Q(t) decreases through time t. This situation can be characterized by 1/[s(t)] = 1/[stm] where m is a positive constant number (m > 1) which yields the solution 1 ðt1m t1m Q ðtÞ ¼ Q ðt0 Þ exp Þ ð8Þ 0 sð1 mÞ In order to understand the properties of this solution for different non-linearity (m), we will determine the ratio of flow for two different time periods.
159
(1) The case of m = 0. When m = 0, the solution satisfies the simple exponential form (5) with time interval t t0. At time interval t t0 = s, Q(t) recesses to Q(t0)exp (1) with the following properties: (a) When t ! 1, then Q(t) ! 0. (b) The ratio of Q(t) at two different points of time, say t1 and t2, is Q ðt1 Þ 1 ¼ exp ðt1 t2 Þ ð9Þ Q ðt2 Þ s This property indicates that the ratio of Q(t1) and Q(t2) depends only on the time interval (t1 t2) and is independent of t1 and t2 (t2 < t1). For a given time interval, for example, t1 t2 = k where k is a constant, the ratio of Q(t1) and Q(t2) will remain constant no matter how much time has expired after the initial time t0 (e.g. if t1 and t2 t0). The parameter 1/s = [log Q(t1) log Q(t2)]/[t1 t2], represents the ratio of log difference between Q(t1) and Q(t2) and the difference between t1 and t2. (2) The case of m = 1. When m = 1 in (6), the decay rate of Q(t) depends not only on Q(t) itself but also on 1/[st] implying the decay rate decreases with time t. The decay of Q(t) follows a power-law relationship with time t. For time ratio t/ t0 = ks, Q(t) is reduced by the factor k to Q(t0)/k with the following properties: (a) When t/t0 ! 1, then Q(t) ! 0. (b) The ratio of Q(t) at two different points of time, say t1 and t2, is 1s Q ðt1 Þ t1 ¼ ð10Þ Q ðt2 Þ t2 This property indicates that the ratio Q(t1)/Q(t2) is proportional to the ratio (t1/t2)1/s, and the parameter 1/s = [log Q(t1) log Q(t2)]/[log t1 log t2], represents the ratio of log difference between Q(t1) and Q(t2) and the log difference between t1 and t2. For comparison purposes we can check the ratio of Q(t2 + k)/Q(t2) at the two different times (t1 and t2) with constant time interval t1 t2 = k. This illustrates that, for any given constant k, when t1 = t2+k ! 1, Q(t1) tends to become the same as Q(t2). (c) According to (10), if t1/t2 = 1/k (constant), then Q(t1)/ Q(t2) = k1/s. For example, if k = 2, then Q(t)/Q(2t) = 21/s, implying that in order to maintain the same reduced rate from Q(t) to Q(2t) one needs to keep constant ratio of t1/t2 = 1/2. Then t2 = 2 days if t1 = 1 day after the initial time t0; and t2 = 10 days if t1 = 5 days after the initial time t0. This indicates that the decay becomes slower moving away from the initial time. (3) The case of m > 1. When m > 1, the solution (8) shows that the decay of Q(t) follows an exponential relationship with time interval t1m t1m . For time interval t1m t1m ¼ sð1 mÞ, Q(t) 0 0 is reduced to simpler forms with the following properties:
160
Q. Cheng
h i t1m 0 (a) When t ! 1, then Q ðtÞ ! Q ðt0 Þ exp sðm1Þ . (b) The ratio of Q(t) at two different points of time, say t1 and t2, is 1m Q ðt1 Þ 1 t1 t1m ð11Þ ¼ exp 2 Q ðt2 Þ sð1 mÞ This property indicates that the ratio of the Q(t) at different times is proportional to the difference t1m t1m , and 1 2 the parameter 1=s ¼ ð1 mÞ½log Q ðt1 Þ log Q ðt2 Þ=½t1m 1 t1m , represents the ratio of log difference between 2 Q(t1) and Q(t2) and (1 m) times the difference between t1m and t1m . The form (11) is the same as the non-linear 1 2 relationship proposed by Horton (1933). It has been shown that the model (6) involves a parameter m determining the non-linearity of the system. The ordinary system with constant-coefficient s becomes the special case of m = 0. Given that the different values of m indicate different decay rates which might be dominated by different types of flow and drainage basins as proposed by many authors, we can explore mixed flow systems by combining the functions with m = 0, 1 and 2 in the preceding expression. Taylor series expansion can be applied to 1/[s(1/t*)] where t*=1/t. In this paper the coefficient function 1/s(t) is approximated by means of Taylor series expansion of 1/[s(1/t*)] resulting in 2 m 1 1 1 1 1 1 1 þ c2 þ þ o ¼ c0 þ c1 sðtÞ t t0 t t0 t t0 ð12Þ where t0 must be positive and t usually is greater than 1 (t0 > 0 and t > 1) so that 1/t 1/t0 is usually negative and much smaller than 1 which guarantees convergence of the summation. The term o( ) is the remainder term and states that there exists a small error for the approximation. If we take the first three terms for approximation, then the form (12) can be rearranged into 1 1 1 1 ¼ a0 þ a1 þ a2 2 þ o 2 ð13Þ sðtÞ t t t where the components a0, a1, and a2 correspond to the coefficients of the first three orders of 1/t. The remainder term states that the error for the approximation approaches zero when t increases infinitely and the order of the error is smaller than the order of (1/t2). Using this approximation one can express the differential model as dQ ðtÞ 1 1 ¼ a0 þ a1 þ a2 2 Q ðtÞ ð14Þ dt t t It can be seen that the Eq. (14) is a combination of Eq. (6) with m = 0, 1 and 2, respectively. Therefore, the model (14) represents the combination of three different types of forms: normal exponential, power-law and exponential of inverse time i a1 h a0 ðtt0 Þa2 1t t1 t 1 0 Q ðtÞ ¼ Q ðt0 Þ ¼ cta1 ea0 tþa2 t e t0 ð15Þ The model (15) gives a new four-parameter model with the parameter c as constant and a0, a1, and a2 representing the coefficients of the first three increment terms,
respectively. In the model (15) both the first term, exp (a0t), and the third term, exp (a2/t), are monotonically decreasing functions of time t representing the decay regulations of flow Q(t) as time increases. An advantage of model (15) in comparison with other multiple parameter non-linear models is that it can be converted into a simple form by applying logarithmic transformation to both sides to give log Q ðtÞ ¼ c a0 t a1 logðtÞ þ a2 t1
ð16Þ
Using this expression, it is straightforward to validate the model graphically and to estimate the parameters using multiple regression. The parameters involved in the combined model can be easily estimated by means of regression applied to log Q(t) against t, log t and 1/t. The multiple correlation coefficient (R2) then can be used as a measure of overall likelihood fitness. This model will be implemented and validated using daily flow data from five gauging stations from the Oak Ridge Moraine area, Canada.
Validation of the new combined power-law and exponential model for characterizing flow recessions The study area on the Oak Ridges Moraine (ORM) is located in southern Ontario, Canada. This is one of the most developed areas in Canada. The Oak Ridges Moraine is one of the most significant landforms in southern Ontario. The moraine is characterized by its rolling hills and river valleys extending 160 km from the Niagara Escarpment to Rice Lake and was formed 12,000 years ago by advancing and retreating glaciers. The moraine contains the headwaters of 65 river systems and has a wide diversity of streams, woodlands, wetlands, kettle lakes, kettle bogs and significant flora and fauna. It is one of the last remaining continuous green corridors in southern Ontario containing 30% forested area. The moraine’s sands and gravel deposits absorb rain and snow melt. The underground water is then stored in layers of sand and gravel (aquifers), filtered and slowly released as cool fresh water to the 65 rivers and streams flowing north into Lakes Simcoe and Scugog and south into Lake Ontario. The greatest threat to the function of the moraine is inappropriate land-use at the surface of the moraine, particularly in headwater areas. Situated on the north shore of Lake Ontario, the Oak Ridges Moraine area experiences a wide range of weather conditions through an average year. Residents of the area have come to expect bitter cold winter and oppressive summer heat. Precipitation includes deep winter snow, heavy spring and summer thunderstorms, and, sometimes, lengthy summer dry spells. The temperature is moderate most of the time and, annually, there is approximately 600– 1000 mm of precipitation including snow in winter season and rainfall in other seasons across the ORM area. Periods of wet or dry weather are not uncommon, often making a year significantly different from normal. Snowfall mainly in December to March accounts for about 15% of the total annual precipitation. In the summer, widespread daily-long rainfalls are rare. Instead, showery precipitation is the gen-
A combined power-law and exponential model for streamflow recessions eral rule. This frequently results in precipitation varying widely from place to place. In addition, this kind of precipitation may come as short but intense downpours that quickly run into drains and streams. A thunderstorm that results in 20 or 30 mm of rain in just one hour does little to relieve the moisture deficit that comes with a long period of dry weather. On occasion, this area receives the effects of the remnants of a tropical storm or hurricane moving through the eastern part of the United States. Occurring from midsummer to autumn, these systems may result in an allday rainfall that can be quite heavy. Precipitation amounts from these systems can also be somewhat variable over relatively short distances. When Hurricane Hazel hit the area in October 1954, 181.6 mm of rain was recorded at Snelgrove station: however, only 38.1 mm at Alton station just 22 km northwest of Snelgrove at the same time. Land-use in the ORM area is largely rural, with forest and agricultural practices dominating the landscape although it includes the Greater Toronto Area that is a highly industrialized area. Increasing population and residential development in this area has led to residential development expanding into the rural areas. Land-use changes, primarily the building of residential subdivisions, the construction of roads and the paving of parking lots,
02HB001
161
increase the imperviousness of the ground surface. Consequently, this surface runoff results in dramatic increases in wet weather flows of the headwater streams on the Moraine causing erosion and degradation of these fragile systems. Various modeling exercises conducted on the study area have demonstrated that runoff volumes observed at the river gauging stations in the drainage basin networks are highly correlated with the baseflow caused by ground water discharge, snow melting, direct flow caused by rainfall as well as the physical properties of drainage basins (Cheng et al., 2006). It has been found that the delay response of surface runoff (river peak flow) to heavy precipitation events varies across the study area depending on the physical properties of drainage networks and stream systems such as the complexity of drainage basin boundaries, glacier sediment types, land-use types, and surface inclination (Ko and Cheng, 2004). It also has been demonstrated that long-term persistency of river flow is related to the characteristics of drainage basins (Ko and Cheng, 2004). The interaction between groundwater and surface runoff also has been investigated so that the influence of groundwater discharge on river flow could be modeled (Cheng et al., 2001, 2006; Cheng, 2004). The conventional SCS method (Soil Conservation Service) of the US Department of Agriculture, Soil Conservation Ser-
DCB
HRB 02HC003
02HC009
CRB
DCB – Duffins Creek Basin 02HC003
02HB002
HRB - Humber River Basin CRB - Credit River Basin
Figure 1 Shaded relief of the digital elevation model of the ORM (DEM data from Kenny, 1997). Lines represent the major basins of the area. White dots represent the locations of five river gauging stations chosen for the study. Labels in the boxes are the IDs of gauging stations.
162 Table 1 Date
Q. Cheng Peak river flow and precipitation observed from five gauging stations in the ORM Precipitation P (mm)
River flow (Q) (m3/s)
t (Day)
Log (t)
Log (Q)
0 1 2 3 4 5 6 7 8 9
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451 0.9031 0.9542
1.3909 0.6590 0.5977 0.4564 0.4564 0.4298 0.3560 0.3096 0.2455
0 1 2 3 4 5
0.0000 0.3010 0.4771 0.6021 0.6990
0.9420 0.6314 0.4564 0.3598 0.2648
Pickering Station on Duffins Creek 02HC006 (201 km2) Dataset 3 for figure (E) and (F) 1954-10-15 78.7 35.7 1954-10-16 0.8 51.8 1954-10-17 1.0 11.5 1954-10-18 0.0 6.17 1954-10-19 0.0 4.84 1954-10-20 0.0 4.28 1954-10-21 0.0 3.96
0 1 2 3 4 5 6
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782
1.7143 1.0607 0.7903 0.6848 0.6314 0.5977
Dataset 4 for figure (G) and (H) 1960-5-12 10.7 1960-5-13 0.0 1960-5-14 0.0 1960-5-15 0.0 1960-5-16 0.0 1960-5-17 0.0
11.4 14.4 7.28 5.64 4.16 3.37
0 1 2 3 4 5
0.0000 0.3010 0.4771 0.6021 0.6990
1.1584 0.8621 0.7513 0.6191 0.5276
Weston Station on East Humber River 02HC003 (820 km2) DataSet 5 for figure (I) and (J) 1954-10-15 125.3 3.79 1954-10-16 1.3 838 1954-10-17 1.3 498 1954-10-18 0 155 1954-10-19 0 87.8 1954-10-20 0 63.7 1954-10-21 0 44.2
0 1 2 3 4 5 6
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782
2.9232 2.4972 2.1903 1.9435 1.8041 1.6454
Erindale Station on Credit River 02HB002 (806 km2) Dataset 6 for figure (M) and (N) 1954-10-15 129.6 1954-10-16 2.3 1954-10-17 1.4 1954-10-18 0 1954-10-19 0 1954-10-20 0 1954-10-21 0 1954-10-22 0
0 1 2 3 4 5 6 7
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451
2.5011 1.8021 1.5092 1.2923 1.0792 1.0000 0.9117
2
Pine Grove Station on East Humber River 02HC009 (260 km ) Dataset 1 for figure (A) and (B) 1956-5-11 23.1 1.95 1956-5-12 3.3 24.6 1956-5-13 3.7 4.56 1956-5-14 0 3.96 1956-5-15 6.8 2.86 1956-5-16 0.9 2.86 1956-5-17 3.1 2.69 1956-5-18 1.9 2.27 1956-5-19 0 2.04 1956-5-20 0 1.76 Dataset 2 for figure (C) and (D) 1960-5-23 11.9 1960-5-24 1.3 1960-5-25 2.5 1960-5-26 0.0 1960-5-27 0.0 1960-5-28 0.0 1960-5-29 0.0
3.23 8.75 4.28 2.86 2.29 1.84 3.48
7.08 317 63.4 32.3 19.6 12 10 8.16
A combined power-law and exponential model for streamflow recessions
163
Table 1 (continued) Date
River flow (Q) (m3/s)
Precipitation P (mm)
t (Day)
Log (t)
Log (Q)
0 1 2 3 4 5 6 7
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451
1.3674 0.9009 0.7284 0.6160 0.6010 0.5378 0.4771
2
Cataract Station on Credit River 02HB001 (174 km ) Dataset 7 for figure (O) and (P) 1954-10-15 69.6 10.2 1954-10-16 1.4 22.2 1954-10-17 3 23.3 1954-10-18 0 7.96 1954-10-19 0 5.35 1954-10-20 0 4.13 1954-10-21 0 3.99 1954-10-22 0 3.45 1954-10-23 0 3 Dataset 8 for figure (Q) and (R) 1974-5-16 29.7 1974-5-17 0.0 1974-5-18 0.0 1974-5-19 0.0 1974-5-20 0.0 1974-5-21 0.0 1974-5-22 1.0 1974-5-23 0.0 1974-5-24 0.0 1974-5-25 0.0 1974-5-26 0.0
4.53 15.3 9.12 5.27 3.74 3.17 2.83 2.89 2.55 2.24 2
0 1 2 3 4 5 6 7 8 9 10
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451 0.9031 0.9542 1.0000
1.1847 0.9600 0.7218 0.5729 0.5011 0.4518 0.4609 0.4065 0.3502 0.3010
Dataset 9 for figure (S) and (T) 1956-5-10 0.0 1956-5-11 30.5 1956-5-12 2.5 1956-5-13 7.1 1956-5-14 0.0 1956-5-15 6.1 1956-5-16 0.0 1956-5-17 0.0
4.62 17.6 8.83 7.48 5.44 4.93 4.56 4.13
0 1 2 3 4 5 6 7
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782 0.8451
1.2455 0.9460 0.8739 0.7356 0.6928 0.6590 0.6160
Dataset 10 for figure (U) and (V) 1960-5-8 44.5 1960-5-9 19.1 1960-5-10 9.1 1960-5-11 5.6 1960-5-12 6.4 1960-5-13 2.0 1960-5-14 0.0
4.7 14 9.03 8.81 7.96 6.68 5.89
0 1 2 3 4 5 6
0.0000 0.3010 0.4771 0.6021 0.6990 0.7782
1.1461 0.9557 0.9450 0.9009 0.8248 0.7701
vice (1972) and IHACRES (Identification of Unit Hydrographs and Component Flows from Rainfall, Evaporation, Stream flow data) developed by Jakeman et al. (1990), Jakeman and Hornberger (1993) and Jakeman et al. (1994) have been applied to predict annual runoff volume in the river system including ungauged basins from observed precipitation records (Cheng et al., 2004). The baseflow component can be separated from total river flow using integrated spatial-frequency analysis. This is essential for establishing a spatial-temporal stochastic on-line model for predicting river flow as multi-scale time series events in the study area (Cheng et al., 2004).
In order to validate the model (16) for fitting the shape of recessions starting at peak flows, we have selected five river gauging stations from the Oak Ridges Moraine (ORM) area (data from HYDAT CD-ROM User’s Manual, 1996). These stations have records of mean daily flow (m3/s) and daily rainfall data since 1950. The locations of the five stream flow gauging stations are shown in Fig. 1 superimposed on a digital elevation model (DEM) showing the landscape of the ORM area (Kenny, 1997). Peak flow events with a recession period longer than 4 days were selected for analysis. Altogether 10 events were selected from five stations and the dates, river flow (Q) and
164
Q. Cheng 100
A
B
100 Q(t) = 11.58e-0.198t, R2 = 0.839
Log Q(t)
Log Q(t)
Q(t) = 14.8t-0.883, R2 = 0.978
10
10
y = 11. 584e
- 0. 1 979 x
2
R = 0. 839
1
1 1
10
1
100
C
D
100 Log Q(t) = 1.531 + 0.06t – 1.77log(t) 2 – 0.40/t, R = 0.988
Log Q(t)
Q(t) = 2.013e2.258/t, R2 = 0.984
Log Q(t)
10
Log t (day)
Log t (day)
10
1
10
1 1
10
Log t (day)
1
10
Log t (day)
Figure 2 Results obtained from dataset 8 showing relationships between flow Q(t) and time t in log–log scale. Dots represent the observed data of log (Q) against log t and lines represent model fits using power-law (A), normal exponential (B), inverse time exponential (C), and the combined model (D), respectively. Log is 10-base logarithmic transformation. The parameters estimated from these models are shown in Table 2.
precipitation (P) are shown in Table 1. The data show that most river flow events reach flow peak after 1 day from the day that heavy rainfall occurs. After reaching the peaks there were no further significant rainfalls; therefore we can use model (16) assuming U(t) is close to zero. Taking the day with heavy rainfall as the time of reference, in most cases the peak flows occur on the first day which can be set as initial time t0. Three additional columns were added to Table 1 to represent the values of t, log (t) and log (Q). Four graphs were used to fit the data (Figs. 2–4) representing log (Q) vs. t, log (t), 1/t, and the combination of a0t + a1log (t) a2(1/t), respectively. The four graphs are used to compare the four models: power-law model (A), normal exponential model (B), exponential model with inverse time (C), and the combined model (D), respectively. The values of the parameters c, a0, a1, and a2 in the model (16) were estimated using multiple regression analysis with the data from the 10 flow events. The calculated results including the regression coefficients and multiple correlation coefficient R2, the correlation coefficient between log (Q) and the calculated value of log (Q) as c a0t a1 log (t) + a2t1, are shown in Table 2. The results illustrate that likelihood fitness for the combined model (16) is significantly better than those obtained with the exponential
model and power-law model alone. It can be concluded that the combined power-law and exponential model (16) is superior to the exponential model and power-law models on their own in characterizing the decay patterns of the peak flows at all chosen stations. Comparing the results in Table 2, one can see that: (1) in general, the parameter (a0) is relatively small (less than 0.3) implying that the contribution of the normal exponential is not very significant; (2) some datasets (datasets 2, 4, 6, and 8) show values of parameters, which are dominated by a1, implying that the model is power-law dominated. Fig. 2 shows the results obtained for dataset 8; (3) some datasets (datasets 3, 5, and 7) show values of parameters that are dominated by a2 implying that the model is dominated by the exponential form with inverse time. This can be seen from the results obtained from dataset 3 (Fig. 3); and (4) some datasets (1, 9, and 10) show values of parameters that are dominated by a1 and a2 implying that the model is dominated by both power-law and inverse time exponential forms. This can be seen from the results obtained from dataset 1 (Fig. 4). These variations might be related to the characteristics of the catchment as well as channel geometry. Further investigation is needed to explain the variations in more detail. The results estimated for those parameters involved in the
A combined power-law and exponential model for streamflow recessions 100
165
100 2
2
Q(t) = 27.52e-0.318t, R = 0.710
Log Q(t)
Log Q(t)
Q(t) = 34.9t-1.259, R = 0.918
10
10
1
1 1
1
10
10
Log t (day)
Log t (day) 100
100 Q(t) = 2.265e3.14/t, R2 = 0.997
Log Q(t) = 0.356 – 0.03log(t) + 1.37/t,
Log Q(t)
Log Q(t)
R2 = 0 998 10
1
10
1 1
10
Log t (day)
1
10
Log t (day)
Figure 3 Results obtained from dataset 3 showing relationships between flow Q(t) and time t in log–log scale. Dots represent the observed data of log (Q) against log t and lines represent model fits using power-law (A), normal exponential (B), inverse time exponential (C), and the combined model (D), respectively. Log is 10-base logarithmic transformation. The parameters estimated from these models are shown in Table 2.
new model using the ordinary linear multivariate regression gives negative values which may be prevented using special type of regression so that only positive regression coefficients can be guaranteed. This will be further investigated.
Concluding remarks It has been discussed in this paper that there are two ways in which non-linearity can be introduced into the flow differential equation. In sections 2 and 3 the implication and implementation of these two types of non-linearity extensions were discussed separately. These two types of non-linearity can be combined into a single equation to derive a solution of combined non-linearity. The differential equation can be expressed as follows: dQ ðtÞ 1 ¼ Q 2n ðtÞ dt nsðtÞ
ð17Þ
where n is the non-linearity index. Using the same approximation with the first three terms of Taylor series expansion applied to the coefficient function s(t) we can get 1 1 n1 Q ðtÞ ¼ c a0 t a1 logðtÞ þ a2 ð18Þ t
In comparison with model (16), model (18) includes an additional parameter (n). This model and its properties will be further discussed in future papers. Two aspects of the river flow recession model can be generalized to incorporate non-linearity: the most commonly used non-linear model (1) involves non-linearity (n) as power transformation applied to the flow, i.e. Qn(t). The other generalization proposed in the paper involves a non-constant monotonically decreasing coefficient function. For the non-linear differential equation with monotonically decreasing coefficient function the solution is a combination of normal exponential, powerlaw and exponential with inverse time forms, which was derived by approximating the coefficient function by the first three terms of a Taylor series expansion. In comparison with other existing multiple parameter flow recession models, the new model involves the combination of components (normal exponential, power-law and exponential with inverse time) which previously had been used as single models for describing flow recession processes. The physical interpretations of the parameters in these earlier models remain approximately valid when the new model is applied. The non-linear model can be expressed as a linear equation by logarithmic transformation and this makes it straightforward to implement the model and validate it graphically for observed river gauging data. The
166
Q. Cheng
A
100
B
100
2
2
Q(t) = 11.38e-0.236t, R = 0.668
Log Q(t)
Log Q(t)
Q(t) = 15.0t-1.027, R = 0.874
10
1
10
1 1
10
1
10
Log t (day)
C
Log t (day)
100
D
100
Q(t) = 1.479e2.736/t, R2 = 0.975
Log Q(t) = -1.198 +0.17t – 2.8log(t) + 2.75/t,
Log Q(t)
Log Q(t)
R2 = 0.995
10
1
10
1 1
10
1
10
Log t (day)
Log t (day)
Figure 4 Results obtained from Dataset 1 showing relationship between flow Q(t) and time t in log–log sacle. Dots represent the observed data of log (Q) against log t and lines represent model fits using power-law (A), normal exponential (B), inverse time exponential (C), and the combined model (D), respectively. Log is 10-base logarithmic transformation. The parameters estimated from these models are shown in Table 2.
Table 2
Results obtained by means of multiple regression
Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset Dataset
1 2 3 4 5 6 7 8 9 10
Constant C
Exponential term a0
Power-law term a1
Exponential term with reversetime a2
Multiple correlation coefficient R2
1.198 0.797 0.356 1.409 2.520 2.210 0.387 1.531 0.915 0.072
0.17 0.00 0.00 0.11 0.20 0.06 0.03 0.06 0.01 0.24
2.80 0.80 0.03 1.86 0.15 2.04 0.21 1.77 0.35 2.62
2.75 0.14 1.37 0.36 0.60 0.23 1.01 0.40 0.34 1.46
0.995 1.000 0.998 0.989 0.995 0.999 0.998 0.988 0.992 0.997
new four-parameter model has been validated with river peak flow data from the ORM area, Canada. The results obtained from 10 datasets from the ORM have clearly demonstrated that the combined model provides better likelihood fitness to the observed flow data.
It should be remembered that the paper deals only with flow recessions in the absence of rainfall. More general results related to continuous flow simulation from rainfall and model validation with simulation will be further investigated.
A combined power-law and exponential model for streamflow recessions
Acknowledgements Thanks are due to Dr. F. P. Agterberg at the Geological Survey of Ottawa, Canada for his constructive comments and improvement of the representation of the manuscript. Four anonymous reviewers are thanked for their constructive comments on reviewing the earlier version of the manuscript. This research has been jointly supported by NSERC Discovery Research ‘‘development of GIS technology for mineral exploration and environmental assessments’’ (ERC-OGP0183993) and Chinese National Foundation of Science Projects ‘‘Distinguished Young Researcher Grant’’ (40525009) and Strategic Research project (40638041).
References Beven, K.J., 2000. Rainfall–runoff Modeling: The Primer. John Wiley and Sons, Chichester, UK, p. 360. Boussinesq, J., 1904. Recherches theoretique sur l’ecoulement des nappes d’eau infiltrees dans le sol et sur le debit des sources. J. Math. Pure Appl. 10 (5th series), 5–78. Brutsaert, W., Nieber, J.L., 1977. Regionalized drought flow hydrographs from a mature glaciated plateau. Water Resour. Res. 13 (3), 637–643. Cheng, Q., 2004. Weights of evidence modeling of flowing wells in the Greater Toronto Area, Canada. Nat. Resour. Res. 13 (2), 77– 86. Cheng, Q., Russell, H., Sharpe, D., Kenny, F., Pin, Q., 2001. GISbased statistical and fractal/multifractal analysis of surface stream patterns in the Oak Ridges Moraine. Comput. Geosci. 27 (5), 513–526. Cheng, Q., Zhang, G., Lu, C., Ko, C., 2004. GIS spatial-temporal modeling of water systems in the Greater Toronto Area, Canada, Earth Sciences. J. China Univ. Geosci. 15 (3), 275–282. Cheng, Q., Ko, C., Yuan, Y., Ge, Y., Zhang, S., 2006. GIS Modeling for predicting river runoff volume in ungauged drainages in the Greater Toronto Area, Canada. Comput. Geosci. 32 (8), 1108– 1119. Coutagne, A., 1948. Etude generale des variations de debits en function des facteurs qui les conditionnent, 2eme partie, les variations de debit en periode non influence par precipitations’, La Houille Blanche, Sept. Oct., pp. 416–436. Croke, B.F.W., 2006. A technique for deriving an average event unit hydrograph from streamflow – only data for ephemeral quickflow-dominant catchments. Adv. Water Resour. 29, 493–502. Hall, F.R., 1968. Base flow recession – a review. Water Resour. Res. 4 (5), 973–983. Horton, R.E., 1933. The role of infiltration in the hydrologic cycle. Trans. Am. Geophys. Union 14, 446–460. HYDAT CD-ROM User’s Manual, 1996. Surface Water and Sediment Data, Atmospheric Environment Program, Version 96 – 1.04 User’s Manual, Environment Canada, 95p. Ishihara, T., Takagi, F., 1965. A study of the variation of low flow. Disaster Prev. Res. Inst. Kyoto Univ. Bull. 15 (2), 75–98. Jakeman, A.J., Hornberger, G.M., 1993. How much complexity is warranted in a rainfall–runoff model? Water Resour. Res. 29 (8), 2637–2649.
167
Jakeman, A.J., Littlewood, I.G., Whitehead, P.G., 1990. Computation of the instantaneous unit hydrograph and identifiable component flows with application to two small upland catchments. J. Hydrol. 117, 275–300. Jakeman, A.J., Post, D.A., Beck, M.B., 1994. From data and theory to environmental model: the case of rainfall–runoff. Environmetrics 5, 297–314. Kenny, F.M., 1997. A chromostereo enhanced digital elevation model of the Oak Ridges Moraine Area, Southern Ontario and Lake Ontario Bathymetry, Geological Survey of Canada, Open file 3423, Scale 1:200,000. Ko, C., Cheng, Q., 2004. GIS spatial modeling of river flows and precipitation in the ORM, Ontario. Comput. Geosci. 30 (4), 379– 389. Kokkonen, T.S., Jakeman, A.J., Young, P.C., Koivusalo, H.J., 2003. Predicting daily flows in ungauged catchments: model regionalization from catchment descriptors at the Coweeta Hydrologic Laboratory, North Carolina. Hydrol. Process. 17, 2219–2238. Littlewood, I.G., 2002. Improved unit hydrograph characterisation of the daily flow regime (including low flows) for the River Teifi, Wales: towards better rainfall–streamflow models for regionalisation. Hydrol. Earth Syst. Sci. 6 (5), 899–911. Littlewood, I.G., Down, K., Parker, J.R., Post, D.A., 1997. IHACRES v. 1.0 User Guide. Center for Ecology and Hydrology, Wallingford, UK & Integrated Catchment Assessment and Management Center, Australian National University, Canberra, 94pp. Maillet, E., 1905. Essai d’hydraulique souterraine et fluviale, Libraire Sci., A, Herman, Paris (Cited by Hall, 1968). Mitchell, W.D., 1972. Model hydrographs. US Geol. Surv. Water Supply Pap. 2005. Moore, R.D., 1997. Storage-outflow modeling of streamflow recession, with application to a shallow-soil forested catchment. J. Hydrol. 198, 260–270. Napiorkowski, J.J., Strupczewski, W.G., 1979. The analytical determination of the Kernels of the Volterra Series describing the cascade of nonlinear reservoirs. J. Hydrol. Sci. 6 (3-4), 121– 142. Sefton, C.E.M., Howarth, S.M., 1998. Relationships between dynamic response characteristics and physical descriptors of catchments in England and Wales. J. Hydrol. 211, 1–16. Tallaksen, L.M., 1995. A review of baseflow recession analysis. J. Hydrol. 165, 349–370. US Department of Agriculture, Soil Conservation Service, 1972. National Engineering Handbook, Section 4, Hydrology, US Government of Printing Office, Washington, DC, 544p. Wagener, T., Wheater, H., Gupta, H.V., 2004. Rainfall–runoff Modeling in Gauged and Ungauged Catchments. Imperial College Press, London, 306p. Werner, P.W., Sundquist, K.J., 1951. On the groundwater recession curve for large watersheds, IAHS General Assembly, Brussels, IAHS Publ., 33, pp. 202–212. Wittenberg, H., 1999. Baseflow recession and recharge as nonlinear storage processes. Hydrol. Process. 13, 715–726. Yuan, Y., Cheng, Q., 2005. Predicting streamflow in ungauged basins with GIS and rainfall–runoff models: a case study of model regionalization in the Oak Ridges Moraine area, Southern Ontario, Canada. In: Qiuming, Cheng, Graeme, Bonham-Carter (Eds.), Proc. of the Int. Ass. for Math. Geol. Annual Conference, 1, pp. 428–433.