Transient Flow in an Open Channel Bound by Two Step Pumping Stations

Abstract

Pumping stations play a great role in open channel flow applications. After starting pump units in the pumping stations, unsteady flows in the open channel are immediately generated. In this paper, the behavior of unsteady flow in a prismatic trapezoidal channel between two step pumping stations is studied. A mathematical model is constructed to simulate one-dimensional, unsteady, gradually-varied open channel flow in the channel. The soil leakage and evaporation from the water surface along the channel are also involved. The Lax diffusive scheme is applied to solve the proposed model. The model is applied on a real open channel located in China. The accuracy of the model is calculated by varying the scheme grid steps. We also calculate and discuss the water surface elevation, the flow rate and the lateral outflow per unit length, as well as several influential factors at various stages along the channel. In this work, the wave propagation mechanism is clearly presented and analyzed. The computed results show that the water surface profile in the open channel varies continuously over time until it becomes almost constant. Thus, the pumping station operation duties change over time. The accuracy of the results is validated by comparing the computed results with measured data.

1. Introduction

The study of open channels is an important topic involving hydraulics and hydrology. The flow conditions in real-life systems usually vary with time, and thus, the flows are unsteady. The analysis of unsteady flow is usually more complex than that of steady flow because unsteady flow conditions may vary with respect to both space and time. After the start-up of pumping station units to meet the flow rate requirements, transients in the open channels are produced and the flows become unsteady; in addition, the depth of transient flow and the pumping station’s operating duty influence each other. Therefore, the situation during the start-up phase is highly complicated.

The governing hydraulic equations, known as the Saint-Venant equations, are derived from the incompressible Navier–Stokes equations [1,2] which consist of continuity and momentum equations [3]. These equations are nonlinear, first-order, hyperbolic, partial differential equations for which closed-form solutions are not available, except for under certain simplified conditions; therefore, these equations are solved using suitable numerical techniques [4,5]. A wide range of numerical schemes based on the finite difference [6,7,8,9,10,11], finite element [12,13,14], and finite volume [15] methods have been applied to solve the open channel flow equations.

Several authors have applied the Saint-Venant equations or the shallow water equations for hydrodynamic simulations in rivers, channels, floodplains, lakes, and estuaries. Ding et al. [16] applied a numerical model based on the coupling of implicit and explicit solution algorithms of the shallow water equations to two large rivers in the reach of the Yangtze river in China and validated the accuracy of the numerical model for measuring both steady and unsteady flows using the prototype’s hydrological data. Aricò and Nasello [17] proposed two numerical models for the solution of the complete and zero-inertia forms of the shallow water equations and applied the proposed models to several unsteady flow routing scenarios. Akbari and Firoozi [18] presented two numerical schemes, namely; the Preissmann implicit scheme and the Lax diffusive explicit scheme, for the numerical solution of the Saint-Venant equations that govern the propagation of flood waves in natural rivers and compared their results with the results of the Hydrologic Engineering Center - River Analysis System (HEC-RAS) model. Machalinska-Murawska and Szydłowski [19] presented and analyzed two explicit schemes for the numerical solutions of one-dimensional Saint-Venant equations describing free surface flow, investigated the applicability of the Lax–Wendroff and McCormack schemes for modeling unsteady, rapid and gradually-varied open channel flow and compared their results to the arbitrary solutions. Viero et al. [20] presented a simplified model for the upstream propagation of a positive surge in a sloped rectangular channel, proposed an approximate analytical solution for situations where the Froude number of the incoming flow is relatively small, and validated the predictive ability of the model by comparing the model’s results with the results of an experimental investigation and with the results of a numerical model that solves the full shallow water equations. Gu et al. [21] applied a smoothed particle hydrodynamics approach to solve shallow water equations to study practical dam-break flows. Liu et al. [22] developed a new lattice Boltzmann approach to solve one-dimensional Saint-Venant equations, where the Taylor and Chapman–Enskog expansion was performed on the lattice Boltzmann equation; this model was verified by three benchmark tests.

To solve the governing hydraulic equations for unsteady open channel flow, Chun and Merkley [23] formulated an ordinary differential equation algorithm to approximate the solution to the characteristic form of the Saint-Venant equations for one-dimensional, gradually-varied, unsteady open channel flow in prismatic irrigation canals. Ferreira et al. [24] evaluated the difference in results between distinct procedures for solving the Saint-Venant equations, verified the Lax diffusive scheme’s results and compared them with the results of the HEC-RAS model. Szydleowski [25] presented explicit and implicit finite volume method schemes of the Roe type to model extreme unsteady, rapidly varied, open channel surface flow described by the Saint-Venant equations. Kalita and Sarma [26] presented the results of two different schemes, the Lax diffusive scheme and Beam and Warming scheme, and compared their results with the results of the MIKE21C modeling tool. They observed that both schemes led to the same results, and the Lax diffusive scheme took a shorter computational time for simulation. Therefore, the Lax diffusive scheme can efficiently solve the Saint-Venant equations, and it will be used in this work.

The purpose of this study is to estimate the water surface profile after start-up of the pumping station units in a trapezoidal channel, where the channel here is bound by two pumping stations. The flow rate and the lateral outflow of the transient flow are also evaluated. The Lax diffusive explicit finite difference scheme is applied to solve the Saint-Venant equations. We apply our model to a realistic channel allocated between Huai’an pumping station and Huaiyin pumping station. To validate the accuracy of the model, we compare the computed results with measured data.

