Positive Surge Propagation in Sloping Channels

A simplified model for the upstream propagation of a positive surge in a sloping, rectangular channel is presented. The model is based on the assumptions of a flat water surface and negligible energy dissipation downstream of the surge, which is generated by the instantaneous closure of a downstream gate. Under these hypotheses, a set of equations that depends only on time accurately describes the surge wave propagation. When the Froude number of the incoming flow is relatively small, an approximate analytical solution is also proposed. The predictive ability of the model is validated by comparing the model results with the results of an experimental investigation and with the results of a numerical model that solves the full shallow water equations.


Introduction
A positive surge in an open channel is a sudden increase of flow depth that looks like a moving hydraulic jump; it can be generated by either natural or artificial causes [1][2][3][4][5].The early contributions on surges include the works of Favre [6], Benjamin and Lighthill [7], and Peregrin [8].Especially notable is the more recent work of Chanson and his co-workers, who mainly focus on tidal bore, i.e., a surge wave generated by a relatively rapid rise of the tide ( [9][10][11][12][13], e.g., among the many).
In the present work, we distinguish tidal bore from positive surge since the former is produced by a fast increase of the downstream water level, whereas the latter is produced by a rapid reduction, often to zero, of the downstream flow rate.The two phenomena share many similarities; however, differences are such that the two waves are not likely to be confused with each other.Here, the focus is on positive surge.
Despite the significant progress in understanding the dynamics of surges and bores in rivers and channels, some problems of practical interest still remain.Most studies of positive surges assume a horizontal channel, whereas many practical applications encompass positive surges propagating upstream on downward sloping channels; e.g., rejection surges in canals serving hydropower stations produced by the sudden decrease or complete interruption of power production, or surges in drainage channels ending at a pumping station when a sudden stopping of pumping occurs.
When a positive surge propagates upstream on a steep slope over an incoming supercritical flow, it progressively decelerates and asymptotically tends toward a stationary hydraulic jump [10]; when, instead, a positive surge propagates upstream on a mild slope over an incoming subcritical flow, the wavefront progressively reduces its height until the surge vanishes.Interestingly, the latter occurrence has never been reported and discussed in the literature and, as it happens, solutions to open channel flow problems never ceases to amaze by showing unexpected behaviours even within simple and schematic frameworks (e.g., [14][15][16][17][18][19][20]).
Numerical models can provide accurate solutions to the problem; however, analytical solutions or conceptual models can help knowing how flow and geometry parameters affect the surge propagation, and can provide quick answers to a number of practical questions.
In this paper, the problem of positive surge propagation against an either subcritical or supercritical uniform flow is faced within the framework of one-dimensional shallow water flow.Section 2 presents an approximate, zero-dimensional model for the description of the surge propagation.The capability of the model in predicting the characteristics of the surge propagation is validated through the comparison of model prediction with the results of laboratory and numerical experiments.Section 3 describes the laboratory setup and the experimental procedures, and the numerical model that solves the shallow water equations.The comparison between model prediction and the results of the laboratory and numerical experiments is discussed in Section 4. Conclusions of this study are shortly drawn in Section 5.

