Estimating the long-term evolution of river bed levels using hydrometric data

The stage-discharge measurements and rating curves accumulated over decades at hydrometric stations are a valuable source of information on the long-term evolution of river bed levels. However, the methodology to extract meaningful geomorphic information from such hydrometric data is not straightforward. We introduce an original method to estimate the parameters of successive rating curves by Bayesian analysis in sequence. These parameters reflect the physical properties of the channel features that control the stage-discharge relation: low-flow riffles, main channel, floodway (bars), floodplain, etc. The dates of rating changes are assumed to be known in existing hydrometric records. The uncertainty interval of each parameter is estimated, assuming, however, that no rating change has been ignored by the station manager. It is thus possible to clearly distinguish overall trends of the channel bed level from the local evolution of riffles and to evaluate whether the observed temporal changes are significant compared to the estimation uncertainties.


Introduction
The multi-decadal evolution of the geometry of river beds is only rarely monitored, typically with topographic surveys (cross-sections, longitudinal talweg or water-surface profiles). Even when such data exist, their analysis poses problems of spatial and temporal representativeness: even when cross-sections are surveyed at the same location, they may not be representative; talweg is mobile and not fully representative of the average bottom level; water-surface profiles are measured for different flows, etc. The long-term morphological evolution of river beds is then difficult to quantify precisely.
Luckily, the gaugings and rating curves accumulated over decades at hydrometric stations are a valuable source of information on the evolution of river beds. A first use is to explore for temporal change in the water stage corresponding to a fixed flow rate, for the series of rating curves available. The choice of the reference flow is important because it determines the nature of the morphological feature the evolution of which is monitored: for a low flow rate, the stage will generally be determined by the width and height of a local cross-section, whereas for an intermediate flow it will be representative of the average level of the main channel. In addition, the stage for an overbank flow will also account for the level of the floodplain. The advantage over topographic surveys is that hydraulically meaningful parameters can be estimated. This paper introduces an original method that has been developed to estimate the parameters of the successive rating curves by Bayesian analysis in sequence. This work is a development of the BaRatin method [1] available for the Bayesian estimation of unique stage-discharge rating curves. These parameters no longer depend on a reference flow to be fixed, but are linked to the various hydraulic controls identified: low-flow control riffle, main channel uniform flow control, floodway (bars), floodplain, etc. The dates of rating changes are taken from existing hydrometric records. The uncertainty interval of each parameter is estimated, assuming, however, that no rating change has been ignored by the station manager. It is thus possible to clearly distinguish evolution of the reach of channel adjacent to stage recorder (e.g., a degrading bed) from that of local riffles and to evaluate whether the observed temporal changes are significant compared to the estimation uncertainties. The method is illustrated with two unstable rivers in New Zealand, each subject to different natural or anthropogenic pressures.