The rest of the paper is organized as follows. Section 2 describes the two-step pumping station system. In Section 3, we present a mathematical model for one-dimensional, gradually-varied, open channel flow, with specification of boundary conditions and initial conditions. The Lax diffusive finite difference scheme is used to solve the mathematical model in Section 4. A prismatic trapezoidal channel bounded by Huai’an pumping station and Huaiyin pumping station in Jiangsu Province China is taken as a case study in Section 5. In Section 6, the computed results are presented in adequate figures and analyzed. To test the accuracy of the model, a comparison between the computed results and measured data is also made in Section 6. We end with conclusions and future works in Section 7.

2. Two-Step Pumping Station System

2.1. Definition of Two-Step Pumping Station System

In this work, we refer to a system consisting of two pumping stations—step 1 and step 2 pumping stations—and an open channel bound between them by a two-step pumping station system. Figure 1 illustrates this system in more detail; $y_1$ and $y_2$ denote the altitudes of the water levels from the sea level for the step 1 pumping station upstream and the step 2 pumping station downstream, respectively (in m), while $h_1$ and $h_2$ denote the water level of the step 1 pumping station downstream and the step 2 pumping station upstream, respectively (in m). In this study, $h_1$ and $h_2$ are considered to be constants. $H_1$ and $H_2$ are the heads of the step 1 and step 2 pumping stations, respectively (in m), while $q_v$ and $q_s$ are the evaporation and seepage losses per unit length along the channel (in m² s⁻¹).

In this system, the two pumping stations are connected in series. Each step pumping station consists of a number of identical pump units working in parallel mode. The role of these pump units is to pump water from a low area to a high position, and transfer water by a channel to the next step pumping station.

2.2. Background of the Case Study

In our case study, we considered a prismatic trapezoidal channel bound by Huai’an pumping station (Figure 2) and Huaiyin pumping station (Figure 3) in Jiangsu Province, China.

Huai’an pumping station is located at the intersection of the Grand Canal and Northern Jiangsu Irrigation Canal. It works with the Yunxi power station and irrigation diversion hub and other projects to form a Huaian water control project. The construction of Huai’an was started in December 1972, and it was put into operation in 1974. The auxiliary system includes a water supply and drainage system, a lubrication oil and pressure oil system and a compressed air system. The automation of the pumping station is high, and automatic monitoring, automatic recording and an automatic alarm have been implemented. This not only ensures the safe operation of the station, but also improves its economic benefit. It has been appraised as a high-quality project in the province.

Table 1. Descriptions of Huai’an and Huaiyin pumping stations.

Items	Huai’an Pumping Station	Huaiyin Pumping Station
Pump type	4500ZLQ60-4.89	ZL30-7-S
Designed head (m)	4.89	4.7
Number of pump units	2	4
Single-machine flow rate (m3 s−1)	60.0	30.0
Total flow rate (m3 s−1)	120.0	120.0
Single-machine capacity (kW)	5000	2000
Total capacity (kW)	10,000	8000
Motor type	TL5000-64	TL2000-48/3250
Rated speed (rpm)	100	125
Impeller diameter (m)	4.5	3.1
Vertical transmission mode	Vertical direct transmission
Inlet passages, discharge passages	Elbow type, syphon type

provides a brief description of the Huai’an pumping station.

Huaiyin Station is the third-stage pumping station for the Water Transfer Project from Yangtze River to the north. It is located at the north embankment of the Irrigation canal, downstream of the Gaoliang Sluice Gate, which is 600 m away from the gate. In addition to pumping water to the north, it can also be combined with storm drainage from the north canal of the west area. It was started in November 1983 and completed in September 1986. The Huaiyin Station was designed by the Provincial Water Conservancy Survey and Design Institute and was built by the Provincial Water Conservancy Construction Engineering Company. More details are given in Table 1.

3. Mathematical Model

3.1. One-Dimensional Open Channel Flow Equations

Unsteady flow problems in open channels are described by the Saint-Venant equations or their approximations, which are kinematic waves, non-inertia waves, gravity waves, and quasi-steady dynamic waves. These approximations are defined depending on the relative importance of the local inertia, pressure gradient, and gravity and friction effects involved in the physical mechanisms. One-dimensional, unsteady, gradually-varied open channel flow in prismatic channels is expressed by the Saint-Venant equations as follows [27]:

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=q,$$

(1)

$$k_l\frac{\partial Q}{\partial t}+k_c\left[\frac{\partial \left(Q^2/A\right)}{\partial x}-Vq\right]+k_{p}gA\frac{\partial y}{\partial x}-k_{f}gA\left(S_0-S_f\right)=0,$$

(2)

where Q represents the flow discharge (m³ s⁻¹); y is the water surface elevation (absolute water level above the sea level) (m); A is the flow cross-sectional area, (m²); q is the lateral outflow per unit length along the channel (m² s⁻¹); V is the mean flow velocity (m s⁻¹); g is the acceleration due to gravity (m s⁻²); $S_0$ is the channel bed slope; $S_f$ is the friction slope; x is the distance in the longitudinal direction (m); and t represents time (s). In Equation (2), the terms $k_l$, $k_c$, $k_p$, and $k_f$ represent the local acceleration, convective acceleration, pressure force, and gravity and friction forces, respectively. The values of the coefficients $k_l$, $k_c$, $k_p$, and $k_f$ are either one or zero, depending on the wave approximation considered. Yen and Tsai [27] determined the coefficient values as shown in

Table 2. Model coefficients for different wave approximations.

