A dynamic rating curve approach to indirect discharge measurement

The operational measurement of discharge in medium and large rivers is mostly based on indirect approaches by converting water stages into discharge on the basis of steady-flow rating curves. Unfortunately, under unsteady flow conditions, this approach does not guarantee accurate estimation of the discharge due, on the one hand, to the underlying steady state assumptions and, on the other hand, to the required extrapolation of the rating curve beyond the range of actual measurements used for its derivation. Historically, several formulae were proposed to correct the steady-state discharge value and to approximate the unsteady-flow stage-discharge relationship. In the majority of these methods, the correction is made on the basis of water level measurements taken at a single cross section where a steady state rating curve is available, while other methods explicitly account for the water surface slope using stage measurements in two reference sections. However, most of the formulae available in literature are either over-simplified or based on approximations that prevent their generalisation. Moreover they have been rarely tested on cases where their use becomes essential, namely under unsteady-flow conditions characterised by wide loop rating curves. In the present work, an original approach, based on simultaneous stage measurements at two adjacent cross sections, is introduced and compared to the approaches described in the literature. The most relevant feature is that the proposed procedure allows for the application of the full dynamic flow equations without restrictive hypotheses. The comparison has been carried out on channels with constant or spatially variable geometry under a wide range of flood wave and river Correspondence to: F. Dottori (francesco.dottori@unibo.it) bed slope conditions. The results clearly show the improvement in the discharge estimation and the reduction of estimation errors obtainable using the proposed approach.


Introduction
Discharge measurement is an issue of major importance for the evaluation of water balance at catchment scale, for the design of water-control and conveyance structures, for rainfallrunoff and flood routing model calibration and validation.
Although several direct measurement approaches exist, only indirect approaches tend to be used operationally in medium and large rivers.Usually, discharge estimates are based on a one-to-one stage-discharge relationship, or steady-flow rating curve, which is derived on the basis of a number of simultaneous stage and discharge measurements.A measure of stage is then directly converted into discharge by means of the developed rating curve.
Such an approach can be considered adequate for all rivers under steady-flow conditions, and also under unsteady-flow conditions, when flood waves show a marked kinematic behaviour, which generally corresponds to rivers with steep bed slopes (>10 −3 ).In all other cases the variable energy slope, associated with the dynamic inertia and pressure forces relevant to the unsteady flow discharge, lead to the formation of a hysteretic rating curve also known as the loop-rating curve (Jones, 1916;Chow, 1959, p. 605;Henderson, 1966, p. 392;Fread, 1975).This implies that the steady-flow rating curve is no longer sufficient and adequate to describe the real stagedischarge relationship.Recently, with a numerical study on the River Po, Di Baldassarre and Montanari (2008) showed that the use of the steady-flow rating curve may lead to major F. Dottori et al.: Dynamic rating curve approach to indirect discharge measurement errors in discharge estimation when significant flood waves occur, which may be greater than 15%, and that another significant error is produced by the extrapolation of the rating curve beyond the range of measurements used for its derivation.
In addition to the water balance error induced by the hysteretic effect and the extrapolation, another error occurs that may strongly affect the calibration and the verification of hydrological models: if calibration is made using discharge values derived from a steady-flow rating curve, then the estimated time of peak discharge will be wrong, because, under unsteady flow conditions, the peak discharge occurs before the maximum water stage, and this delay can be significant (several hours) in very mild river slope conditions.Schmidt and Garcia (2003) described different methods historically used to overcome this problem; these methods mainly consist of empirical adjustments of the rating curve, derived from experimental data, while, less frequently, especially in river reaches affected by backwater effects, estimations are adjusted using a reference value of water surface slope, computed as the "fall", or difference in water level between the section of interest and a second reference section, where stage is known (Herschy, 1995).
Aside from more or less empirical approaches, several formulae based on full or simplified dynamic flow equations have been developed to account for the observed hysteresis in stage-discharge relationship.Unfortunately, very few comparisons of the different methods can be found in the literature (Perumal and Moramarco, 2005).Some of these methods explicitly account for the water surface slope using stage measurements in two reference sections, while a larger number of authors, for practical reasons, convert the surface slope into a time derivative, which is then estimated using successive water stage measurements at the same cross section where a rating curve is available.However, most formulations were, on one hand, obtained under restrictive hypotheses on flow and river bed geometry, thus reducing the possibilities for operational applications (Schmidt and Yen, 2002;Perumal et al., 2004), and, on the other hand, rarely tested to verify their validity range.
In this paper, after an excursus on the historically proposed approaches, an alternative methodology is introduced, which explicitly accounts for the longitudinal variation of the water surface slope, through the use of couples of simultaneous water stage measurements at two adjacent cross sections.Such procedure, which also requires the geometrical description of the two cross sections, allows for the application of the full dynamic flow equations without restrictive hypotheses.
The proposed methodology is fully described and compared, on the basis of several test bed experiments, to the wide variety of existing approaches that can be found in the literature.