The Model
In a sense, the proposed model is a zero-dimensional model since it consists of a set of equations that depend only on time.The model, hereinafter referred to as 0D model, considers a rectangular channel and assumes the bottom slope, s small enough so that cos θ ≈ 1, sin θ ≈ tan θ = s, with θ the angle of the channel bed to the horizontal (Figure 1).
The incoming flow is assumed uniform with water depth Y 0 and bulk velocity U 0 .Downstream of the wave front, energy dissipation is neglected and the water surface is assumed horizontal.These assumptions are likely accurate during the early stage of the surge propagation when most of the incoming water goes back upstream transported by the surge and hence the flow velocity downstream of the front is negligibly small.The volume of water per unit width stored downstream of the wave front (the shaded area in Figure 1) is where ∆Y = Y 1 − Y 0 is the wave front height, L the front position, and s the bed slope.Since a = dL/dt, with a the surge celerity for an observer standing on the bank, positive upstream, and U 0 Y 0 = dW/dt, we find and The continuity and momentum equation across the surge, in the frame reference moving with the wave front and when neglecting the front acceleration and the local inertia, yields [21]: The set of Equations ( 2)-( 4) can be integrated in time to yield ∆Y, a, and L (a simple numerical integration procedure is given in Appendix A).In addition, the water depth, h(t), at the most downstream section (see Figure 1) is given by and it quickly grows approximately with sL since the reduction with time of ∆Y is relatively slow.
In order to make the model analysis more robust, we introduce suitable time and length scales for model variables and parameters.Both a and U 0 are scaled with the small wave celerity c 0 = gY 0 , the wave front height, as well as all vertical lengths, are scaled with the undisturbed flow depth Y 0 , whereas the front position L is scaled with Y 0 /s.With this, from Equation (2), a suitable timescale, τ, is found It is worth noting that τ actually measures the front propagation time since it is also given as the ratio of the scale length for the front position (i.e., Y 0 /s) to the scale velocity for a (i.e., c 0 ).Interestingly, with the above scaling, the solutions a(t/τ)/c 0 , ∆Y(t/τ)/Y 0 , and L(t/τ)s/Y 0 , turn out to depend only on the Froude number of the incoming flow, F 0 = U 0 / gY 0 .
When the incoming flow is subcritical, the surge progressively reduces its height while propagating upstream until the front vanishes.The distance, L M , travelled by the surge before vanishing can be estimated with Equation (2) that is here rewritten as When the wave front vanishes, at t = t M , ∆Y reduces to zero and the front velocity reduces to the absolute small wave celerity a = −U 0 + gY 0 so that The numerical solution of model equations for a wide range of values of model parameters confirms that the term between brackets at the denominator of Equation ( 8) depends only on F 0 .Therefore, Equation (8) can be written as An approximate equation for Φ(F 0 ) that fits well the results of the model (see Figure 2) is When the incoming flow is supercritical, the surge progressively decelerates and asymptotically tends toward a stable stationary hydraulic jump [10,22], the surge celerity reduces to zero gradually so that the distance L grows to infinity.

An Approximate Analytical Solution
When F 0 is relatively small (i.e., approximately F 0 < 0.6), which is the case in many practical problems, the time behaviour of the wave front velocity is extremely well approximated by a second order polynomial as where a 0 is the wave front velocity at t = 0, and k 1 , k 2 are two coefficients that depend only on F 0 .Since the front position at t = 0 is L = 0, a 0 can be computed by combining Equation ( 2) with Equation ( 4) to eliminate ∆Y.By solving the 0D model for a large set of F 0 values, the following approximate equations for k 1 and k 2 are found (see Figure 3a): Once a(t/τ) is computed with Equation (11), ∆Y(t/τ) and L(t/τ) can be computed with Equations ( 4) and (3), respectively.Figure 3b,c compares the numerical solution of the 0D model to the approximate analytical solution for F 0 = 0.65 and F 0 = 0.6, respectively.When F 0 = 0.6, the differences between the numerical and analytical solutions are negligibly small; at smaller Froude number, the numerical and analytical solutions are indistinguishable.

Numerical and Laboratory Experiments
In order to check the predictive capability of the 0D model, the model results are compared to the results of an experimental investigation as well as with the results of a numerical model that solves the full shallow water equations, for a wide range of flow and geometric conditions.
Laboratory experiments are performed in a flume that is relatively short if compared to the distance that a surge can normally cover; the results of these experiments are used to check the capability of the model in predicting the first stage of formation and propagation of the surge.The numerical experiments are then used to test the model over the longer distances.