Wave Approximations	Model Coefficients
	${\mathit{k}}_{\mathit{l}}$	${\mathit{k}}_{\mathit{c}}$	${\mathit{k}}_{\mathit{p}}$	${\mathit{k}}_{\mathit{f}}$
Dynamic waves	1	1	1	1
Quasi-steady dynamic waves	0	1	1	1
Non-inertia waves	0	0	1	1
Kinematic waves	0	0	0	1

.

3.2. Hypotheses of the Saint-Venant Equations

The following assumptions are classically made when deriving the Saint-Venant equations [28],

• The flow is one-dimensional, i.e., the velocity is uniform over the cross-sectional area, and the water level across the section is horizontal.
• The streamline curvature is small and vertical accelerations are negligible; hence, the pressure is hydrostatic.
• The effects of boundary friction and turbulence can be accounted for through resistance laws analogous to those used for steady-state flow.
• The average channel bed slope is small so that the cosine of the angle it makes with the horizontal may be replaced by unity.
• The variation of the channel width along x is small.

The Saint-Venant equations are valid as long as the hypotheses stated are fulfilled.

3.3. Boundary Conditions

The discharges of the two pumping stations are derived from head-capacity (H-Q) curves. An H-Q curve can be expressed as a parabola. Therefore, mathematically, the H-Q curves can be represented by second-degree polynomial equations, as follows [29]:

$$H_1=a_0+a_1{Q}_1+a_2{Q}_1{}^2,$$

(3)

$$H_2=b_0+b_1{Q}_2+b_2{Q}_2{}^2,$$

(4)

where $Q_1$ and $Q_2$ are the discharges of step 1 and step 2 pumping stations, respectively, (m³ s⁻¹); $a_0$, $a_1$, $a_2$ and $b_0$, $b_1$, $b_2$ are constants based on step 1 and step 2 pumping stations, respectively. The constants $a_0$ and $b_0$, $a_1$ and $b_1$, and $a_2$ and $b_2$ depend on the pump geometry, impeller diameter, and pump speed, respectively. $Q_1$ and $Q_2$ are expressed as follows:

$$Q_1=\frac{-a_1\pm \sqrt{a_1^2-4a_2\left(a_0-y_1+h_1\right)}}{2a_2},$$

(5)

$$Q_2=\frac{-b_1\pm \sqrt{b_1^2-4b_2\left(b_0+y_2-h_2\right)}}{2b_2}.$$

(6)

At the free surface and along the wetted surface of the channel, a considerable part of the flow is lost via seepage and evaporation from the channel. The evaporation loss is proportional to the area of free surface. This loss depends on the supply of energy that can provide the latent heat of vaporization, and the ability to transport the vapor from the evaporating surface, which, in turn, depends on the wind velocity over the surface and the specific humidity gradient in the air above the water surface. The evaporation loss can be calculated by

$$q_v=T{E}_v,$$

(7)

where T is the top width of the channel’s flow (m); and $E_v$ is the evaporation discharge per unit free surface area (m s⁻¹). $E_v$ is given by $E_v=\left(e_s-e_d\right)f_w$, where $e_s$ is the saturation vapor pressure of the air at the temperature of the water surface (Pa); $e_d$ is the saturation vapor pressure of the air at the dew point (Pa); and $f_w$ is the wind function (m(s Pa)⁻¹). The difference between the saturation vapor pressure of the air at the temperature of the water surface and at the dew point $\left(e_s-e_d\right)$ was given by Cuenca [30] as

$$e_s-e_d=610.78\left[\exp \left(\frac{17.27{\theta }_w}{237.3+{\theta }_w}\right)-R_h\exp \left(\frac{17.27{\theta }_a}{237.3+{\theta }_a}\right)\right],$$

(8)

where ${\theta }_w$ is the water surface temperature (°C); ${\theta }_a$ is the mean air temperature (°C); and $R_h$ is the relative humidity expressed as a fraction. The wind function for a flowing channel was given by Fulford and Sturm [31] as

$$f_w=3.704\times {10}^{-11}\left(1+0.25u_2\right),$$

(9)

where $u_2$ is the wind velocity at 2 m above the free surface (m s⁻¹).

The seepage loss depends on the water in the channel, the channel bed-soil characteristics and the ground water level; for the estimation of seepage in a channel, Davis and Wilson [32] suggested the following formula:

$$q_s={10}^{-3}\times rS{y}^{1/3}P\left(\frac{1}{1+9.125\times {10}^{-4}\sqrt{V}}\right),$$

(10)

where $r=\min \left(1;0445\sqrt{z}{Q}^{-0.185}\right)$. z is the altitude difference between the water surface and the underground water (m); S is a soil parameter that depends on the soil type; and P represents the wetted perimeter of the cross section of the channel (m).

3.4. Initial Conditions

At the initial time point, the water level along the channel has the same elevation, which is known, and it is taken as initial condition, i.e.,

$$y=k,$$

(11)

where k is a constant. The flow velocity along the channel is calculated using the following equation:

$$V=Q/A.$$

(12)

Initially, the water has no velocity: $V=0$.

4. Finite Difference Method

Equations (1) and (2) are nonlinear, first-order, partial differential equations of hyperbolic type. Since these equations cannot be solved analytically, the finite difference method will be used to obtain approximate solutions. In this technique, the partial derivatives in the equations are replaced by their corresponding finite difference approximations, resulting in a series of algebraic equations that can be solved on a computer. The solution to the discrete problem represents an approximation of the solution to the continuous problem.

4.1. Lax Diffusive Scheme