Data and methods
All the methods presented in this paper derive from the 1-D shallow water momentum equation, by disregarding one or more terms: where ∂z ∂x = ∂z 0 ∂x + ∂y ∂x is the water surface slope, composed of the channel bottom slope and the pressure force term; K 2 is the friction slope; and t, the time coordinate [s]; x, the longitudinal distance along the reach [m]; z, the water surface elevation above a horizontal datum [m]; z 0 , the river bed elevation above a horizontal datum [m]; y, the water depth [m]; A, the wetted area [m 2 ]; Q, the river discharge [m 3 s −1 ]; g, the gravity acceleration [m s −2 ]; and β, the Boussinesq momentum coefficient.
Please note that in the present paper, following the findings of Gasiorowski and Szymkiewicz (2007), Eq. ( 1), is written in conservative form, by keeping under the derivative sign the non-linear terms appearing in the convective and the local acceleration terms.
As previously stated, all the formulae which try to provide improved discharge estimates from water level measurements, basically derive from Eq. (1).A comprehensive bibliography of existing methods for unsteady flow discharge estimation can be found in Fenton andKeller (2001), andPetersen-Øverleir (2006).As discussed in the sequel, several other authors have also proposed original formulae or modifications to previous formulae (see Sect. 2.1 and 2.2), while others have carried out comparisons or evaluations of existing formulae, using both numerical simulation or measured data in natural rivers (see Sect. 2.3).
In this paper, given the large number of available formulae, it was preferred to present them in chronological order.Moreover, for reasons of space, the analytical derivation of each formula from the original dynamic flow equation has been omitted, but can be easily found in the referenced works.Finally, for the sake of clarity some of the original symbols were converted in order to preserve consistency all throughout this paper, and the complete list of symbols used is provided in Appendix A.

Jones formula
Among the formulae existing in the literature, the Jones formula is, without any doubt, the most well-known.
Hydrol.Earth Syst.Sci., 13, 847-863, 2009 www.hydrol-earth-syst-sci.net/13/847/2009/ Jones (1915) used the Chezy equation of friction slope and a geometric analysis to estimate the water-surface slope based on the surface velocity and the rate of change of the stage at the gauge.Consequently, unsteady-flow discharge can be computed as: in which, Q r is a "reference" discharge for the given stage, S r is a "reference" water-surface slope and U s is the surface velocity, which Jones defined as the mean velocity U divided by 0.90 for large streams and divided by 0.85 for smaller streams.According to Jones, Eq. ( 2) needs to be calibrated by measuring S r under different flow conditions and in the presence of a constant discharge.Therefore, Q r may be set equal to the steady flow discharge Q s ,and S r is the water surface slope under steady flow condition S s .
Since its publication, the Jones formula has been the subject of many research works, either as the starting point for obtaining more accurate equations, or for establishing a general applicability criterion; a number of these works are herein reviewed.Chow (1959, p. 532) found that the Jones formula can be derived from the momentum Eq. ( 1) by neglecting the convective and local acceleration terms and assuming that the flood wave moves downstream with a constant velocity and without changing its profile (the so-called "uniform-progressive" flow); if the steady flow condition is also the uniform flow, in which S s =S 0 , Eq. ( 2) may be rewritten as: where the kinematic wave celerity c can be approximated according to the definition given by Henderson (1966, p. 367): Equation ( 3) is the standard form of Jones formula and has been used in almost all of the successive works.Perumal and Ranga Raju (1999) pointed out that if Eq. ( 3) may be regarded as the approximate convection-diffusion (ACD) equation, which has the same form as that of the kinematic wave equation, it can therefore be used to describe the flood wave characterized by a narrow loop.(1966, p. 393) proposed a modification of Eq. (3) through the introduction of a term, based on a parabolic approximation of flood wave, which accounts for wave subsidence:

Henderson
where r=S 0 / (∂y/∂x) is the ratio of the channel bottom slope to the entering wave slope.According to Henderson (1966), the term r can be approximated from the characteristics of a typical flood event for the concerned reach; r is therefore given by the ratio of wave height to its half-length, the latter computed from the product of average wave celerity c and the time to peak stage (the typical wave is assumed to be kinematic).

Di Silvio formula
Di Silvio (1969) used a triangular approximation of the flood wave and the hypothesis of narrow loop in the rating curve to obtain a formula for discharge estimation; for the rising limb the relation is: while in the receding limb it is: In Eqs. ( 6) and (7) ϕ= 2 m p 1 − R ∂C ∂A ; Q b and Q p are the base and the peak discharge of the flood wave; T rise and T receding are the duration of the rising and receding limb respectively; A p and R p are the area and hydraulic radius values corresponding to peak discharge; A is the mean between area values corresponding to base and peak discharge; m is the exponent of the hydraulic radius in the friction law used (for instance, when using Chézy expression, m=1/2); p is the exponent of the wetted area in the friction law used (for instance, when using Manning expression, p=2).

Fread formula
Fread (1975), using the same approximation for the space derivative introduced by Henderson (1966), derived a hysteretic rating curve model from the full one-dimensional unsteady flow equations.The proposed formulation allows the computation of either the discharge Q or the water stage y as a function of the time derivative of the other variable: ∂B ∂y and r is the ratio of the channel bottom slope to the entering wave slope, which can be computed F. Dottori et al.: Dynamic rating curve approach to indirect discharge measurement using the following expression, similar to that proposed by Henderson for Eq. ( 5): where Q b , Q p and Ā are defined as for Eqs. ( 6) and ( 7); T rise is the duration of rising limb expressed in days; h b and h p are the water stages corresponding to base and peak discharge.
The underlying hypotheses of Eq. ( 8) are correct in kinematic, or quasi kinematic flow conditions, in wide channels with approximately constant width (Fread, 1975).Also note that Eq. ( 8) is implicit and therefore must be solved via iterative methods.Marchi (1976) proposed the following simplified relationship for estimating the unsteady-flow rating curve:

Marchi formula
The derivation of Eq. ( 10) can also be found in Perumal and Moramarco (2005).Equation ( 10) can be obtained through a kinematic approximation of momentum equation, which leads to the following differential equation: where η is practically a constant.

Faye and Cherry formula
The method developed by Faye and Cherry (1980) combines both momentum and continuity equations to obtain a single expression in which the pressure gradient∂y/∂x is substituted using the kinematic wave approximation, under the hypothesis of stable wave profile during the downstream routing.Evaluation of such expression by finite-difference approximation leads to the following quadratic equation: In Eq. ( 12), ∂U and ∂t are evaluated by backward difference aproximations and ∂y is evaluated by a central-difference approximation.The subscripts indicate the time step at which variables have to be computed.
2.1.7Lamberti and Pilati formulae Lamberti and Pilati (1990) developed two equations designed to compute the difference between steady and unsteady flow rating curves: where ; a and b are first and second order incremental ratios defined as a= ; S z is the water surface slope ∂z/∂x, which can be approximated by bed slope S 0 .Both formulae can be solved without iterations, using the terms T 1 and T 2 computed in the previous time step.
Equations ( 13) and ( 14) can be applied in kinematic or quasi kinematic conditions, that is, in presence of narrow loops of the rating curve, with a maximum difference between unsteady and steady flow rating curve of about 10%.The authors provided a quantitative criterion to establish the ratio Q/Q 0 from channel and wave characteristics: where T 1 is the characteristic channel time, as defined for Eqs. ( 13) and ( 14), and T rise is the duration of the wave rising limb.

Fenton formulae
Fenton (1999) worked on the original unsimplified shallow water equations, reducing them to a single expression in which space derivatives are substituted by time derivatives using a polynomial Taylor series.The complete procedure can be found in Fenton and Keller (2001).Fenton proposed two formulations: the first one may be regarded as an extension of Jones formula, which includes a diffusive term: where D= Q 0 2BS 0 is the diffusion coefficient.The second proposed expression includes a third order time derivative of the water stage: where the term G is given by: 2.1.9Perumal formulae Perumal andRanga Raju (1999), andPerumal et al. (2004) refined the time derivative of the Jones formula by incorporating expressions for the inertial forces of the onedimensional momentum equation.They obtained two different expressions for unsteady-flow rating curve: the first one has the following form: where m is the exponent of the hydraulic radius in the friction law used (as in Sect.2.1.3Eqs.6 and 7); the second equation is: Perumal et al. ( 2004) also identified a criterion to establish the suitability of Eqs. ( 3), ( 19) and ( 20), as a function of bed and wave slopes; according to the authors, the estimation given by these methods may be considered good if the following condition holds:

Chow formula
Many authors presented a general discussion over the magnitude of the different terms composing Eq. ( 1) (see Henderson, 1966, p. 364;Todini and Bossi, 1986;Lamberti and Pilati, 1996;Schmidt and Yen, 2003).In most rivers, during a flood event the local and the convective acceleration terms in Eq. ( 1) can be neglected because their values range from one tenth to one hundredth of the other terms appearing in the equation.By neglecting the convective and the local acceleration terms, a parabolic approximation of the full de Saint Venant equations can be obtained, which leads to the expression: Using Eq. ( 22) under both hypotheses of prismatic channel and uniformly progressive wave (that is, (1959, p. 532) obtained an expression for unsteady flow computation: As previously stated in Sect.2.1.1,when the flood wave behaves as a kinematic wave the longitudinal gradient of water stage can be directly related to the time derivative of the stage, by means of the kinematic celerity: from which the Eq. ( 3), that is the Jones formula, can be derived.Therefore, the Jones formula can be regarded as an approximation of the parabolic assumptions used by Chow, and is valid when approaching the kinematic conditions expressed by Eq. ( 24).

Fenton and Keller formula
Chow formula can only be correctly applied for prismatic channels where for which the uniform flow is meaningful.To overcome this limitation, Fenton and Keller (2001) pointed out that Eq. ( 22) can be directly used to estimate discharge, without the need of introducing kinematic conditions expressed by Eq. ( 24), by writing it in the following form: Therefore, it is worth noting that Eq. ( 25) may be seen as an extension to a more general version of Chow formula (Eq.23), which overcomes the limitations of prismatic channels and of quasi kinematic flow conditions.

The stage-fall-discharge method
This method is described in detail by Herschy (1995), and is a way to correct the steady flow rating curve in channels influenced by backwater conditions; the method requires direct measurements of the "fall", that is, the difference in water surface level measured between two sections.The expression is: where S z is the water surface slope and S r is a reference fall.According to Herschy, S r may be assumed either as constant where backwater effects are always present, or variable if this is not the case, thus corresponding to steady flow water surface slope S s .The author also describes a procedure for the calibration in both the cases, using measured values of stage, fall and discharge.Although Herschy does not provide a theoretical background for Eq. ( 26), Fenton and Keller (2001) point out that the equation derives from the Chezy law, and the same expression may also be found in Jones (1915) as the starting point for the development of Eq. (3).
Actually, the stage-fall-discharge method may be directly derived from Eq. ( 22) by dividing the unsteady flow discharge equation Q=K √ S z , by the corresponding steady F. Dottori et al.: Dynamic rating curve approach to indirect discharge measurement flow discharge equation Q s =K √ S s computed at the same water depth, which implies the same hydraulic conveyance K.
A more comprehensive expression would involve the convective acceleration for friction slope computation.Therefore the stage-fall-discharge method should be rewritten as in Schmidt and Yen (2002): where the subscript "s" indicates that computation must be referred to the steady flow rating curve, taken as the reference.

The proposed DyRaC formula
Given that parabolic approximation may be not suitable in non prismatic channels, where the effect of longitudinal variation of cross sections may mean that the convective acceleration term in Eq. ( 1) becomes of the same magnitude as the other terms (Schmidt and Yen, 2002), Eq. ( 1) can be used directly to derive the dynamic rating curve.Alternatively, when the parabolic approximation is sufficiently accurate, by neglecting the local acceleration term, it is possible to derive the dynamic rating curve from the following equation: as given in Aricò et al. (2008), Dottori et al. (2008).
In any case, the proposed approach basically neglects the continuity of mass equation between the two cross sections on which Eq. (1) (or Eq. 28) is applied, by assuming (1) that no significant discharge enters (or leaves) the reach between the two adjacent sections, and (2) the two cross sections are close enough to accept the hypothesis that ∂Q/∂x ∼ =0.On these grounds it is then possible to discretise Eq. ( 1) between the upstream and the downstream cross sections, to obtain: or, by neglecting the local acceleration as in Eq. ( 28), to obtain: Please note that Eq. ( 30) is none other than the standard equation used for the estimation of the water surface profile under steady (but non uniform) flow assumptions.In Eqs. ( 29) and ( 30), the upstream and downstream conveyance values K u andK d can be computed assuming a constant energy slope along the cross section (Chow, 1959, p. 139).Each cross section is divided in k subsections, each with conveyance K j , and the total conveyance can be expressed as a function of the corresponding subsection conveyances, as: while the Boussinesq momentum coefficient can be estimated as: where v(a) represents the velocity at an infinitesimal element of the wetted area of the cross section.
Please note that the distance between the two adjacent sections must be sufficiently small to allow for the constant flow rate assumption to be realistic, but at the same time it must be sufficiently large to allow the difference in water stage to be greater than the measurement instrument sensitivity and the water elevation fluctuations.
Equation ( 30) can be solved explicitly with respect to Q, to give: Once the water levels in the upstream and downstream sections are measured, for a given roughness all the terms of Eq. ( 33) are known as a function of the water stage.Therefore the equation can be used, similarly to a standard rating curve, to dynamically estimate the discharge values as a function of the water level as well as of the water surface slope, which continuously varies in time.This differs from the use of the classic steady-flow rating curve, which only depends on the water depth by implicitly assuming an average, but constant in time, water surface slope.Therefore, due to its dynamic nature, in this paper the new approach will be called the Dynamic Rating Curve (DyRaC).Whenever needed, namely when the local acceleration term is not negligible, the DyRaC expression can be expanded by re-deriving it from Eq. ( 29), to give: Hydrol.Earth Syst.Sci., 13, 847-863, 2009 www.hydrol-earth-syst-sci.net/13/847/2009/ where the time derivative ∂ (Q/A)/∂t, can be approximated using the incremental ratio Ū / t, where t is the sampling time step and Ū is the average velocity within the reach, which can be estimated as Ū =2Q/ (A u + A d ), which leads to: where Ūt− t =2Q t− t / (A u + A d ) t− t is the average velocity computed at the previous time interval.
As opposed to Eq. ( 33), which is explicit in terms of discharge, Eq. ( 35) must be solved iteratively.This can be easily done using a simple Newton-Raphson approach, which converges to the required accuracy in a very limited number of iterations (∼5-6).
Nonetheless, it will be shown that the results obtained using Eq. ( 33) are already adequate to accurately estimate the discharge in natural rivers.

Design and preparation of numerical experiments
As described in Sect.2.1 and 2.2, several methods for unsteady flow discharge have been developed and are available in the literature; however, the literature does not offer useful criteria for a comprehensive evaluation of methods, nor for identifying the most appropriate ones, for the different application conditions.
What mainly emerges from the literature is a lack of publications dealing with a comprehensive comparison of the different methods; Perumal and Moramarco (2005) addressed this issue but their analysis was limited to few methods, either well-known or developed by the authors themselves.
In terms of application conditions and ranges, it appears that most methods have been designed to provide unsteady flow discharge estimation in kinematic or quasi kinematic conditions; in such conditions, due to the limited amplitude of the unsteady flow loop, the proposed correction formulae produce limited improvements with respect to what is obtained using the steady flow rating curve.Other methods also use the restrictive hypothesis of constant width channel, which limits their operational use in natural rivers.
Moreover, the issue of practical application of these methods appears to be seldom addressed in the literature.Few works present an extensive operational use of discharge estimation formulae in natural rivers.Barbetta et al. (2002) and Franchini and Ravagnani (2007) carried out formulae applications in quasi kinematic flow conditions, resulting, as mentioned earlier, in limited improvements.In a few cases, formulas have been applied to flood waves characterized by wide loop rating curves (Fread, 1975;Faye and Cherry, 1980;Petersen-Øverleir, 2006) while other authors address the issue of formulae evaluation using numerical experiments (Lamberti and Pilati, 1990;Fenton, 1999;Perumal et al., 2004).
Therefore, one of the objectives of the present paper is to compare the existing methods and to test their reliability under different application conditions.Thus, a number of numerical experiments were set up, to simulate a wide range of flow conditions over channels with different bed slope and geometry.These experiments, which attempt to reproduce hydraulic conditions observed in natural rivers, are summarised in Table 1.
The values of bed slope used in the experiments vary from 10 −3 (steep slope) to 2.5×10 −5 (very mild slope), including the intermediate values of 5×10 −4 , 2×10 −4 , 10 −4 and 5×10 −5 ; three types of wave were used in the simulations: a fast wave with a rising time of 24 h and a peak discharge of 900 m 3 s −1 , a medium wave with a rising time of 72 h and a peak discharge of 900 m 3 s −1 and a slow wave with a rising time of 168 h and a peak discharge of 10 000 m 3 s −1 .The  choice of the bed slope values and the flood wave characteristics was made bearing in mind the results of numerical experiments carried out in previous works (Lamberti and Pilati, 1990;Perumal et al., 2004), in order to analyse not only typically kinematic or quasi kinematic flow conditions, but also to explore the range between kinematic and parabolic flow conditions.
In addition, the values of peak discharge, flood wave duration and channel geometry were chosen as a function of bed slope values, in order to recreate flow conditions close to those which usually take place in natural rivers; for example, the channels with a bed slope of 5×10 −5 and 2.5×10 −5 have a section width of 400 m, much larger than the other channels with steeper bed slopes.The geometry of channels used in the numerical experiments is described in more detail in the sequel: cases from 1 to 8 relate to a channel with rectangular cross sections and constant width; cases 9, 10 and 11 were introduced to assess the different expressions under variability of cross sections, and in particular case 9 is characterised by a cross section changing from rectangular to trapezoidal, while cases 10 and 11 relate to a channel with irregular cross sections (Fig. 1), as one might expect from natural water courses.
The flood waves were generated in all the cases using the following expression: where Q b is the base flow discharge (equal to 100 m 3 s −1 in all cases), T p the time to peak flow, Q p the peak discharge and γ a coefficient assumed to be equal to 16. Please note that the term in square parenthesis in Eq. ( 36) is raised to power 16.This produces waves that, although the time to peak is T p , will grow infinitesimally for many hours, and will appear as waves with a raising limb of a time duration significantly shorter than T p , as can be noticed from the resulting figures.
All the simulations were made using two well-known 1-D hydraulic models, Hec-Ras (HEC, 2001) and Mike11 (DHI, 2003), in order to assess the reliability of the results.The results of the simulations using the two models proved to be very similar in all cases both in terms discharge and stage values.These results were thus taken as the "true" values in order to assess the validity of the different formulae.

Formulae assessment
Before comparing the different approaches, the suitability of each method was assessed according to the criteria established by the authors or by other researchers in successive works.
The suitability of the Jones formula (Eq. 3) and of derived formulae presented by Perumal et al. (2004, Eqs. 9 and 20) were evaluated using the criterion expressed by Eq. ( 21), computed for all the simulation time steps; applications show that these formulae should provide acceptable discharge values in cases 1 (fast wave over steep river bed slope), 2 (fast wave over medium river bed slope) and 3 (medium wave over medium-mild river bed slope); in cases 4 (fast wave over medium-mild river bed slope) and 5 (medium wave over mild river bed slope) the values obtained using criterion of Eq. ( 21) are occasionally greater than the threshold value, which means that estimation could be locally inexact, while in the remaining cases the results from formulae are expected to be unreliable.According to the analysis made by Perumal and Moramarco (2005), the same results may also be considered valid for the Marchi and Fenton formulae (Eqs. 10,16 and 17).
The Di Silvio formula (Eqs.6 and 7) was not evaluated, since it requires the knowledge of flood wave peak and duration in order to be applied, that is, it can not be operationally used in real time, although it could be used after the flood peak has passed.
An analogous procedure was applied for the two equations developed by Lamberti and Pilati (1990, Eqs. 13 and 14); according to the criterion given by Eq. ( 15), these formulae should correctly estimate the discharge in cases 1, 2 and 3, while in remaining cases the results would be incorrect.Henderson (1966) did not provide a quantitative criterion for his formula, neither did Fread (1975) and Faye and Cherry (1980); however, these three equations (Eqs.5, 8 and 12, respectively), share the hypotheses of kinematic wave approximation and stable wave profile during the downstream translation.Therefore when such hypotheses are no longer valid significant errors are expected: namely in cases 4 (fast wave over medium-mild river bed slope), 6 (fast wave over mild river bed slope) 7 and 8 (slow wave over very mild river bed slope).
Given that Chow equation Eq. ( 23) and Fenton and Keller equation Eq. ( 25) are identical in prismatic channels, only Chow equation is referenced to in Figs. 3 to 8 showing the comparison on cases 1 to 8, while Fenton and Keller equation is referenced in Fig. 9, which deals with a comparison on non-prismatic channels.
Finally, the DyRaC formulae (Eqs.33 and 35) are theoretically reliable under all flow conditions, in particular Eq. ( 35) is needed when the influence of the local acceleration term in Hydrol.Earth Syst.Sci., 13,[847][848][849][850][851][852][853][854][855][856][857][858][859][860][861][862][863]2009 www.hydrol-earth-syst-sci.net/13/847/2009/ Eq. ( 33) is not negligible, since this term may become significant in channels and rivers with very mild slopes subject to fast rising flood waves (hyperbolic flood wave conditions).Several simulations were carried out in order to assess the relevance of the local acceleration term in the numerical experiments used to evaluate the effectiveness of the different equations.The results are presented in Fig. 2a and b, in terms of R a , the ratio of the local acceleration term and the hydraulic head slope.The figures relate to cases 6 and 8, where R a reaches its maximum values.As can be seen from the figures, the local acceleration term is always negligible, since R a , which is plotted versus the hydraulic head slope reaches at most 1% of the latter.Therefore, due to the very small magnitude of the local acceleration term in all the reported experiments, which were chosen close to natural flood wave conditions in rivers, Eq. ( 33) was always used instead of Eq. ( 35) since it provides the same results, without requiring an iterative solution.Moreover, it should be noted that the waves simulated in cases 6 and 8 are significantly faster than flood waves generally taking place in natural rivers with similar bed slopes.For example, the bed slope of the final reach of the River Po in Italy is around 510 −5 , while at the same time, the rising time of the flood waves is generally longer than one week (168 h), with a rate of change in stage of few cm h −1 .Therefore, although Eq. ( 35) can always be used when the inertial term becomes significant, it must be stressed that Eq. ( 33) can probably be applied, for practical operational purposes, on all types of natural rivers and under all flood conditions.

Operational estimation of discharge in natural rivers
Another topic of major relevance is the reliability of the reviewed methods under operational conditions.Generally, the formulae presented in this paper were tested using high precision data from numerical or laboratory experiments, by assuming perfect water stage measurements, whereas operationally water stage measurements in natural rivers are generally affected by measurement errors (typically around ±1 cm) in terms of instrument precision, while local oscillations of the water surface can add additional uncertainty; as a consequence it is not possible to arrive at a correct estimate of the real discharge using single instantaneous measurements.An alternative methodology to provide reliable estimates can be applied by installing gauge stations with sensors capable of carrying out a number of discharge estimates in a limited amount of time (a few minutes), during which discharge can be considered as constant.This permits an iterative computation of the expected value of discharge µ(Q), using the following equation: Where i is the number of measurements; Q i is the i-th computed discharge value; µ i (Q) and µ i−1 (Q) are the mean val- ues computed using (i) and (i-1) measurements.The standard deviation of the computed values may also be estimated as: Where µ i (Q 2 ) is the mean of square values of Q, estimated using the following recursive equation: The accuracy of the estimated mean value of discharge is given by the standard deviation of the mean, defined as: As can be seen from Eq. ( 40), the uncertainty of the estimation of Q improves at each new measure, so that the procedure can be iterated until the error of estimation falls below a required precision.
The effectiveness of the proposed methodology needs then to be tested by showing the actual number of iterations required to reach an acceptable level of precision, which, for practical purposes, must be limited.In the present paper, the methodology has been assessed by applying the following procedure: the reference values of water stage (computed by the hydraulic model as stated in Sect.2.3) were perturbed by adding a random error, computed using a Gaussian distribution with zero mean and a standard deviation of √ 5 cm, roughly comparable with an error deriving from the accuracy of water stage sensors (±1 cm) and from the water surface oscillations (±2 cm).For each time step a set of perturbed stage values was produced to simulate a series of continuous sensor measurements.Then the procedure, starting from a minimum number of 10 and 20 couples of simultaneous stage measurements was iterated until the standard deviation σ i (µ i (Q)) reached a value smaller than 5% with respect to the mean µ i (Q).This was done by defining the following indicator I σ/µ , which interrupted the calculation when <0.05: A similar approach can also be used to estimate a roughness-depth relationship given a series of discharge measures and, for each discharge measured value, several simultaneous couples of water stage measurements; it is worth noting that a number of investigations pointed out the need of relating roughness value with water depth by calibrating a continuous relationship (Simons and Richardson, 1961;Fread, 1975); such procedures have been used with very good results in a previous application of DyRaC methodology (Dottori et al., 2008).
Once the roughness-depth relationship is established, the conveyance can be easily computed for each water depth values, given the knowledge of cross section geometry from which area, width and wetted perimeter can be derived.

Analysis of results
The estimated discharge values produced by the different formulae were evaluated by comparing the mean error and the error variance with respect to the discharge "true" values, namely the ones computed by using the hydraulic model (Sect.2.3), taken as "true".In a first set of experiments (Sect.3.1 and 3.2) the water stage measurements were considered as "perfect", namely not affected by measurement errors.The effect of measurement errors was then assessed and is discussed in Sect.3.3.

Comparison on channels with prismatic constant section
Figures 3 and 4 show the mean error and the error variance of the succession in time of the discharge estimates produced by the alternative formulae for cases 3, 4, 5 and 6; the values obtained for the other cases were not represented to allow a clearer representation of results since the values were either very low (for cases 1 and 2) or very high (for cases 7 and 8) with respect to those presented in the two graphs.In addition, some of the formulae were omitted because of the strong similarities existing among them: for instance Eq. ( 19) parameters almost coincide with those of the Jones formula (Eq.3), which was also found in a previous analysis work by Perumal et al. (2004).
As expected, in all cases from 1 to 8 the Chow and Fenton and Keller formulae (Eqs.23 and 25, respectively) provide identical results, since the steady flow condition coincides with the uniform flow, and the results from stage-falldischarge method (Eq.26) are always equal to those of the Eqs.( 23) and ( 25).
Finally, in all cases from 1 to 8 Eqs. ( 23) and ( 25) and the DyRaC formula (Eq.33) give the same results, which is not surprising given the use of prismatic cross sections.
As expected, the ability of the different equations to estimate discharge strongly depends on the channel and flood wave characteristics.
In cases 1 and 2 (fast wave over steep and medium river bed slope), the mean error is always below 2 m 3 s −1 for all the formulae and the percentage errors at peak are less than 1.2%, which means that the "true" values are all very well reproduced.However this is also true for the values given by the steady-flow rating curve: the discharge-level hydrograph (Fig. 5, left) and the comparison between steady and unsteady flow rating curves (Fig. 5, right) for case 2 shows the absence of a significant loop, which implies that flow conditions can be considered quasi kinematic.
In cases 3 (medium wave over medium-mild river bed slope), 4 (fast wave over medium-mild river bed slope), 5 (medium wave over mild river bed slope) and 7 (slow wave over mild river bed slope), the degree of accuracy is more variable since incoming waves become progressively steeper with respect to bed bottom slope; nonetheless, it can be seen that the DyRaC formula (Eq.33) maintains a very low error rate, and that Perumal 2 (Eq.20), Henderson (Eq. 5) and Fread (Eq. 8) formulae perform slightly better than other ones (see the hydrographs of case 4 presented in Fig. 6).
Hydrol.Earth Syst.Sci., 13,[847][848][849][850][851][852][853][854][855][856][857][858][859][860][861][862][863]2009 www.hydrol-earth-syst-sci.net/13/847/2009/The performances given by Henderson and Fread formulae (Eqs.5 and 8) are strongly dependent on the corrective coefficient r, which is a function of a so-called "typical" or reference wave for the concerned reach (see Eq. 6); since it is not possible to set a reference wave for the channels used in the simulations, r was computed for each case from the incoming wave characteristics.Such a procedure, although it produces good results in theoretical cases, can only be applied in natural rivers to reconstruct the flood hydrograph after the event has passed, and not for an operational on-line discharge measurement.
The improved performance of Perumal 2 formula (Eq.20) with respect to the others was also found by Perumal and Moramarco (2005), using similar numerical experiments.
In case 6 (fast wave over mild river bed slope), the accuracy of formulae based on single section measurements decreases significantly, as one can see from the observation of mean error values (Fig. 3) and from the hydrographs (Fig. 7); lastly, analysis of case 8 (slow wave over very mild river bed slope) shows that, in reaches with a very mild bed slope, none of the formulae using single water stage measurement is able to correctly estimate the discharge (Fig. 8).The results of   On the left: estimated and "true" discharges hydrograph for case 9 (channel with bed bottom slope 10 −4 , and variable prismatic cross section); on the right: estimated and "true" discharges hydrograph for case 11 (channel with bed bottom slope 10 −4 , and variable irregular section).
Fenton and Perumal formulae (Eqs.16 and 20) show a high level of noise.Nonetheless, in order to show a fair comparison of all the tested approaches, it was not felt appropriate to filter out noise because, to the other equations, the results were obtained using "perfect measurements".On the contrary, as will be discussed in Sect.3.3, the application of a filter to the DyRaC results, is an essential prerequisite for operational installations in order to eliminate the "measurement errors", which will inevitably affect the discharge estimates.
On the other hand, even in the presence of fast flood waves, formulae using simultaneous couples of water stage measurements, like the Chow (Eq.23) and DyRaC (Eq.33) formulae, provide accurate estimation, with a maximum error of the order of 1%.

Comparison on channels with spatially variable sections
The analysis of results in cases 1 to 8 shows that the Chow equation Eq. ( 23) and DyRaC equation Eq. ( 33) provide almost coincident results when dealing with prismatic channels.As mentioned earlier, in Figs. 1 to 8 Fenton and Keller formula Eq. ( 25) was not explicitly mentioned because it is identical to Chow equation in prismatic channels.However, as pointed out by Schmidt and Yen (2002), in natural rivers even Eq.( 25) may become incorrect if the longitudinal sec-tion variation makes the convective acceleration terms relevant.The magnitude of this term has been evaluated using both a channel with varying prismatic sections (case 9) and a channel with irregular sections (cases 10 and 11); Fig. 9 illustrates flood hydrographs for cases 9 (left) and 11 (right) and, as can be seen, only the Jones, Fenton and Keller and DyRaC formulae have been represented, along with the exact discharge and the steady flow rating curve.In both cases, unlike the DyRaC formula (Eq.33), the Fenton and Keller approximation (Eq.25) is not able to return the correct discharge hydrograph.Hence, it may be inferred that the parabolic approximation (Eq.22) used for the derivation of the Fenton and Keller equation Eq. ( 25), which implies neglecting both the convective and the local acceleration terms, can seldom be applied to discharge estimation in natural rivers unless the river reach involved is characterised by constant cross sections.

Influence of measurement accuracy on discharge estimation
The methodology described in Sect.2.5 has been applied to case 10, which uses irregular cross sections, to simulate a typical operational use of the DyRaC formula (Eq.33).
Figure 10 shows the resulting hydrograph compared to the "true" value and to the one derived from the steady-flow rating curve (top left); the values of I σ/µ the cut-off indicator, obtained at each time step (top right), the error rate (bottom left) and the number of measurements needed to reach the required precision of 5% of I σ/µ (bottom right).As can be seen, even when initialising the estimation process with a minimum number of samples (10) the required precision is automatically reached; only in a limited number of cases are more measurements necessary.Please note that in order to "filter" the water oscillations, measures should be taken at random time intervals, with an average delay ranging from 1 to 5 s.Therefore, the results obtained imply that even in the worst cases a discharge value can be operationally estimated in a couple of minutes.Also note that the estimation accuracy can be improved by increasing the number of initial samples; the graphs in Fig. 11 show the results obtained using a minimum of 20 samples for each time step: as can be seen, the error rate significantly decreases with respect to the previous example shown in Fig. 10.Although the described procedure should be operationally verified in real world applications, the results presented in this work are very promising and it is reasonable to believe that the DyRaC approach can be successfully applied in most natural rivers.

Operational discharge measurements under difficult conditions
As is evident from Sect.2.5, the application of most reviewed methods for real time discharge evaluation require a considerable amount of information regarding channel geometry, along with stage-discharge measurements for calibration.This means that the monitored river must have good channel stability, without significant deposition and erosion processes, and ease of access, in order to install and apply the necessary instrumentation.However, flow estimation in rivers located in impervious areas or characterized by strong sediment transport and braided channels appears to be a more complicated task, since the correct application of methods reviewed is no longer possible.In such situations, discharge estimation methodology needs to be based on simplifying assumptions and a limited amount of data.An example of simplified methodology is given in Petersen-Øverleir (2006); this study provides a method based on the Jones formula and nonlinear regression, which requires only stage-discharge measurements and a stage hydrograph.The regression model is developed by applying the monoclinal rising wave approximation and the generalized friction law for uniform flow, along with simplifying assumptions regarding the hydraulic and geometric properties of the river channel in conjunction with the gauging station.Experimental application of this methodology has provided good results.

Conclusions
Results obtained in the present work confirm the need to estimate discharge by means of expressions accounting for water surface slope, as stated by several authors (Henderson, 1966;Fenton, 2001;Schmidt and Garcia, 2001).Formulae not explicitly accounting for water surface slope can provide good estimations in kinematic or quasi kinematic conditions and, generally speaking, in channels with a steep bed slope (approximately 5×10 −4 or greater), while they perform poorly in other conditions, especially in the presence of fast flood waves over mild bed slopes.In these cases, particularly in reaches with variable or irregular cross sections, it is necessary to measure the water surface slope directly and use a methodology like the proposed Dynamic Rating Curve.Results obtained by this procedure have proven to be accurate and reliable in all the numerical experiments; however, it is important to stress that the application of formulae using simultaneous stage measurements is slightly more demanding, in that, apart from the knowledge of the stage in two adjacent cross sections, it also requires the description of two river cross section geometries and the use of a small piece of code.
Nonetheless, the DyRaC approach offers many advantages in contrast to the steady-flow rating curve: it not only takes into account the loop characterising unsteady flow, but it also drastically reduces the steady flow rating curve extrapolation errors at higher flow regimes.The steady-flow rating curve is generally fitted using measurements taken during low or medium flow regimes.Extrapolation beyond the range of measurements is essentially dominated by one parameter, an exponent, which controls the curvature of the rating curve; this produces a significant uncertainty in the extrapolation with large discharge estimation errors.On the other hand, in the DyRaC approach the curvature of the rating curve is correctly driven by the cross section geometry, which is known, while the evaluation of the roughness coefficient, which is the only required parameter, has a limited influence since it may be considered more or less constant at high flow regimes.This is why the DyRaC approach allows for an accurate calibration even when using stage-discharge measurements taken at low and medium flow conditions.
Finally, as found in previous works (Dottori et al., 2008), the DyRaC methodology also allows for accurate discharge estimation in sections affected by backwater effects, which is taken into account during the experimental stage-discharge measurements, used for roughness calibration.Application of the DyRaC approach to natural rivers will be presented in a forthcoming paper by the same authors.
At present, a measurement instrument based on DyRaC is under development and will be operationally installed and tested on several rivers presenting different hydrological characteristics and conditions.

Fig. 1 .
Fig. 1.Case 10: upstream (left) and downstream (right) cross sections in the channel reach where discharge has been estimated; distance between the two section is 1 km.

Fig. 2 .
Fig. 2. Case 6 (left) and 8 (right); time evolution of Ra, the ratio between local acceleration term (1/g) [∂ (Q/A) /∂t] and hydraulic head slope ∂H /∂x, expressed as a function of the rate of change in water depth ∂y/∂t.

Fig. 4 .
Fig. 4.Comparison of standard deviation of discharge estimation error for cases 3, 4, 5 and 6.Values of standard deviation error for Eq.(16) in cases 5 and 6 are larger than 200 m 3 s −1 .

Fig. 7 .
Fig. 7. Case 6 (channel with bed bottom slope 10 −4 , wave with a 24 h rising time period): estimated and "true" discharges hydrograph.Chow and DyRaC values coincide almost exactly with the "true" curve.

Fig
Fig.9.On the left: estimated and "true" discharges hydrograph for case 9 (channel with bed bottom slope 10 −4 , and variable prismatic cross section); on the right: estimated and "true" discharges hydrograph for case 11 (channel with bed bottom slope 10 −4 , and variable irregular section).

Fig. 10 .
Fig. 10. 10 with error affected stage measurements; top-left: estimated and "true" discharges hydrograph; top-right: computed values of the cut-off indicator I σ/µ ; bottom-left: normalised discharge estimation error (estimation error divided by "true" value); bottom-right: number of measurement samples needed to reach the required accuracy: the minimum number for each time step is set to 10.

Fig.
Fig. Case 10 with error affected stage measurements; top-left: estimated and "true" discharges hydrograph; top-right: computed values of the cut-off indicator I σ/µ ; bottom-left: normalised discharge estimation error (estimation error divided by "true" value); bottom-right: number of measurement samples needed to reach the required accuracy: the minimum number for each time step is set to 20.

Table 1 .
Characteristics of numerical experiments.In all the experiments, Manning's roughness has been set equal to 0.035 m −1/3 s.