Morphogenic changes in hydraulic controls
Commonly in natural rivers, the stage-discharge relation is controlled by critical sections over downstream riffles at low flows and by friction in more or less uniform channels at higher flows. At sites where a two-segment rating curve applies, a simple model relating discharge to stage ℎ is therefore The water elevation at which the hydraulic control changes is called the breakpoint. This parameter is estimated through the numerical solution of the continuity condition applied to the stage-discharge model: ( ) = 1 ( − 1 ) 1 = 2 ( − 2 ) 2 .
If the low-flow segment of the rating curve is determined by a critical-flow section that can be approximated by a rectangular weir equation, then 1 = √2 with [-] a discharge coefficient, [m] the weir perpendicular-to-flow width, [m/s²] the gravity acceleration, 1 [m] is the mean elevation of the weir crest and = 3/2 [-] is the exponent.
The second segment of the rating curve can typically be modelled using the Manning equation for a wide rectangular channel, for which the hydraulic radius is equal to the flow depth. The parameters of this channel control are 2 = √ 0 where [m 1/3 /s] is the Strickler flow resistance coefficient (inverse of Manning's ), [m] is the channel width, is the mean bed elevation of the controlling reach across the stage-recorder section and 2 = 5/3.
Morphogenic floods will cause changes in the stage-discharge model and in the values of some of its parameters [2]. We focus on the common situation in which the number and types of controls does not change, i.e. Eq. 1 remains valid over all the successive periods divided by the occurrence of morphogenic floods, and the values of exponents 1 and 2 are constant. It may happen that the values of coefficients 1 and 2 are varying because of changes in the width of the cross-sections. However, the most common changes affect the values of the offsets 1 and 2 , which is often reflected by "shifts" applied by the field hydrologists to their rating curves. The offset 2 is modified by overall changes in the river bed level, i.e. degradation or aggradation of the whole controlling reach (cf. Figure 1b). In addition to overall changes, the offset 1 is also modified by local changes, i.e. reshaping and relocation of the corresponding critical cross-section (cf. Figure 1c). Therefore, the time evolutions of offsets 1 and 2 are linked but the interpretation of 2 changes is simpler as it is representative of the reach-average bed evolution.

Bayesian estimation of rating curves
The parameters of the stage-discharge model (cf. Eq. 1) have to be estimated for each period. The traditional methods for developing stage-discharge rating curves include graphical fits against stage-discharge measurements, called "gaugings", combined with more or less explicit expert knowledge on the active hydraulic controls. Such methods would be impractical for the re-analysis of large sets of rating curves in the long-term over several hydrometric stations. On the other hand, automated calibration techniques relying on data only would be insufficient as observational data sets (gaugings) usually do not contain enough information to constrain the problem. Bayesian techniques allow the automated combination of formalised expert knowledge and observational data, accounting for their respective uncertainties and contents of information. Typically, the BaRatin method [1] is composed of three main steps (cf. Figure 2): the establishment of the stage-discharge model and the prior distribution of its parameters, based on the available information (except the gaugings) on hydraulic controls; the inference of the posterior distribution through Bayesian analysis, using both the prior distribution and the likelihood computed from the gaugings and their individual uncertainties; and Markov chain Monte Carlo (MCMC) sampling of the posterior distribution of the rating curve parameters and computation of probabilistic rating curves and streamflow time series. The uncertainty intervals at any probability level can be computed from the samples drawn in the last step.
The existing BaRatin method assumes a stable stage-discharge relation. To apply it to a sequence of periods between morphogenic floods with parameter changes at known times, the following procedure was set up. The parameter results obtained from the estimation of the previous rating curve were recycled to be used as priors for the estimation of the next rating curve. The prior distributions were assumed to be Gaussian and interactions between parameters were ignored. The variances of the prior distributions of the next rating curve were set equal to the variances of the posterior distributions computed for the previous rating curve, plus an additional variance for parameters that may have changed (the offsets, typically). The additional variance was constant for all the periods and it was estimated a priori to reflect the largest possible changes to be expected. In this approach, information on the rating curve parameters is transferred to a given period from all the previous ones. The intention of this is to constrain the values of stable parameters using the information of the gaugings of previous periods, while the dynamic parameters are relatively free and can be estimated using the gaugings of the new period.

Shotover River at Bowens Peak, New Zealand
The Shotover River is a fast-flowing river of the Southern Alps, in the South Island of New Zealand. At Bowens Peak, its last hydrometric station before the junction with the Kawarau River, the catchment area is 1079 km² and the altitude (at gauge) is 320 m a.s.l. The Shotover channel at Bowens Peak (cf. Fig. 3) is constricted in a canyon and its stagedischarge rating has constantly degraded in steps, totalling approximately 2.5 m over 50 years , as shown by 735 gaugings (14.7 gaugings per year on average). The station was included in the rating shift occurrence study by Ibbitt and Pearson [3]. The station is still actively monitored, but was relocated a short distance downstream on 1 st July 2012.
The cause for the constant degradation has not been firmly established yet. The alluvial fan formed near the confluence of the Shotover and the Kawarau, which drains Lake Wakatipu, possibly evolved, lowering the base level of the Shotover. Another possible, not exclusive, explanation might be that the river is still adjusting after the large amounts of sediment introduced from hydraulic sluicing by gold miners. Indeed, after it was discovered that the Shotover was one of the richest gold-laden rivers in the world, an intensive gold rush started in 1862.
Our field observations and information from the station managers suggest that a single channel control is active for all flows. This takes the form of a gravel bedded normal flow reach that extends approximately 200 m upstream and 200 m downstream of the stage recorder and has been slowly degrading. The effective width of this channel control has been fairly constant over time due to the very steep, hard-rock banks of the canyon. Therefore, a single-segment stage-discharge model was considered for that site: (ℎ) = 2 (ℎ − 2 ) 2 , (cf. Eq. 1) where offset 2 , i.e. the mean bed elevation, is the only variable parameter of the rating curve. The rating curves for the 269 periods (1967-2014) defined by the station's managers were re-estimated using the sequential application of BaRatin explained above. Prior distributions of the parameters were specified according to the BaRatin method [1,2] using available information on the geometry and roughness of the channel. The 95% discharge uncertainty of the gaugings was assumed to be ±7%. The computation lasted for about 1 hour on a normal laptop. The time evolution of the 2 estimates and their uncertainty intervals is presented in Figure 4. The estimated rating curve of each period fits well with the gaugings and the uncertainty intervals are much narrower than the quantified changes. As expected, the river bed shows a degrading trend but several periods of a few years each show substantial aggradation. Further investigation would be necessary to clarify if these aggradation periods may be related to enhanced sediment inputs from the catchment, e.g. floods or debris flows happening in active tributaries or headwater catchments.

Pohangina River at Mais Reach, New Zealand
The Pohangina River, in the North Island of New Zealand, is a highly dynamic gravelbed river with wandering planform [4]. This is a transitional pattern between multi-thread braided and single-thread meandering channels [5]. The slope is roughly 4 mm/m. The main channel periodically changes its path in constant interaction with secondary channels, abandoned channels, chutes and in-stream bars. Avulsion and dramatic changes can occur during morphogenic floods, as was the case during the February 2004 flood in the Manawatu catchment due to a '150-year' storm event [4].
The rating shifts of the Pohangina at Mais Reach station (operated by Horizons Regional Council) were also studied by Ibbitt and Pearson [3]. The hydrometric archive provides a very dense set of gaugings (717 gaugings between 1969 and 2015, 15.6 gaugings per year on average) generally showing riverbed degradation over time (cf. Fig. 5). Very high flows were not gauged after the 1970's and it is unclear how much the high-flow control may have changed. The transition between low-flow and high-flow segments appears to occur around 2 m at gauge, but this breakpoint elevation is likely to change substantially due to the complex morphological changes of the bedforms and riverbed in that river. As a first approach, the two-segment stage-discharge model (Eq. 1) was considered and the controlling channel widths were assumed to remain constant over time. Such approximations would ideally need to be elaborated to better reflect the complexity of such an avulsing river. However, the estimated rating curves were found to be in acceptable agreement with the gaugings for most of the periods. It was enough to demonstrate the ability of the method to estimate the parameters of the two-segment rating curves of the 273 periods  found in the hydrometric records, in a reasonable computing time (about one hour). Again, prior distributions of the parameters were specified according to the BaRatin method [1,2] using available information on the geometry and roughness of the channel, and the 95% discharge uncertainty of the gaugings was assumed to be ±7%.
The results are presented in Figure 6. As expected, offsets 1 and 2 show correlated time evolutions with a general degradation trend, but the trend is much less pronounced than in the Shotover case. Again, the hydrometric records provide a detailed view of continuous bed changes that are much more complex than the overall degradation trend and might be related to morphogenic flood events and local processes like bank retreat, gravel bar migration or channel avulsion. The channel control elevation ( 2 ) is generally slightly lower than the low flow control elevation ( 1 ). This is because the bed rises towards the critical cross-section, as typically happens transitioning out of a run/pool to a riffle. As it reflects both overall and local channel changes, the low flow control elevation ( 1 ) logically shows more variability around the general trend than 2 .
Last, the breakpoint ( ) between the two segments of the rating curve (cf. Eq. 1) shows a much larger scatter. The large variations of seem to magnify the scatter of 1 estimates. There are difficult to interpret physically but they may reflect the limited accuracy of the twosegment stage-discharge model for a such a complex river system. Using a third channel control to account for the floodway over the alternate bars might improve the results. Also, changes in coefficients 1 and 2 could be allowed since the widths of the riffles and of the main channel may substantially change at such a dynamic site.

Conclusions and perspectives
Hydrometric data provide a continuous and reliable source of information on the morphological evolution of river beds at many hydrological stations over decades. The proposed approach allows an automated estimation of identifiable morphological parameters, with their uncertainties. This type of analysis could be applied to large sets of hydrometric stations to provide a synoptic analysis of the morphological evolution of river bed levels throughout a watershed or a country. However, the accuracy of the results will be lower for hydrometric stations that come with less frequent gaugings and possibly more missed rating changes than the typical New Zealand sites presented in this study. For such cases, it would be even more useful to combine the analysis of hydrometric data with information from topography surveys and geomorphological studies along the river network.
It would be interesting to compare the results with those of more traditional approaches, e.g. the time evolution of the stage corresponding to a characteristic flow (the median flow, for example). A limitation of the sequential application of BaRatin is that parameter information transfer is time-oriented, from previous periods to the next, which deprives the first periods of information from future gaugings. This could be solved by using the BaRatin-SPD multi-period stage-discharge model developed by Mansanarez [2]. Also, we are still working on the detection of rating changes, instead of assuming they are perfectly known in hydrometric records.
NIWA (Kathy Walter) provided the data and information for the Shotover at Bowens Peak. Horizons Regional Council (Brent Watson) provided the data and information for Pohangina at Mais Reach.