The Lax diffusive scheme [4] is an explicit scheme that is first-order accurate in time and second-order accurate in space. This scheme is one of the simplest to program and yields satisfactory results for typical engineering applications [33], and it can efficiently solve the Saint-Venant equations [26]. This scheme divides the $x-t$ plane into a grid in which the grid interval along the $x$-axis is $\Delta x$ and the grid interval along the t-axis is $\Delta t$. For the distance axis, use $i$ for the $i\Delta x$ grid point and $i+1$ for the $\left(i+1\right)\Delta x$ grid point. For the time axis, use $j$ for the $j\Delta t$ grid point and $j+1$ for the $\left(j+1\right)\Delta t$ grid point. To refer to different variables at these grid points, use the number of the distance grid as a subscript and that of the time grid as a superscript. This scheme considers that the flow variables are known at time level $j$ and their values will determine the time level, $j+1$. Thus, the partial derivatives and the other variables for dynamic wave ($k_l=1$, $k_c=1$, $k_p=1$, $k_f=1$) are approximated as follows:

$$\frac{\partial A}{\partial t}=\frac{A_i^{j+1}\left(\left(A_{i+1}^j+A_{i-1}^j\right)/2\right)}{\Delta t},$$

(13)

$$\frac{\partial Q}{\partial t}=\frac{Q_i^{j+1}-\left(\left(Q_{i+1}^j+Q_{i-1}^j\right)/2\right)}{\Delta t},$$

(14)

$$\frac{\partial Q}{\partial x}=\frac{Q_{i+1}^j-Q_{i-1}^j}{2\Delta x},$$

(15)

$$\frac{\partial \left(Q^2/A\right)}{\partial x}=\frac{\left[{\left(Q_{i+1}^j\right)}^2/A_{i+1}^j\right]-\left[{\left(Q_{i-1}^j\right)}^2/A_{i-1}^j\right]}{2\Delta x},$$

(16)

$$\frac{\partial y}{\partial x}=\frac{y_{i+1}^j-y_{i-1}^j}{2\Delta x},$$

(17)

$$A=\frac{A_{i-1}^j+A_{i+1}^j}{2},$$

(18)

$$V=\frac{Q_{i-1}^j+Q_{i+1}^j}{A_{i-1}^j+A_{i+1}^j},$$

(19)

$$S_f=\frac{{\left(s_f\right)}_{i-1}^j+{\left(s_f\right)}_{i+1}^j}{2}.$$

(20)

In this work, the lateral outflow per unit length along the channel is considered in the evaporation loss and the seepage loss as $q=-(q_v+q_s)$,

$$\begin{array}{l}q=-q_v-{10}^{-3}\times S\times \left(\min \left(1;0.445\sqrt{z}{\left(\frac{Q_{i-1}^j+Q_{i+1}^j}{2}\right)}^{-0.185}\right)\right)\\ \times {\left(\frac{y_{i-1}^j+y_{i+1}^j}{2}\right)}^{1/3}\times \left(\frac{P_{i-1}^j+P_{i+1}^j}{2}\right)\times {\left(1+9.125\times {10}^{-4}\sqrt{\frac{Q_{i-1}^j+Q_{i+1}^j}{A_{i-1}^j+A_{i+1}^j}}\right)}^{-1}\end{array}.$$

(21)

Substituting Equations (13)–(21) into Equations (1) and (2) and simplifying, we obtain

$$\begin{array}{l}{A}_i^{j+1}=\frac{\left(A_{i+1}^j+A_{i-1}^j\right)}{2}-\frac{\Delta t\left(Q_{i+1}^j+Q_{i-1}^j\right)}{2\Delta x}-\Delta t{q}_v-\Delta t\times {10}^{-3}\times S\\ \times \left(\min \left(1;0.445\sqrt{z}{\left(\frac{Q_{i-1}^j+Q_{i+1}^j}{2}\right)}^{-0.185}\right)\right)\times {\left(\frac{y_{i-1}^j+y_{i+1}^j}{2}\right)}^{1/3}\\ \times \left(\frac{P_{i-1}^j+P_{i+1}^j}{2}\right)\times {\left(1+9.125\times {10}^{-4}\sqrt{\frac{Q_{i-1}^j+Q_{i+1}^j}{A_{i-1}^j+A_{i+1}^j}}\right)}^{-1}\end{array}$$

(22)

and

$$\begin{array}{l}{Q}_i^{j+1}=\frac{\left(Q_{i+1}^j+Q_{i-1}^j\right)}{2}-\Delta t\frac{\left[{\left(Q_{i+1}^j\right)}^2/A_{i+1}^j\right]-\left[{\left(Q_{i-1}^j\right)}^2/A_{i-1}^j\right]}{2\Delta x}-g\Delta t\frac{\left(A_{i-1}^j+A_{i+1}^j\right)}{2}\left(\frac{y_{i+1}^j+y_{i-1}^j}{2\Delta x}-S_0\right)\\ -g\Delta t\frac{A_{i-1}^j{\left(s_f\right)}_{i-1}^j+A_{i+1}^j{\left(s_f\right)}_{i+1}^j}{2}-\Delta t\left(\frac{Q_{i-1}^j+Q_{i+1}^j}{A_{i-1}^j+A_{i+1}^j}\right)q_v-\Delta t\left(\frac{Q_{i-1}^j+Q_{i+1}^j}{A_{i-1}^j+A_{i+1}^j}\right)\times {10}^{-3}\times S\\ \times \left(\min \left(1;0.445\sqrt{z}{\left(\frac{Q_{i-1}^j+Q_{i+1}^j}{2}\right)}^{-0.185}\right)\right)\times {\left(\frac{y_{i-1}^j+y_{i+1}^j}{2}\right)}^{1/3}\times \left(\frac{P_{i-1}^j+P_{i+1}^j}{2}\right)\times {\left(1+9.125\times {10}^{-4}\sqrt{\frac{Q_{i-1}^j+Q_{i+1}^j}{A_{i-1}^j+A_{i+1}^j}}\right)}^{-1}\end{array}.$$