Laboratory Setup and Experimental Procedures
The experiments are performed in a 20 m long, 0.38 m wide, tilting flume with Perspex walls and bottom.Water is recirculated via a constant head tank that maintains steady flow conditions, and the flow rate is measured by a magnetic flowmeter (accuracy of about 0.2%).In order to achieve uniform flow, a downstream weir can adjust water depth when the flow is subcritical, whereas an upstream, vertical sluice gate adjusts the supercritical water depth.A vertical sluice gate, approximately 0.5 m upstream of the downstream gate, is used to generate the positive surge by quasi-instantaneously stopping the flow (frame by frame analysis of some gate closures recorded with a camera at 12.5 Hz frame rate allowed to estimate the gate closure time that is less than 0.2 s).
We use five ultrasonic probes, calibrated against a pointer gauge, to measure both uniform and unsteady water depth (accuracy of about 0.13 mm, acquisition frequency = 100 Hz) [23].The five probes (S 1 to S 5 ) are located at the channel centreline and at X S = 0.15 m , 3.0 m, 5.96 m, 8.96 m, and 11.89 m, upstream of the sluice gate used to generate the surge.The bed slope is varied between s = 0.0002 and s = 0.0052 in order to vary the Froude number of the incoming uniform flow in the range 0.3 < F 0 < 1.4.Each run was repeated three times with negligible differences in the recorded water depth.Some preliminary experiments, aimed at gaining confidence with the phenomenon, are also performed in a smaller flume, 0.3 m wide and 6 m long, equipped the same as the the larger flume.Experimental conditions are summarized in Table 1 for both the flumes.
Table 1.Summary of experimental conditions.Runs 1 to 41 are performed in the long flume, runs 42 to 73 are performed in the short flume; F S is the surge Froude number defined by Equation (13).

The Numerical Model
Although the problem is basically one-dimensional, since available, we use a two-dimensional horizontal (2DH), shock-capturing hydrodynamic model based on the Finite Volumes technique.
The numerical model solves the depth-averaged shallow water equations, written in conservative vector form, using a Godunov-type method on unstructured triangular grids [24][25][26].To properly deal with sloping channel, the model uses a second-order accurate description of the channel bed, i.e., bottom elevations are defined at the grid nodes and are assumed to vary linearly within each element of the mesh [27,28].The model is made adaptive inasmuch as data at cell faces are reconstructed selecting either primitive or conservative variables according to the local Froude number [29], thus avoiding the generation of unphysical discontinuities at cell faces even in the case of smooth flows.First-order adaptive schemes, which utilize a second-order accurate description of terrain, were shown to be robust, efficient, and accurate for many engineering applications [29].In addition, the model performance is increased in terms of computational cost by using the Local Time Stepping method as described in [30].
Simulations are performed over a rectangular domain, formally 1.0 m wide, having a length, L M , in the range 100-2000 m depending on the distance travelled by the surge before vanishing or before turning into a quasi-stationary hydraulic jump.The number of triangular computational elements is kept constant and equal to 4000; hence, the size of grid elements in the longitudinal direction is L M /2000 m.