(23)

The values of $A_0^t$, and $A_N^t$ need to be imposed in this scheme; we consider them as

$$A_0^t=A_1^t,$$

(24)

$$A_N^t=A_{N-1}^t,$$

(25)

where $N=\frac{L}{\Delta x}$; L represents the channel length (m).

Obviously, once $A_i^{j+1}$ has been evaluated, $y_i^{j+1}$ and $P_i^{j+1}$ can be determined from the cross-sectional geometry. In this paper, C# programming language is applied for the calculations, which is intended to be simple, modern, and provide support for software engineering principles.

4.2. Stability Conditions

The stability of the Lax diffusive scheme is provided under Courant conditions [34],

$$\Delta t\le \frac{\Delta x}{\left|V\mp c\right|}$$

(26)

where c is the water wave celerity (m s⁻¹), which is evaluated in terms of the hydraulic depth as

$$c=\sqrt{gD},$$

(27)

where D is the hydraulic depth (m) and $D=A/T$. Equation (26) is satisfied for all nodes at all of the time stages.

5. Case Study

To study the unsteady flow behavior in an open channel between two pumping stations, we considered a prismatic trapezoidal channel bound by Huai’an and Huaiyin pumping stations in Jiangsu Province, China. The channel length is L = 28.47 km, and the bottom width is b = 50 m. The channel bed and banks consist of sand and clay; the bed surface is flat, and the bank’s shape is uniform. The channel soil is clay, while the side slope can be represented by M = 1:3; and the longitudinal bed slope of the channel is equal to 1:10,000. The friction slope can be calculated from Manning’s equation as

$$S_f=\frac{n^2{V}^2}{R^{4/3}},$$

(28)

where n is the Manning roughness factor—the value of n depends mainly upon the surface roughness, the amount of vegetation, and the channel irregularity, and, to a lesser degree, upon stage, scour and deposition and channel alignment [4]. In our case study, we used the value of n computed by Wu et al. [35], which is n = 0.0225; In addition, R represents the hydraulic radius (m), which is defined by the following equation: $R=A/P$, where $P=b+2y\sqrt{1+M^2}$.

Initially, the water level of Huai’an pumping station (downstream) was $h_1$ = 5.6 m, and the water level of Huaiyin pumping station (upstream) was $h_2$ = 13.6 m. The Huai’an pumping station consists of two pump units, and the Huaiyin pumping station consists of four pump units. The water level values, 8.5 m, 9.13 m and 10 m, were taken as the initial elevation values. For different blade angles, the measured data of the head and discharge for one pump unit of the Huai’an and Huaiyin pumping stations were obtained from the pumping stations’ data. Figure 4 and Figure 5 show the H-Q curves of the two pumping stations.

The pumping station operation mode was selected as shown in

Table 3. The pumping station operation mode.

Name of Pumping Station	Number of Running Pump Units	Blade Angle
Huai’an pumping station	2	−6°
Huaiyin pumping station	4	4°

. Then, the least squares method [29] was applied to estimate the values of coefficients $a_0$, $a_1$, $a_2$, $b_0$, $b_1$, and $b_2$ in Equations (3) and (4), which uses H-Q curves for estimation.

Thus, the equivalent equations of the pumping station operation mode are

$$H_1=16.9148-0.1425Q_1+2.05429\times {10}^{-4}{Q}_1{}^2,$$

(29)

$$H_2=3.2917+0.1257Q_2-1.03529\times {10}^{-3}{Q}_2{}^2,$$

(30)

$$Q_1=\frac{0.1425-\sqrt{0.000821716y_1+0.0018055}}{0.000410858}$$

(31)

$$Q_2=\frac{0.1257+\sqrt{0.00414116y_2-0.0268878}}{0.00207058}.$$

(32)

The measured data of the evaporation from the channel can be represented by $q_v$ = 4.777 × 10⁻⁶ m² s⁻¹, the soil parameter is S = 0.0034, and the average value of z is 1.1 m. By substituting the values of the factors into Equation (21), the lateral out flow per unit length is

$$q=-4.777\times {10}^{-6}-\min \left(1;0.467Q^{-0.185}\right)\left(\frac{y^{1/3}\left(170+21.5035y\right)\times {10}^{-6}}{1+9.125\times {10}^{-4}\sqrt{\frac{Q}{50y+3y^2}}}\right).$$

(33)

6. Results and Discussion

In this work, the Lax diffusive scheme was used to calculate the water level, the flow rate and the lateral outflow per unit length along the channel at different times. All results were computed using a Microsoft Windows 7.0 computer with a 64-bit operating system, and an Intel(R) Xeon(R) CPU E5-2690 v4 processor at a speed of 2.60 GHz and with a memory of 32.0 GB.

6.1. Grid Size, Accuracy and Computational Time

For an initial water level value of 9.13 m, Figure 6 shows the results of water surface elevations for different grid sizes, i.e., different time and distance steps. We can observe that distance and time steps of the grid points determine the accuracy of the simulation. More grid points correspond to smaller step sizes, thereby leading to smaller errors and a more accurate solution. However, increasing the number of grid points requires more computational time. Therefore, the main challenge is to choose a sufficient grid size that can efficiently solve the model. This choice can be done by increasing the number of grid points until the change or the error in the solution is either negligible or sufficiently small enough to be acceptable within a reasonable computational time limit.

Figure 6a shows water surface profiles after 20 min (t = 20 min), while Figure 6b shows them after three and six days (t = 3 days, t = 6 days) for different grid sizes. From Figure 6a we can observe that the water surface profile errors were relatively large compared with the profile errors in Figure 6b. After 20 min, the maximum errors along the channel relative to the most accurate results in terms of grid size (Δx = 0.5 m, Δt = 0.005 s) were also calculated; these were 47.535 mm and 28 mm for grid sizes (Δx = 100 m, Δt = 1 s) and (Δx = 10 m, Δt = 0.1 s), respectively, while the error was only 1.367 mm for grid size (Δx = 1 m, Δt = 0.01 s). From Figure 6, it can be seen that the water surface profiles of grid sizes (Δx = 1 m, Δt = 0.01 s) and (Δx = 0.5 m, Δt = 0.005 s) overlapped, while the computational time required for grid size (Δx = 1 m, Δt = 0.01 s) was much less than the computational time required for the grid (Δx = 0.5 m, Δt = 0.005 s) (see

Table 4. Computational time for water surface elevations for different grid sizes and their maximum errors at t = 20 min, t = 3 days, and t = 6 days.

Time	Computational Time and Maximum Error	Distance and Time Step
		Δx = 100 m Δt = 1 s	Δx = 10 m Δt = 0.1 s	Δx = 1 m Δt = 0.01 s	Δx = 0.5 m Δt = 0.005 s
20 min	Computational time	735 ms	23 s	29 min 6 s	1 h 37 s
	Maximum error(mm)	47.535	28.000	1.367	—
3 days	Computational time	2 min 8 s	1 h 14 min 43 s	3 days 35 min 56 s	12 days 12 h 39 min 43 s
	Maximum error(mm)	0.396	0.105	0.013	—
6 days	Computational time	4 min 26 s	2 h 29 min 20 s	6 days 2 h 59 min 51 s	24 days 23 h 43 min 7 s
	Maximum error(mm)	0.242	0.033	0.004	—

). Therefore, in order to guarantee stability and accuracy, and limit the computational time, we fixed the grid size at (Δx = 1 m, Δt = 0.01 s). This will be used in the rest of the work. Figure 6b shows water surface profiles after three and six days (t = 3 days, t = 6 days). We observed that the water surface profile errors were small enough for the grid size (Δx = 100 m, Δt = 1 s). Table 4 presents the computational time required for computing water surface profiles for different grid sizes, and its maximum error along the channel relative to the most accurate results for a grid size of (Δx = 0.5 m, Δt = 0.005 s).

6.2. Water Surface Elevation Along the Channel

Figure 7 represents the water surface profiles at different time points for different initial water level values: 8.5 m, 9.13 m, and 10 m. Figure 7a,c,e reveals that the water surface profile changed rapidly and violently with time over the first two hours after start-up of the pump units. In this period of time, Huai’an pumping station was pumping water into the channel on one side, while Huaiyin pumping station was pumping water from the channel on the other side. Thus, the water surface profile rises at Huai’an pumping station and decreases at Huaiyin pumping station.

From day 1 until day 6, as shown in Figure 7b, the whole water surface profile increased with time until it became almost constant, while it decreased with time in Figure 7d,f until it became almost constant as a result of the water surface profile attempting to reach the last steady state situation. In addition, the pumping station discharges were influenced by the varying water level values. For the initial water level value of 8.5 m, the water levels at Huai’an pumping station (upstream) and Huaiyin pumping station (downstream) increased with time, whereas they decreased with time for the initial water level values of 9.13 m, and 10 m. The variation in the water surface profiles along the channel could be given as a function of distance, such as $y=y\left(x\right)$.

For different initial water level values (8.5 m, 9.13 m, and 10 m), the water level values were evaluated for the Huai’an pumping station (upstream) and the Huaiyin pumping station (downstream) at different times.

Table 5. Water level values of Huai’an pumping station (upstream) and Huaiyin pumping station (downstream) at different times for different initial water level values (m).

Time	Initial Water Level Value
	8.5 m		9.13 m		10 m
	Huai’an Pumping Station (Upstream) Water Level	Huaiyin Pumping Station (Downstream) Water Level	Huai’an Pumping Station (Upstream) Water Level	Huaiyin Pumping Station (Downstream) Water Level	Huai’an Pumping Station (Upstream) Water Level	Huaiyin Pumping Station (Downstream) Water Level
8 h	8.6375	8.5739	9.1465	9.0877	9.8714	9.8181
1 day	8.7926	8.7299	9.1288	9.0736	9.6303	9.5845
2 days	8.9299	8.8704	9.1141	9.0587	9.4006	9.3508
3 days	9.0043	8.9465	9.1059	9.0503	9.2677	9.2154
4 days	9.0452	8.9882	9.1014	9.0457	9.1921	9.1383
5 days	9.0677	9.0113	9.0989	9.0431	9.1495	9.0948
6 days	9.0802	9.0240	9.0975	9.0417	9.1257	9.0705
7 days	9.0871	9.0310	9.0967	9.0409	9.1124	9.0569
8 days	9.0909	9.0350			9.1050	9.0493
9 days	9.0930	9.0371			9.1009	9.0451
10 days	9.0942	9.0383			9.0986	9.0428
11 days	9.0949	9.0390			9.0973	9.0415
12 days					9.0966	9.0408
t days (error < 1 mm)	9.0949	9.0390	9.0967	9.0409	9.0966	9.0408
8 h − t days	−0.4574	−0.4651	0.0498	0.0468	0.7749	0.7773