Comparison between the 0D Model Prediction and the Laboratory Experiments
Laboratory experiments are used to check the 0D model assumptions and its capability in predicting the first stage of formation and propagation of the surge.
Figure 4 shows two examples of free surface elevation (rather than flow depth) over time recorded by the five ultrasonic probes; free surface elevation is constant before the surge reaches the generic flume section S i , equipped with a probe, and, after the surge passage, all free surface elevations follow the same temporal path and overlap.This behaviour, which is found in all of the experiments, confirms that the free surface downstream of the surge front is actually nearly horizontal, as assumed by the 0D model.Figure 5 compares the measured dimensionless flow depth, Y/Y 0 , as it varies with dimensionless time t/τ, with the results of the 0D model for some of the experimental runs.The results of this comparison are actually very good, especially considering that the proposed model has no calibration parameters.The model accurately predicts the arrival time of the surge front at the gauged cross sections although, as a general trend, it slightly underpredicts the arrival time, and hence the surge celerity.Discrepancies between predicted and measured surge celerity are, however, very small.
The average height of the surge front is also sufficiently well predicted even when the surge Froude number, defined to be [31] F S = U 0 + a gY 0 (13) is small, and the surge is an undular wave [32,33].In addition, the rate of water level raise behind the front is satisfactorily predicted by the model.The greater differences between model predictions and experimental results are found when the incoming flow rate and the uniform flow depth are very small (e.g., Figure 5e); in this case, the Reynolds number, Re = U 0 Y 0 /ν, with ν the kinematic viscosity, is small (e.g., Re = 4730 for flow conditions of Figure 5e) and viscosity is likely to affect the flow.
It is also interesting to observe in Figure 5e that the front height reduces significantly during the surge propagation (∆Y/Y 0 ≈ 1.5 at X S1 , ∆Y/Y 0 ≈ 0.6 at X S5 ) and this behaviour is correctly predicted by the model.
As a general trend, the model slightly underpredicts surge wave celerity and overpredicts the surge height compared to experimental results; indeed, differences between experimental results and model predictions increase with the increasing of channel bed slope.From the comparison between model predictions and experimental results, we can conclude that the model performs satisfactorily at least during the early stage of surge formation and propagation.
We also use experimental data to assess the type and characteristics of the surge.Figure 6, which also includes data from the literature, shows the relative maximum depth at the first wave crest, Y max /Y 0 , as a function of the surge Froude number, F S .The present experimental results use the data measured by the third (S 3 ) and fourth (S 4 ) water level probes, in the large flume, and by one probe, located approximately 1.0 m upstream of the operating sluice gate, in the small flume.
At large F S (F S 1.6), Y max /Y 0 grows approximately as the relative conjugate depth Note that Equation ( 14) slightly underestimates Y max /Y 0 possibly because it neglects the effect of gravity [16].Breaking surge typically occurs in this range of Froude numbers.
At smaller F S (approximately, F S < 1.25 − 1.30), Y max /Y 0 grows nearly linearly with F S ; interestingly, if we denote with η = Y max − Y 0 the maximum wave front height and assume η/∆Y = 1.8, with ∆Y = Y conj − Y 0 , the resulting curve plotted in Figure 6, fits very well with the experimental data.In this range of Froude numbers, undular surge is observed.Relative maximum depth at the first wave crest, Y max /Y 0 , as a function of surge Froude number, F S .The black curve denotes the conjugate depths ratio as given by Equation ( 14); the dashed curve is from Equation ( 15), the dot-dashed line denotes the breaking criterion for a solitary wave (Y max − Y 0 )/Y 0 = 0.78.
In the range 1.25-1.30< F S < 1.6, the transition from undular to breaking surge is observed; close to the lower boundary of this interval, the surge looks undular but some spilling can be observed down the crest of the wave; at larger F S , close to the upper boundary of the interval, the surge looks breaking, but a train of secondary waves is formed behind the front.These results agrees with most of the observations reported in the literature (e.g., [6,11,12,16,34,35]).

Comparison between Model Prediction and the Numerical Experiments
Numerical experiments are mainly used to check the predictive capability of the 0D model for surges propagating over long distances (i.e., longer than the available experimental flumes).Conditions of numerical simulations are summarized in Table 2.It is worth saying that the 0D model only allows evaluating the timing of the surge front and the average water level increase downstream of the front; water waves at (and behind) the front, typical of undular surges, cannot be captured by the equations that assume hydrostatic pressure and neglects free surface slope and curvature.
Examples of the comparison between the 0D model prediction and the numerical solution are shown in Figure 7.The position of the front with respect to time, as provided by the 0D model (black solid lines) and the numerical model (red dashed lines) respectively, compare favourably for most of the tested conditions (Figure 7a-e).
A general good agreement is found also in terms of front height (Figure 7f-h) and celerity (Figure 7i-m).Minor deviations (Figure 7f,i,j) can be observed for the smaller Froude numbers, F 0 .Specifically, the model underpredicts the distance travelled by the surge before the front vanishes; however, differences between model predictions and numerical results are non-negligible only at the very late stage of the surge propagation, when the front height is largely reduced and close to vanishing.At the larger Froude numbers (Figure 7g,h,k,l,m), the agreement between model predictions and numerical results is satisfactory.
The reason for the discrepancy that is observed at low Froude numbers and when the front is close to vanishing is related to one of the main assumptions of the 0D model, i.e., that the free-surface remains horizontal downstream the front.Typical free surface profiles at different times during the surge propagation are shown in Figure 8.The free surface behind the front is actually horizontal during the first stage of surge propagation, when the front height is relatively large and the flow rate required to fill the volume of water behind the front is relatively small so that the flow velocity just behind the front, and hence energy dissipation, are negligibly small.As the front progresses, its height and celerity reduce and, importantly, the flow rate just behind the front gradually increases to eventually meet the incoming flow rate.In this way, energy dissipation behind the front approaches energy dissipation of the uniform flow and, for relatively small Froude numbers (Figure 8a), the free surface slope increases toward bed slope.Because of this behaviour, when the incoming flow is subcritical and the Froude number is relatively small, the front height reduces to zero asymptotically.The results of the numerical simulations show that the deviation from horizontal of the free-surface downstream of the front decreases with bed slope increasing, and finally reverses when F 0 is relatively large (see, e.g., Figure 8b).The concavity of the free-surface profile downstream of the surge is likely related to flow unsteadiness; as water depth increases in time, flow rate per unit width decreases from U 1 Y 1 , just downstream of the front, to zero at the gate section.This condition is somehow similar to the case of a distributed outflow from a channel when the flow is steady and subcritical; in such a case, the free-surface profile is known to have a positive gradient in the flow direction and a negative curvature [21].
Interestingly, when the incoming flow is supercritical and the surge tends to a stationary hydraulic jump (Figure 7h,m), the main model assumption (i.e., horizontal free surface downstream the front) always holds.Indeed, the front height does not reduce to zero and the velocity behind the front is thus considerably smaller than the incoming flow velocity.