shows that the water level values at the two pumping stations at the eighth hour had great differences compared with the last values; however, they converged to give the same results. The time periods required to reach a stable water surface elevation with an error of less than 1 mm for the initial water surface elevation values (8.5 m, 9.13 m, and 10 m) were 11 days, 7 days, and 12 days, respectively. In addition, regardless of whether the water surface elevation increased or decreased, the water surface elevation difference compared with the last value along the channel at different times was almost constant after one day.

Table 6. The times required to reach different water surface elevation error values along the channel over two consecutive days.

Water Surface Elevation Error Values	Initial Water Level Value
	8.5 m	9.13 m	10 m
5 mm	8 days	4 days	9 days
1 mm	11 days	7 days	12 days

shows the times to reach water surface elevation with errors of 5 mm and 1 mm along the channel over two consecutive days. When the initial water level values were closer to the last value, less time was required to converge.

From the calculated results for the different initial water level values, Huai’an pumping station’s last water level value was 9.096 m, and Huaiyin pumping station’s last water level value was 9.040 m; these values were compared with the measured data [36], as shown in

Table 7. A comparison between measured and computed water level results.

Items	Water Level of Huai’an Pumping Station (Upstream) (m)	Water Level of Huaiyin Pumping Station (Downstream) (m)	Lateral Outflow (m3 s−1)
Measured	9.100	9.000	2.210
Computed	9.096	9.040	2.164
Comparisons	−0.004	0.04	−0.046
% Error	0.044	0.444	2.081

. From the comparison we show that our model results are in accordance with field test measured results.

For the initial water level values (8.5 m, 9.13 m, and 10 m), the average water level values along the channel were computed. Figure 8 shows the differences between these average values and the last value at different times. Over time, the difference decreases until the steady state is reached.

6.3. Flow Rate Along the Channel

Figure 9 shows the flow rate variation with time along the channel for different initial water level values. In the first 2 h after start-up of the pumping station units, the flow rate changed violently as a result of the kinetic energy being imparting continuously to fluid that had no velocity. Subsequently, all of the water was moving at various velocities along the channel. From the computational results, the last discharge and head for the Huai’an pumping station were 112.4 m³ s⁻¹ and 3.496 m, respectively, while those for Huaiyin pumping station were 110.3 m³ s⁻¹ and 4.560 m, respectively. In the range of 40 to 50 min after start-up of the pump units, the flow rate in the channel reached its maximum value, i.e., for the initial water level values of 8.5 m, 9.13 m, and 10 m, the maximum flow rates were 190.218 m³ s⁻¹, 194.081 m³ s⁻¹, and 197.627 m³ s⁻¹, respectively. These flow rates were 1.71~1.78 times the speed of the last flow rate. The higher initial water values had higher maximum flow rates during the transient period. Measurements at the cross section of the channel of the pumping station flow rate that are performed too soon after starting the pump units will cause substantial error. Over time, the velocity change becomes increasingly small until a situation is reached where Huai’an pumping station’s discharge satisfies the summation of Huaiyin pumping station’s discharge and the lateral outflow along the channel (which is comprised of seepage and evaporation losses). In addition, we can observe that the flow rate in the channel depends on the pumping stations’ discharges, which depend on the heads of the pumps or the two end water level values of the channel.

The discharge variation and the head variation over time for the two pumping stations for different initial water level values are shown in Figure 10. After start-up of the pumping station units, the p discharge changed as a result of variation in the water level values of the channel. Huai’an pumping station’s discharge decreased over time for the initial water level value of 8.5 m, because the upstream pumping station’s water level increased, and thus, the water head increased. In contrast, Huai’an pumping station’s discharge increased over time for the initial water level values of 9.13 m and 10 m, because the upstream pumping station’s water level decreased. Huaiyin pumping station’s discharge increased over time for an initial water level value of 8.5 m, because the downstream pumping station’s water level increased; while the discharge decreased over time for the initial water level values of 9.13 m and 10 m, because the downstream pumping station’s water level decreased.

After 8 h, as shown in

Table 8. Huai’an pumping station’s discharge and Huaiyin pumping station’s discharge after 8 h and 1 day for the different initial water level values (m³ s⁻¹).

Time	Initial Water Level Value
	8.5 m		9.13 m		10 m
	Huai’an Pumping Station Discharge	Huaiyin Pumping Station Discharge	Huai’an Pumping Station Discharge	Huaiyin Pumping Station Discharge	Huai’an Pumping Station Discharge	Huaiyin Pumping Station Discharge
8 h	117.229	105.487	111.893	110.720	104.496	117.333
1 day	115.589	107.138	112.077	110.584	106.931	115.306
t days (water level difference <1 mm)	112.435	110.242	112.401	110.274	112.403	110.272
8 h − t days	4.794	−4.755	−0.508	0.446	−7.907	7.061

, the magnitude of variation of Huai’an pumping station’s discharge was approximately equal to the magnitude of the variation of Huaiyin pumping station’s discharge, but with a different sign, i.e., $\Delta Q_{\mathrm{I}}\approx -\Delta Q_{\mathrm{II}}$, where, $\Delta Q_{\mathrm{I}}$ and $\Delta Q_{\mathrm{II}}$ are the variations in Huai’an pumping station’s discharge and Huaiyin pumping station discharge’s, respectively.