Conclusions
In this work, a simple zero-dimensional model to simulate the propagation of a positive surge in a sloping channel is proposed and its predictive capability is checked against experimental and numerical data.The studied surge is generated by the instantaneous closure of a downstream gate.
When comparing the 0D model results with experimental data, small discrepancies are observed, which are mainly caused by the strongly nonlinear effects produced by free surface slope and curvature.As a general trend, the model slightly underpredicts the surge wave celerity and overpredicts the surge height; differences between experimental results and the 0D model predictions increase with the increasing of channel bed slope.
From the comparison between the numerical results and the predictions of the proposed model, it emerges that the 0D model can accurately simulate the initial stage of the surge formation and propagation.During the late stage of surge propagation, and when the Froude number of the incoming flow is small, so that the fate of the front is to vanish, the assumption of horizontal free surface behind the front stops holding true and the surge front reduces its height slower than that predicted by the model.However, this discrepancy between model predictions and numerical results is possibly unimportant from a practical point of view since it occurs during the late stage of surge propagation when the front height is largely reduced.
As a whole, the proposed 0D model, and even more so its analytical solution for F 0 < 0.6, are simple tools that allow for a quick preliminary estimate of flow conditions following the development of a positive surge in a sloping channel and are of help in the design process as they clearly show the impact of channel geometry and flow conditions on the characteristics of the surge.

Figure 1 .
Figure 1.Schematic of the surge propagating upstream against a uniform flow of depth Y 0 and velocity U 0 , with notation.

Figure 2 .
Figure2.Plot of Φ(F 0 ) as a function of F 0 .White circles denote the results of the 0D model, the black line is the approximation(10).The inset shows the same plot, with Φ(F 0 ) on the log-scale.

1 =Figure 3 .
Figure 3. (a) coefficients k 1 and k 2 as they vary with F 0 (symbols denote the numerical solution of the 0D model, full lines are the approximation given by Equation (12)); (b) wave front velocity, a, and height, ∆Y, and the position of the front, L, as a function of time for F 0 = 0.65 (full lines denote the numerical solution of the 0D model, dashed lines denote the analytical solution); (c) same as (b) for F 0 = 0.60.
Figure 6.Relative maximum depth at the first wave crest, Y max /Y 0 , as a function of surge Froude number, F S .The black curve denotes the conjugate depths ratio as given by Equation (14); the dashed curve is from Equation (15), the dot-dashed line denotes the breaking criterion for a solitary wave (Y max − Y 0 )/Y 0 = 0.78.

Table 2 .
Summary of numerical conditions.k St is the roughness coefficient in the Strickler formula for uniform flow.