6.4. Lateral Outflow per Unit Length Along the Channel

Figure 11 shows the lateral outflow per unit length in regard to its variation over time along the channel for different initial water level values. The evaporation loss depended on the area of the free surface, and the seepage loss depended on the wetted perimeter, the water velocity and the underground water level. As shown in Figure 11a,c,e the lateral outflow per unit length was very high because the water had no velocity for 20 min after the start-up of the pumping station units. Later, it decreased with time when the water gained velocity. The lateral outflow per unit length decreased to 20.4% when the water was moving compared with when the water was at rest. As shown in Figure 11b the lateral outflow per unit length rose with time until it became almost constant as a result of the area of free surface and the increase in the wetted perimeter as the water surface elevation rose, whereas it decreases with time in Figure 11d,f until it became almost constant as a result of the area of free surface and the decrease in the wetted perimeter as the water surface elevation decreased. The calculated results for different initial water level values showed a lateral outflow along the channel of 2.164 m³ s⁻¹. This value was compared with the measured data which was 2.210 m³ s⁻¹, and there was no significant error in our result. Therefore, our model can effectively evaluate the lateral outflow along the channel.

6.5. Transient Flow Formation Mechanism (Wave Propagation Mechanism)

Any change of flow in an open channel causes a wave to be propagated from the point where the change starts. Therefore, when starting the step 1 pumping station and the step 2 pumping station at the same time, two waves are propagated (positive and high wave, and negative and low wave), as shown in Figure 7 and Figure 12. The positive wave moves in the flow direction, whereas the negative wave moves in the opposite direction. Over time, these two waves meet each other near the middle of the channel, thereby increasing the water surface slope. As a result, the flow rate rises in this area and generates two waves that travel in opposite directions to the channel ends. Afterwards, the waves reflect back to the middle of the channel to meet each other again, but in this process, the flow rate decreases gradually from the two ends to the middle of the channel. Then, the two waves travel in opposite directions to the channel ends. Subsequently, they reflect back to the middle of the channel to meet each other again, and the process repeats itself. The amplitude of these two waves when they meet each other becomes increasingly small until a steady state is reached. From the results, these two waves take approximately 30 min to meet each other for the first time, with an average speed of 7.908 m s⁻¹ and an amplitude of 0.14 m. In addition, we can see that for an initial water level value of 9.13 m, from 20 to 80 min after start-up of the two pumping stations, the flow rate was much bigger than the pumping station’s flow rate, and at 1 h, the flow rate in the channel rose to 186 m³ s⁻¹. The wave velocity and amplitude are affected by many factors, such as time, flow cross sectional area, distance along the channel, flow rate, channel bed slope and depth of water. One can observe that the higher the initial water value is, the faster the wave speed is, and vice versa.

7. Conclusions

In a two-step pumping station system, a transient flow in the open channel is immediately generated after starting the pump units. In this paper, the unsteady flow behavior in a prismatic trapezoidal open channel between two step pumping stations was studied. We constructed a mathematical model to simulate one-dimensional unsteady, gradually-varied, open channel flow in the channel. The soil leakage, and evaporation from the water surface along the channel were also involved. The Lax diffusive scheme was applied to solve the proposed model. The model was applied on a real open channel allocated between the Huai`an and Huayin pumping stations in China. The summary and limitations of this study are as follows.

(1)

The accuracy of the model was calculated by varying the scheme grid steps. In order to guarantee the stability and the computational accuracy at all times, the optimal distance step obtained was Δx = 1 m, while the optimal time step was Δt = 0.01 s.

(2)

We calculated and discussed the water surface elevation, the flow rate and the lateral outflow per unit length at various times along the channel. In Huai’an pumping station (upstream), the last water level value was 9.096 m, and the pumping station’s discharge was 112.4 m³ s⁻¹ with a head of 3.496 m. In contrast, in Huaiyin pumping station (downstream), the last water level value was 9.040 m, and the pumping station’s discharge was 110.3 m³ s⁻¹ with a head of 4.560 m. The lateral outflow along the channel was 2.164 m³ s⁻¹.

(3)

The computed results showed that the water surface profile, the flow rate and the lateral outflow along the channel are influenced by the initial water level value at which the unsteady flow occurs. Different initial water levels tend towards the same results in terms of water surface profile and flow rate. The time required to reach the almost steady state after the transient produced in the open channel also depends on the initial water level value.

(4)

The wave propagation mechanism is clearly presented and analyzed. Interestingly, two waves were generated. A positive wave moves in the flow direction, whereas a negative wave moves in the opposite direction. Over time, these two waves meet each other near the middle of the channel, thereby increasing the water surface slope. As a result, the flow rate rises in this area and generates two waves that travel in opposite directions to the channel ends. Afterwards, the waves reflect back to the middle of the channel to meet each other again, but in this process, the flow rate decreases gradually from the two ends to the middle of the channel. Then, the two waves travel in opposite directions to the channel ends. Subsequently, they reflect back to the middle of the channel to meet each other again, and the process repeats itself. The amplitudes of these two waves when they meet each other become increasingly small until the steady state is reached.

(5)

The accuracy of the results was assessed by comparing the computed results with measured data.

Our future work based on open channel transient flow will be as follows:

(1)

We will study the transient flow in an open channel between two step pumping stations under abnormal situations, such as when one pumping station cannot start up, or when one pump unit breaks down suddenly.

(2)

We will study the transient flow for a system with three and more step pumping stations.

(3)

We will study the optimal operation of a step pumping station system.