One of the main problems in water management of irrigation systems is the control of the equitable distribution of water among different orifice offtakes. The difficulty of managing a canal is partly caused by the lack of knowledge of the canal state because the scheduled demand is often not fulfilled, since farmers extract more water than is scheduled and it is impossible for the watermaster to determine the canal state. However, an innovative developed algorithm called CSE is proposed in this paper. This algorithm is able to estimate the real extracted flow and the hydrodynamic canal state (that is, the water level and velocity along the irrigation canal). The algorithm solves an inverse problem implemented as a nonlinear optimization problem using the Levenberg–Marquardt method. The algorithm is tested, taking into account several numerical examples, and a practical implementation is made for a real case study in the PAC-UPC canal, a 220 m laboratory canal especially designed for research into irrigation canal control area and irrigation canal modelling. This useful algorithm evaluates the real extraction flow and the canal state and could be a useful tool for a feedback controller.

INTRODUCTION

Frequently, irrigation managers are worried about the equitable distribution of water among different orifice offtakes. Some farmers may demand much more water than was scheduled, and this normally involves several problems such as important disadjustment on canal management. In such cases, we propose solving this problem with the CSE algorithm, which is able to estimate the real extracted flow and the hydrodynamic canal state (that is, the water level and velocity along the irrigation canal). The water saved may then be used to increase the irrigated area or save money. Alternatively, reduction of water consumption and peak requirement may allow a reduction in the water withdrawal from rivers or reservoirs, reduction in the capacity, and hence cost, of storage works, and of conveyance and distribution works such as canals and pumping stations. The potential sophistication of on-farm water management is highly dependent upon the level of water delivery service provided to individual farms, which in turn, depends upon the conveyance manageability within the complete water distribution system (Plusquellec 1988). The Food and Agriculture Organization of the United Nations indicates that for the year 2030 agricultural production will have to be increased by +80% to fulfil food demand, but it will have to be done without the possibility of increasing water withdrawals by more than +12% (Faurès et al. 2010). This can be done by reducing spillages along canals, which is technically possible since studies indicate that in the Mediterranean region, by the year 2025, water savings can be 65% for irrigation, 22% for industrial use and 13% for domestic use (Institut Méditerranéen de l'Eau 2007). These savings can be obtained through the implementation of several measures allowing better management of water demand at several institutional and technical levels. The survey carried out in the south of France in 1997 confirmed the figures found in other sources (ASCE 1993; Mareels et al. 2005), showing that water losses at the distribution level were around 50% on average and could be reduced to less than 10% by the modernization of canal real-time operation including automation (Malaterre & Rogers 2008). In addition, the investment in automation allows the ‘simplification of operational control, and makes possible to reduce and/or optimize the number and qualifications of the operating staff’ (Goussard 1993).

In canal operations, the scheduling of gate operations to satisfy known changes in water demands is a necessity established by some authors. All our effort is focused on controlling and managing the flow deliveries of the system because in that way the profitability of the system can be increased. One of the obstacles to scheduling the gate operations for an irrigation period is providing real input data to the control algorithms, because it is possible to predict the initial conditions in a canal from the water level measures and introduce the backwater profile in the algorithm, but the problem is due to unknown flows extracted by farmers during the canal operation. In that sense, the input data introduced in the control algorithm is far from the reality and the challenge is to calculate a scheduling of gate operations to satisfy the water demands.

The management of water resources in a canal is sometimes a difficult task, since the delivery scheduled the day before could be modified by the users during the irrigation cycle without notice. For instance, the demand deliveries could be increased because a farmer extracted a greater flow than the scheduled delivery, so the watermaster would not know the reason why the water level decreased at the target points. In that sense, the proposed algorithm (CSE) would be a useful tool for the watermaster. CSE is able to obtain first the unknown flows withdrawn by farmers and, additionally, the hydrodynamic canal state, and thus the water level and velocity at each cross-section of the canal. On the other hand, if a farmer pumps water from an unknown point in the canal, CSE can obtain the approximate location of this point and the volume of water withdrawn by the farmer, so this can be a useful tool in the management of any canal.

There are a few authors who have developed algorithms involved in calculating the offtake discharges in irrigation canals, including Delgoda et al. (2013) and Van Overloop et al. (2008).

The performance of this algorithm is quite simple. In a case where we introduce a flow change in a particular section of a canal, the perturbations associated with this flow change can be estimated from a computer model based on the Saint-Venant equations. CSE solves the inverse problem from the registered perturbations (water level measurements at checkpoints), and it calculates the flow changes associated with them.

The input variables for the CSE algorithm are the water level measurements at certain points of the canal and the gate trajectories during a past time horizon. The output variables of CSE are the extracted flow during this time horizon and the hydrodynamic canal state at the current time obtained by implementing an inverse problem.

To achieve this objective, the CSE algorithm uses the hydraulic influence matrix (HIM), which establishes a relationship between the real extraction flow and the water depth and velocity at any point in the canal during a past time horizon. This matrix was defined by several authors (Soler 2003; Soler et al. 2013; Bonet 2015). In order to make the meaning of the HIM more understandable, the next example is given.

When we pump in a section of the canal, we modify the canal state, water level and velocity at any time in particular points (Figure 1). For instance, at t = K + 1, we modify the water level and velocity in cross-sections close to the pump location (cross-section i), but we do not modify the hydrodynamic variables in a far cross-section (i + 1). Instead, at t = K + 2, we modify the hydrodynamic variables in all sections between i to i + 1. In this sense, the HIM considers the range of influence of a flow change at every section of the canal during a time interval (the theoretical definition of the HIM matrix is introduced in the next section).
Figure 1

Pumping in a cross-section i.

Figure 1

Pumping in a cross-section i.

First, the algorithm was tested in a canal with two pools introducing an unknown extraction flow. In a second example, the CSE was tested with the test cases (Clemmens et al. 1998) introduced by the ASCE Task Committee on canal automation algorithms. The last example was done in a laboratory canal, PAC-UPC, especially designed to develop basic and applied research in irrigation canal control area located in Barcelona's School of Civil Engineering.

METHODS

The CSE solves an inverse problem to obtain the flow disturbances from the water level changes at several points in the canal (1), usually next to the canal offtakes: 
formula
1
where ΔY represents the changes in water level measurements at selected points of the canal, ΔQb represents the flow disturbance, HIM′(Qb) is the simplified hydraulic influence matrix that represents the influence of an extraction flow on the water level at different points of the canal, and HIM(Qb) is the hydraulic influence matrix that represents the influence of an extraction flow on the water level and velocity at the canal.

The HIM matrix is a square matrix and positive definite matrix (see Bonet 2015). The method used to solve the nonlinear optimization problem is the Levenberg–Marquardt method, which is a robust method with easy implementation and is a special method to solve an ill-conditioned matrix such as the HIM matrix. Thus, the algorithm solves a nonlinear optimization problem using the Levenberg–Marquardt method to evaluate the last expression (1).

The HIM matrix

The HIM matrix is based on the full Saint-Venant equations, which describe the free surface flow in canals. In partial derivatives this system is of the hyperbolic, quasi-linear and second-order type. The two equations are based on mass and momentum conservation.

Like any hyperbolic system, it can be transformed into its characteristic form. Such transformation of the Saint-Venant equations gives an ordinary system of four Equation (2) (Walker & Skogerboe 1987): 
formula
2
where y is the water level measured from the bottom of the canal, v is the weighted average velocity in a cross-section, x is the space, t is the time, S0 is the canal bottom slope, Sf (y, v) is the friction slope and c is the celerity of a gravity wave, where A(y) is the area of the wetted surface or a cross-section, x+ is the position of the upstream characteristic curve, x is the position of the downstream characteristic curve and T(y) is the top width of the free surface.
The system (2) has no analytical solution, so the use of numerical techniques is necessary. There are many methods that can be used. In order to have the longest possible integration time period with a minimum loss of accuracy, a discretization with second-order finite differences has been adopted, known as ‘the method of characteristic curves’ (Crandall 1956; Strelkoff 1970; Ames 1977; Gómez 1988). If this method is applied to Equation (2) and the characteristic curves that contain the points P-R and Q-R (Figure 2), respectively, are taken into account, we obtain the following equations: 
formula
3
where SfR=Sf(yR,vR), SfP=Sf(yP, vP), SfQ=Sf(yQ, vQ) and 0 ≤ θ ≤ 1 is the weighting coefficient that indicates the type of numerical scheme used. When θ = 1, the numerical scheme is implicit, if θ = 0 it is explicit, and when θ = 1/2 the numerical scheme is in central differences, the method of the characteristic curves.
Figure 2

The steps for the interpolation onto a structured grid.

Figure 2

The steps for the interpolation onto a structured grid.

If the flow conditions at points P′ and Q′ are known, xP′, tP′, yP′, vP′ and xQ′, tQ′, yQ′, vQ′ are also known, so xR, tR, yR, vR remain as unknowns (Figure 2), which can be found by calculating the last four Equation (3), using any of the methods to solve nonlinear equations, such as the Newton–Raphson method. The way of calculating the influences shown in this section is closely linked to the numerical scheme of characteristic curves. However, usually this scheme is not exactly used because it gives the solution at a point R whose coordinates (xR, tR) are unknown a priori. These coordinates are part of the solution, and normally it is more important to know the solution of the flow conditions at a specific point and at specific time instants. To solve this problem, first interpolate and then solve.

A structured grid like this one (Figure 2) creates a new nomenclature. Indeed, every variable will have a double index, where k refers to time and i to space. Thus, yik and vik represent the values for water level and average velocity at the coordinates xi=iΔx and tk=kΔt where Δx and Δt are selected by the user.

Obviously, the same system (3) is solved, but now with the new unknowns xP, yR, vR and xQ, where the values of yp and yQ are dependent on the value of xP and xQ evaluated using an interpolation function of second order too, (to be coherent with the numerical procedure used) we have used the Lagrange factors (a way of representing quadratic splines). For a dummy variable z, the result is: 
formula
4
In this way the variables yP, vP, yQ and vQ become functions of xP and xQ, as follows: 
formula
5
On the other hand, there are many control structures in canals such as gates, orifice offtakes, lateral weirs, etc. which allow flow control according to the specification of the watermaster. The individual study of each one is impossible in this work, so for this reason the most common structures are introduced. A common one found is a checkpoint, a target point where the water level is measured with a depth gauge, and it includes a sluice-gate, a lateral weir outlet, offtake orifice or a pump, as can be seen in Figure 3. The interaction of this control structure with the flow can be described according to the mass and energy conservation Equation (6): 
formula
6
where:

S(ye) is the horizontal surface of the reception area in the checkpoint position

A(ye)*ve is the incoming flow to checkpoint, defined in terms of water level and velocity

A(ys)*vs is the outgoing flow to checkpoint which continues along the canal, described in terms of water level and velocity

Cd is the discharge coefficient of the sluice-gate and ac is the sluice-gate width

d is the checkpoint drop and u is the gate opening

qb is the pumping offtake

qs(ye) is the outgoing lateral flow through the weir where Cs is the discharge coefficient, as is the weir width and y0 is the weir height measured from the bottom, called weir equation

qofftake(ye) is the outgoing offtake orifice flow where Cd is the discharge coefficient, A0 is the area of the offtake orifice, called orifice offtake equation.

The presence of checkpoints or control structures in the canal leads to the sub-division into canal pools, in such a way that there is always a canal pool between two checkpoints, and there is a checkpoint between two pools. If we discretize the control structure in a structured grid (Figure 4), linked with the characteristics of Equation (3), and then change the nomenclature, we should rewrite the control structure Equation (6) arriving at the following system of six equations (Equation (7)).
Figure 3

Diagram of a checkpoint with gate, lateral weir and pump.

Figure 3

Diagram of a checkpoint with gate, lateral weir and pump.

Figure 4

Graph with discretization of the control structure equations.

Figure 4

Graph with discretization of the control structure equations.

Thus, represents the water level at node n in the section upstream of the control structure at time k+1, that is, the incoming water level ye. In the same way, is defined as the existing water level at the first node of the downstream pool from the checkpoint at the same time k+1, and ys the outgoing water level at the control structure (Figure 4). The same can be said for the velocities and : 
formula
7
where: 
formula
and xQ are the unknown variables.
In order to continue with the calculation of the influences of a particular variable on the canal flow, we define the extraction flow qb as our particular variable. Thus, applying the first derivative of the functions (7) with respect the extraction flow, we obtain: 
formula
8

In Equation (8), for the first time, the extraction flow qb explicitly appears in the description. Despite the fact that the specific form of this function is still unknown, Equation (8) shows that the influence of the parameter qb on flow conditions at time k+1 is the sum of the indirect influence of the conditions at instant k and the direct influence at instant k+1 through the term ‘L’, which represents the variation in the extraction flow.

As a result, the method of characteristics is applied to the Saint-Venant equations in order to obtain algebraic equations to establish a relation between the influence parameter qb and the hydrodynamic canal state, and all the influences are lumped together in a global matrix, which is referred to as HIM(Qb). Based on this system of equations, and using the first derivative (∂y/∂qb, ∂v/∂qb) on an analytical process, we can establish the changes in flow behaviour (water level and velocity) due to a flow change at a point at a certain time instant.

The optimization problem

The inverse problem (1) is formulated as an unconstrained optimization problem. The optimization problem applied is the classical nonlinear problem without constraints and the method used to solve it is the Levenberg–Marquardt method (Press et al. 1992). Some authors have used optimization problems to solve inverse problems, e.g., gate stroking (Wylie 1969) and CLIS (Liu et al. 1998), but all of them implement inverse problems in control algorithms.

To introduce the optimized problem, some vectors applied in the problem have to be evaluated. As explained before, the CSE algorithm needs as input data the water level measurements at some points (checkpoints). Now, let us consider a vector, which contains the water level measurements at the checkpoints at time instant k (9), whose dimension is nc, where nc is the number of checkpoints: 
formula
9
Finally, every vector (9) defined for each time instant of the past time horizon is combined to define the ‘measured water level vector’, whose dimension is ny, where ny = kF × nc, where kF is the final instant of the past time horizon. We define this vector as: 
formula
10
We can check the measured water level vector values in a computational grid in Figure 5.
Figure 5

Sketch of a numerical grid of a canal with two pools controlled by two checkpoints downstream of each pool. There are pump stations close to each checkpoint. Pump flow trajectories are defined with four operation periods. Also, it shows the x/t-dots where the flow behaviour is defined. Notice that ‘K’ with a capital letter denotes time interval of control and ‘k’ with a small letter denotes time instant of simulation.

Figure 5

Sketch of a numerical grid of a canal with two pools controlled by two checkpoints downstream of each pool. There are pump stations close to each checkpoint. Pump flow trajectories are defined with four operation periods. Also, it shows the x/t-dots where the flow behaviour is defined. Notice that ‘K’ with a capital letter denotes time interval of control and ‘k’ with a small letter denotes time instant of simulation.

In another way, we can obtain the ‘state vector’ x(k), which is defined as the vector containing the numerical solution at the time instant k of all the discretization points: 
formula
11
where yi(k) and vi(k) = water depth and mean velocity at point i; and ns = number of cross-sections in which the canal is discretized. In this way, the vector x(1) is the known initial condition.

The state vector at the current time defines the current hydrodynamic state. The state vector is shown in a computational grid in Figure 5 (triangles).

All state vectors (11) may be included for each k-instant during a past time horizon into a single vector that is called ‘prediction vector’ (12). The dimension of this vector is nx = (2 × ns) × kF: 
formula
12
We are only interested in the water level at target points where we also obtain the water level measurements. We define a new vector that contains the water depth values given at a prescribed number of points (nc) at the time instant k: 
formula
13
This vector is constituted by a subset of values of the state vector (11).
All vectors shown in Equation (13) for all the time instants during the past time horizon are lumped in the so-called ‘prediction output vector’: 
formula
14
The dimension of the prediction output vector is nY = kF × nc. The vector (14) contains all water depth values and is shown in Figure 5. If you look closely at this figure, the position of the elements of the vector (14) in the grid domain coincides with the elements of the measured water level vector.
The prediction output vector is clearly related to the prediction vector (12) in the form: 
formula
15
where C is a matrix, called a discrete observer matrix by Malaterre (1994), a matrix of dimension nY × nX and whose components are only ‘zeros’ or ‘ones’. This matrix defines the direction of the control logics along a canal pool: downstream level control, upstream level control or control of intermediate water levels.
As previously introduced, CSE calculates the extracted flow at several points (for instance, pump stations) during a past time horizon. In that case, as illustrated in Figure 6, the pump stations are operating with an operation period K. Then, the extracted flow trajectories can be approached with piecewise functions. The extracted flow vector is defined by lumping together all the extracted flows during the past time horizon, as follows: 
formula
16
where the dimension of this vector is nQ = nP × KF, nP is the number of pump stations and KF is the final operation period of the past time horizon.
Figure 6

Mathematical representation of a pump flow trajectory.

Figure 6

Mathematical representation of a pump flow trajectory.

In this way, only Qb determines canal behaviour along the past time horizon. When the extracted flow trajectories are implemented in the canal, the flow response given by the model will be unique. Inversely, one flow behaviour is caused by only one set of extracted flow vectors, as a flow change is also responsible for water level disturbances.

We can check the extracted flow vector values in a computational grid in Figure 5.

Once CSE has estimated the extracted flow vector, the algorithm can also estimate the state vector at the current time, considering all variables that control the flow in the canal such as the scheduled demands, the initial conditions, and the gate trajectories during the past time horizon. Thus, the CSE algorithm calculates the flow disturbance (ΔQb), which better explains the changes between the measured and predicted water levels (ΔY) (1).

If we focus on the optimization problem, the objective is to make the prediction output vector more similar to the measured water level vector by manipulating the extracted flow vector (see Gill et al. 1982; Fletcher 1987). In mathematical terms, the objective is to obtain the extracted flow vector that minimizes the following performance criterion: 
formula
17
where Q′ matrix is a weighing matrix and the dimension of the matrix is nY×nY. This matrix could be used to weight the water level error at a particular checkpoint. This matrix is defined as the identity matrix in CSE. Qb contains the extracted flow trajectories (16).

RESULTS AND DISCUSSION

Numerical example: a canal introducing a single disturbance

Several scenarios are proposed to test the CSE algorithm in a canal that has two pools separated by sluice-gates (Figure 7). The flow is controlled by a gate downstream from the reservoir. Water is delivered through gravity outlets at the downstream end of each pool, where the checkpoints are located. There are pumping stations at the end of each pool that can introduce disturbances in the system in space and time.
Figure 7

Canal profile of a canal with a single pool.

Figure 7

Canal profile of a canal with a single pool.

The canal, with a trapezoidal section, is represented in Figure 7, and the general data are shown in Table 1. The characteristics of the checkpoints, sluice-gates, pump stations and orifice offtakes are shown in Tables 2 and 3.

Table 1

General canal features

Pool number Pool length (km) Bottom slope (%) Side slopes (H:V) Manning's coefficient (n) Bottom width (m) Canal depth (m) 
2.5 0.1 1.5:1 0.025 2.5 
II 2.5 0.1 1.5:1 0.025 2.5 
Pool number Pool length (km) Bottom slope (%) Side slopes (H:V) Manning's coefficient (n) Bottom width (m) Canal depth (m) 
2.5 0.1 1.5:1 0.025 2.5 
II 2.5 0.1 1.5:1 0.025 2.5 
Table 2

Sluice-gate features (canal structures)

Number of control structure or checkpoint Gate discharge coefficient Gate width (m) Gate height (m) Step (m) 
0.61 5.0 2.5 0.6 
0.61 5.0 2.5 0.6 
Number of control structure or checkpoint Gate discharge coefficient Gate width (m) Gate height (m) Step (m) 
0.61 5.0 2.5 0.6 
0.61 5.0 2.5 0.6 
Table 3

Pump station/orifice offtake features (canal structures)

Number of control structure or checkpoint Discharge coef./diameter orifice offtake (m) Orifice offtake height (m) Lateral spillway height (m) Lateral spillway width (m)/discharge coefficient 
– – – – 
0.6/0.77 1.0 2.3 10/1.99 
Number of control structure or checkpoint Discharge coef./diameter orifice offtake (m) Orifice offtake height (m) Lateral spillway height (m) Lateral spillway width (m)/discharge coefficient 
– – – – 
0.6/0.77 1.0 2.3 10/1.99 

In these examples, an upstream large reservoir is considered, whose water level Hreservoir is 3 m constantly throughout the test. This is the upstream boundary condition. At the end of the last pool, there is a control structure with orifice offtake and a pump station. The flow through the orifice offtake depends on the upstream water level of the orifice, and the disturbance is introduced by the pump station. This is the downstream boundary condition. This example starts from an initial steady state (Tables 4 and 5), with a specific and constant demand delivery at the end of the pools (5 m3/s through the orifice offtake), and the disturbance is not introduced initially.

Table 4

Initial conditions in the canal

Control structure Initial flow rate (m3/s) Control structure Initial water level (m) 
Gate 1 10.0 Checkpoint 1 2.0 
Gate 2 5.0 Checkpoint 2 2.0 
Control structure Initial flow rate (m3/s) Control structure Initial water level (m) 
Gate 1 10.0 Checkpoint 1 2.0 
Gate 2 5.0 Checkpoint 2 2.0 
Table 5

Flow delivered by an orifice offtake at the initial time step

Control structure Flow delivered through an orifice offtake (m3/s) 
Gravity outlet 1 5.0 
Gravity outlet 2 5.0 
Control structure Flow delivered through an orifice offtake (m3/s) 
Gravity outlet 1 5.0 
Gravity outlet 2 5.0 

The disturbance

In order to test the algorithm, a disturbance in the canal is introduced, which is unknown for CSE. There will be differences between the measured and expected water levels after the disturbance. The water level measurements (Figure 8) are obtained with a flow disturbance of 2 m3/s (pump station 1 for 15 minutes, from minute 40 to 55). This disturbance is introduced to the numerical model, based on the model of characteristics introduced before, as a flow change, and in this way we get the water level values. Once the water levels are introduced in the CSE algorithm, it will propose the pump flow trajectories that describe with the best accuracy the variation of water level at the checkpoints during the past time horizon.
Figure 8

Water level measured at checkpoint 1 and 2 during the past time horizon in the first example.

Figure 8

Water level measured at checkpoint 1 and 2 during the past time horizon in the first example.

Flow disturbance reduces the water level at checkpoint 1 from 2 m to 1.60 m, and at checkpoint 2 from 2 m to 1.92 m (see Figure 8). In that sense, a flow change of 2 m3/s at pump station 1 has a significant impact on the canal profile, although the water level at checkpoint 1 and 2 recovered to the desired water level (2 m) in just 160 minutes and 150 minutes, respectively, due to the short period of time that the disturbance was applied. As soon as the water level at the checkpoints recovers to the desired water level at these points, the flow through the orifice offtakes returns to 5 m3/s.

Results

The results of the pump flow obtained by CSE are shown in Figures 9 and 10. The results of the real pump flow (disturbance) and the simulated pump flow calculated by CSE can be compared. Both curves match exactly, and differences between them are practically nonexistent. The maximum error is around 0.005 m3/s, so considering that the real pump flow is 2 m3/s, the percentage maximum error between real and estimated flow is 0.25%.
Figure 9

Pump flow at checkpoint 1 during the past time horizon.

Figure 9

Pump flow at checkpoint 1 during the past time horizon.

Figure 10

Pump flow at checkpoint 2 during the past time horizon.

Figure 10

Pump flow at checkpoint 2 during the past time horizon.

CSE computes the unscheduled offtake changes introduced to the canal and writes these offtake changes as a pump flow trajectory (temporal function) associated with a certain pump station. The algorithm assigns to the pump stations all flow changes between the initial steady state to the current state (Bonet 2015).

CSE is also able to get the hydrodynamic state of a canal during an interval of time and the current canal state, which is also very useful for a feedback controller that has to know the hydrodynamic canal state in real time to calculate the optimum gate trajectories. In that sense, we show the hydrodynamic canal state calculated by CSE in this example at 2,700 seconds (5 minutes after introducing the disturbance). We can check the water profile simulated by a model and the water profile simulated by CSE in Figure 11.
Figure 11

Real water profile obtained by a computer model V.S. Water profile obtained by CSE at 2,700 seconds.

Figure 11

Real water profile obtained by a computer model V.S. Water profile obtained by CSE at 2,700 seconds.

Both water profiles are similar, as the accuracy of CSE in calculating the extracted flow vector is so high. An error of 5 L/s in the extracted flow is equivalent to a water level error of 1 mm and a velocity error of 0.001 m/s in cross-sections close to the extracted point.

Numerical example: a canal with multiple disturbances at the same time

Test cases

In this numerical example, we want to demonstrate that CSE is able to obtain the real extraction flow in a canal with several pools with multiple flow extractions at the same time. In that case, a test case proposed by the ASCE Task Committee to evaluate control algorithms, which we had already used to test CSE, is introduced.

Two canals are considered by the ASCE Task Committee (Clemmens et al. 1998) in several scenarios, but here only one is considered, the Maricopa Stanfield canal. The canal has eight pools separated by undershot sluice-gates. All the canal pools have been discretized and numbered in the direction of flow from upstream to downstream. Geometric characteristics of canal 1 (Maricopa Stanfield canal) are shown in Table 6 and Figure 12, and the control structures in Table 7. In the Maricopa Stanfield canal, there are gravity outlet orifices at the downstream end of each canal pool.
Table 6

Features of Maricopa Stanfield canal pools

Pool number Pool length (km) Bottom slope Side slopes (H:V) Manning's coefficient (n) Bottom width (m) Canal depth (m) 
0.1 2 × 10−3 1.5:1 0.014 1.1 
II 1.2 2 × 10−3 1.5:1 0.014 1.1 
III 0.4 2 × 10−3 1.5:1 0.014 1.0 
IV 0.8 2 × 10−3 1.5:1 0.014 0.8 1.1 
2 × 10−3 1.5:1 0.014 0.8 1.1 
VI 1.7 2 × 10−3 1.5:1 0.014 0.8 1.0 
VII 1.6 2 × 10−3 1.5:1 0.014 0.6 1.0 
VIII 1.7 2 × 10−3 1.5:1 0.014 0.6 1.0 
Pool number Pool length (km) Bottom slope Side slopes (H:V) Manning's coefficient (n) Bottom width (m) Canal depth (m) 
0.1 2 × 10−3 1.5:1 0.014 1.1 
II 1.2 2 × 10−3 1.5:1 0.014 1.1 
III 0.4 2 × 10−3 1.5:1 0.014 1.0 
IV 0.8 2 × 10−3 1.5:1 0.014 0.8 1.1 
2 × 10−3 1.5:1 0.014 0.8 1.1 
VI 1.7 2 × 10−3 1.5:1 0.014 0.8 1.0 
VII 1.6 2 × 10−3 1.5:1 0.014 0.6 1.0 
VIII 1.7 2 × 10−3 1.5:1 0.014 0.6 1.0 
Table 7

Maricopa Stanfield control structures

Target points Gate discharge coefficient Gate width (m) Gate height (m) Step (m) Length from gate 1 (km) Orifice offtake height (m) 
0.61 1.5 1.0 1.0 – 
0.61 1.5 1.1 1.0 0.1 0.45 
0.61 1.5 1.1 1.0 1.3 0.45 
0.61 1.5 1.0 1.0 1.7 0.40 
0.61 1.2 1.1 1.0 2.5 0.45 
0.61 1.2 1.1 1.0 4.5 0.45 
0.61 1.2 1.0 1.0 6.2 0.40 
0.61 1.0 1.0 1.0 7.8 0.40 
– – – – 9.5 0.40 
Target points Gate discharge coefficient Gate width (m) Gate height (m) Step (m) Length from gate 1 (km) Orifice offtake height (m) 
0.61 1.5 1.0 1.0 – 
0.61 1.5 1.1 1.0 0.1 0.45 
0.61 1.5 1.1 1.0 1.3 0.45 
0.61 1.5 1.0 1.0 1.7 0.40 
0.61 1.2 1.1 1.0 2.5 0.45 
0.61 1.2 1.1 1.0 4.5 0.45 
0.61 1.2 1.0 1.0 6.2 0.40 
0.61 1.0 1.0 1.0 7.8 0.40 
– – – – 9.5 0.40 
Figure 12

Maricopa Stanfield profile. The arrows mark the position of checkpoints. The first pool is number I and the first checkpoint is number 1.

Figure 12

Maricopa Stanfield profile. The arrows mark the position of checkpoints. The first pool is number I and the first checkpoint is number 1.

Scenario: test case 1-2 (Maricopa Stanfield)

After the first 2 hours, the initial backwater profile is changing, so an unscheduled flow change is introduced in the canal modifying some offtake flow, as shown in Table 8.

Table 8

Initial and unscheduled offtake changes on test case 1-2

Pool number Offtake initial flow (m3/s) Check initial flow (m3/s) Unscheduled offtake changes at 2 hours (m3/s) Check final flow (m3/s) 
Heading – 2.0 – 2.0 
0.2 1.8 – 1.8 
0.0 1.8 0.2 1.6 
0.4 1.4 −0.2 1.4 
0.0 1.4 0.2 1.2 
0.0 1.4 0.2 1.0 
0.3 1.1 −0.1 0.8 
0.2 0.9 – 0.6 
0.9 0.0 −0.3 0.0 
Pool number Offtake initial flow (m3/s) Check initial flow (m3/s) Unscheduled offtake changes at 2 hours (m3/s) Check final flow (m3/s) 
Heading – 2.0 – 2.0 
0.2 1.8 – 1.8 
0.0 1.8 0.2 1.6 
0.4 1.4 −0.2 1.4 
0.0 1.4 0.2 1.2 
0.0 1.4 0.2 1.0 
0.3 1.1 −0.1 0.8 
0.2 0.9 – 0.6 
0.9 0.0 −0.3 0.0 

The unscheduled deliveries are more significant at target 8 in test case 1-2, where the flow rate changes from 0.9 m3/s to 0.6 m3/s (33%).

Results

We show the results obtained by CSE divided by eight graphs, one for every pool.

Every graph shows the scheduled and unscheduled offtake changes (demanded extracted flow) by the farmers, the real value delivered by the offtake and the real extracted flow calculated by CSE.

In this test 1-2, the canal is in steady state during the first 2 hours. After the first 2 hours, unscheduled water deliveries are introduced to the system at 7,200 s, although the algorithm takes no notice until the next regulated period, once the water level is measured at the checkpoints. It is important to note that the unscheduled water deliveries are relevant in all targets but especially at target 8, because in just one regulation period, the water level at checkpoint 8 increases from 0.8 m to 1.08 m, due to the introduction of a water delivery change of 0.3 m3/s at checkpoint 8.

The extracted flow is calculated by CSE, which establishes the result with a maximum numerical error of 2% between the real extracted flow and the value by CSE, as shown in Figure 13.
Figure 13

Demanded extracted flow in pools 1–8: extracted flow by CSE and real extracted flow in test case 1-2 (Maricopa Stanfield).

Figure 13

Demanded extracted flow in pools 1–8: extracted flow by CSE and real extracted flow in test case 1-2 (Maricopa Stanfield).

Numerical example in a laboratory canal

In this example, the CSE algorithm has been tested in a real canal. We want to verify that CSE is able to assess with accuracy a disturbance introduced in a laboratory canal, and in this way, some tests were done to check the good results obtained by CSE in the last examples.

The PAC-UPC is a laboratory canal specially designed to develop basic and applied research in irrigation canal control area and in all subjacent areas like irrigation canal instrumentation, irrigation canal modelling, water measurements, etc. The canal is located in the Laboratory of Physical Models inside the North Campus of the UPC (Barcelona's School of Civil Engineering).

General description

The original idea was to build a canal that could show notorious transport delays in order to behave as closely as possible like a real irrigation canal. The canal was constructed with a zero longitudinal slope and sufficient length. Due to the lack of space inside the laboratory, the canal was designed with a serpentine shape to achieve the maximum canal length in the available space. With this particular design, the result was a 220 m long, 44 cm wide and 1 m deep rectangular canal. A detailed scheme of this laboratory canal is presented in Figure 14 and Table 9.
Table 9

Features of the PAC-UPC canal in the example

Pool Pool length (m) Canal depth (m) Manning's coefficient (n) Width (m) 
1 (from G1 to G242 1.00 0.016 0.44 
2 (from G2 to G345 1.00 0.016 0.44 
3 (from G3 to G445 1.00 0.016 0.44 
4 (from G4 to G545 1.00 0.015 0.44 
5 (from G5 to W443 1.00 0.016 0.44 
Pool Pool length (m) Canal depth (m) Manning's coefficient (n) Width (m) 
1 (from G1 to G242 1.00 0.016 0.44 
2 (from G2 to G345 1.00 0.016 0.44 
3 (from G3 to G445 1.00 0.016 0.44 
4 (from G4 to G545 1.00 0.015 0.44 
5 (from G5 to W443 1.00 0.016 0.44 
Figure 14

Schematic diagram of the top view of the PAC-UPC canal.

Figure 14

Schematic diagram of the top view of the PAC-UPC canal.

The canal is supplied by a small reservoir at the upstream end. The objective of this element is to provide the canal with enough water to feed the canal. Water comes from the reservoir through gate 1 (G1), which is normally in submerged conditions. This gate regulates the inflow by adjusting the gate opening.

The water that is not used is recirculated to the laboratory pumping system. It is possible for the user to arrange the canal with several pool configurations, i.e., a canal with only one very long pool, a canal with one long pool and one short pool, etc.

Test geometry canal

In this example, the geometrical configuration of the canal is as follows: G1 is the only sluice-gate operating because gates 3 and 5 are out of the water. Only the weirs W2 and W4 are operating, and the algorithm only uses the water level measurements at sensors L1, L6, L10 and L11, although the data of sensor L7 were used to check some results.

Gate 1 is made of methacrylate reinforced with a metal skeleton in order to provide enough stiffness and a low weight. The vertical movement of the gate is guided by metal frameworks embedded in the canal and is executed by three-phase servomotors. The gate features in this example are shown in Table 10.

Table 10

Features of gate 1 of the PAC-UPC canal

Gate Gate discharge coefficient Contraction coefficient L1 water level reservoir Gate width (m) Height of the gate opening (m) Steep (m) 
0.68 0.60 1.257 0.434 0.122 0.0 
Gate Gate discharge coefficient Contraction coefficient L1 water level reservoir Gate width (m) Height of the gate opening (m) Steep (m) 
0.68 0.60 1.257 0.434 0.122 0.0 

Rectangular weirs are used to extract water laterally in order to emulate the effect of offtake discharges in real irrigation canals. These weirs have a width of 39 cm and were constructed starting from a 35 cm canal height, except for the end weir (W4) that starts from the canal invert (Sepúlveda 2007). From this minimum height of 35 cm, it is possible to increase the height of the weir by placing PVC pieces in metal rails, one on top of the other. These pieces are 5 cm, 10 cm, 20 cm and 35 cm. With different combinations, it is possible to cover a broad range of weir heights, and thereby to achieve a broad range of output flow through the weir. This structure allows easy measurement of the extracting flow rate measuring the water level upstream of the weir. The features of the sharp-crested weir are shown in Table 11.

Table 11

Features of weir 4 of the PAC-UPC canal

Weir Weir discharge coefficient (CWWeir height (m) Weir width (m) 
0.5776 0.35 0.39 
Weir Weir discharge coefficient (CWWeir height (m) Weir width (m) 
0.5776 0.35 0.39 

Initial conditions for the example

A particular steady state is the initial condition for the canal. The total flow is 110 L/s through G1, and weirs 1, 2 and 3 are not operative, only W4 works. Table 12 shows the water level measurements, which were measured manually, and the flow rate at particular points at initial time. The water level upstream from G1 and the height of the gate opening were shown in Table 10.

Table 12

Initial conditions (water level at particular points)

Checkpoint Initial water level (m) Flow (m3/s) 
1 (L40.758 0.110 
2 (L60.730 0.110 
3 (L80.689 0.110 
4 (L100.644 0.110 
5 (L110.604 0.110 
Checkpoint Initial water level (m) Flow (m3/s) 
1 (L40.758 0.110 
2 (L60.730 0.110 
3 (L80.689 0.110 
4 (L100.644 0.110 
5 (L110.604 0.110 

Scenario

At the beginning of the test, the flow rate in the canal is 110 L/s. At a particular time (250 s after starting the test), some of the pieces that make up the lateral weir (W2) were removed, so the weir height changed to 55 cm. Later, at time 1950 s, the weir was closed again (Table 13).

Table 13

The disturbance characteristics introduced into the system

Time (s) G1 opening (m) W1 height (m) W2 height (m) W3 height (m) W4 height (m) 
0.122 0.90 0.90 0.90 0.35 
250 0.122 0.90 0.55 0.90 0.35 
1,950 0.122 0.90 0.90 0.90 0.35 
3,410 0.122 0.90 0.90 0.90 0.35 
Time (s) G1 opening (m) W1 height (m) W2 height (m) W3 height (m) W4 height (m) 
0.122 0.90 0.90 0.90 0.35 
250 0.122 0.90 0.55 0.90 0.35 
1,950 0.122 0.90 0.90 0.90 0.35 
3,410 0.122 0.90 0.90 0.90 0.35 

There was no flow meter in the canal, and for this reason the exact value of the flow through W2 was not measured directly. However, it was possible to estimate the flow through W2 because we obtained the water level measurements of the L7 sensor (Figure 15), which is the closest sensor to W2, and the discharge coefficient of W2 was calibrated in previous works (see Horváth 2013).
Figure 15

Water level measurements in sensor L7.

Figure 15

Water level measurements in sensor L7.

Table 13 shows the changes in the weir height to introduce disturbance into the canal.

Figure 15 shows the exact time when the disturbance into the canal was introduced, because the water level at sensor L7 decreases quickly by 3 cm. Obviously, the flow extracted through the weir modified the water level surface and the flow along the canal. From sensors L1, L6, L10 and L11 were obtained the water level measurements used by CSE (see Figures 16 and 17).
Figure 16

Water level measurements in sensor L1 (upstream reservoir).

Figure 16

Water level measurements in sensor L1 (upstream reservoir).

Figure 17

Water level measurements in sensors L6, L10 and L11.

Figure 17

Water level measurements in sensors L6, L10 and L11.

The water level measurements at the reservoir collected by sensor L1 were used by CSE for the establishment of the upstream boundary condition. Although the water level at the reservoir should be constant in time, due to limitations of laboratory installations the water level at the reservoir is variable in time as it depends on the water level downstream of the gate.

Results

The disturbances are introduced into the system by modifying the weir height. Sensors L6, L10, L11 get the water level measurements along the canal and these values are introduced into the CSE algorithm, which calculates a discharge through the weir that generates a variation in the water levels at the checkpoints equal to the water level measured at sensors L6, L10 and L11. The extracted hydrograph explains the evolution of the water level measurements at the sensors during the past time horizon.

The hydrograph obtained by CSE was filtered using a ten time steps moving average (see Figure 18). On the other hand, this figure also shows the flow through the weir (W2) calculated using the equation of a sharp-crested weir (eq. weir) from the water level measurements at sensor L7 and the discharge coefficient calibrated by several authors (Galvis et al. 2011; Horváth et al. 2012; Horváth 2013).
Figure 18

Hydrograph obtained by CSE vs. hydrograph obtained using the equation of a sharp-crested weir.

Figure 18

Hydrograph obtained by CSE vs. hydrograph obtained using the equation of a sharp-crested weir.

After analyzing Figure 18, the following can be added:

  • The algorithm obtains an extracted hydrograph similar to the real flow extracted through the weir (W2), especially at the initial moment of introducing the disturbance.

  • Although the real extracted flow was obtained using the weir equation with the water level measurements of sensor L7 and the discharge coefficient calibrated by Hórvarth, the calculated flow extraction was quite accurate due to the noise error of the measurements not being significant.

  • The flow rate difference between the hydrograph obtained by CSE and the real hydrograph was around 2.5 L/s.

  • Probably, the differences between both hydrographs may be the result of deviations in water level measurements at sensors L6, L7, L10 and L11 and/or little deviation between the real Manning or hydraulic loss coefficients and those simulated, although several tests were done to adjust the coefficients as well as possible in the current conditions of the canal.

CONCLUSIONS

The CSE algorithm is very stable numerically because it calculates the extracted hydrograph more similarly to the real hydrograph that has been extracted by the weir. The result is physically feasible, and CSE does not look for incoherent solutions that could also define unreal extraction flows and canal states. The results obtained in the numerical examples verify the good expectations for the algorithm.

The water levels measurements must be as accurate as possible, since the result (extracted hydrograph) obtained by CSE will be as accurate as the input data (water level measurements) are.

The CSE algorithm is sensitive to some physical parameters such as the Manning roughness coefficient and, by extension, other parameters such as the local energy losses in the canal bend. For this reason, the algorithm should not be used in canals where the physical parameters are not well known. Although it is beyond the scope of this paper (it is difficult to always determine the true conditions of any canal, as opposed to simulation models where they are presumed to be known exactly), we did some tests introducing changes in the Manning coefficient, measurement errors in depth gauges or by missing water level measurements of gauges (see Bonet 2015) to check the robustness of CSE.

The CSE is able to calculate the hydrodynamic canal state as shown in the first example. The accuracy of CSE in calculating the extracted flow vector is essential for calculating with accuracy the hydrodynamic canal state with CSE. The CSE calculates the hydrodynamic canal state from the scheduled demands, gate trajectories, initial conditions and the extracted flow vector. The hydrodynamic canal states obtained by a computer model using the real disturbances or from the results obtained by CSE are very similar.

REFERENCES

REFERENCES
Ames
W. F.
1977
Numerical Methods for Partial Differential Equations
.
Academic Press
,
Orlando, FL
,
USA
.
ASCE
1993
Unsteady-flow modeling of irrigation canals
.
Journal of Irrigation and Drainage Engineering
119
(
4
),
615
630
.
Bonet
E.
2015
Experimental Design and Verification of a Centralize Controller for Irrigation Canals
.
PhD Thesis
,
UPC, Polytechnic University of Catalonia
,
Spain
.
Clemmens
A. J.
Kacerek
T. F.
Grawitz
B.
Schuurmans
W.
1998
Test cases for canals control algorithms
.
Journal of Irrigation and Drainage Engineering
124
(
1
),
23
30
.
Crandall
S. H.
1956
Engineering Analisys-A Survey of Numerical Procedures
.
Wiley Interescience
,
New York
,
USA
.
Delgoda
D. K.
Saleem
S. K.
Halgamuge
M. N.
2013
Multiple model predictive flood control in regulated river system with uncertain inflows
.
Water Resources Management
27
(
3
),
765
790
.
Faurès
J. M.
Hoogeveen
J.
Bruinsma
J.
2010
The FAO irrigated area forecasted for 2030
.
(accessed 22 November 2015)
.
Fletcher
R.
1987
Practical Methods of Optimization
,
2nd edn
.
John Wiley & Sons
,
Chichester
,
UK
.
Galvis
R. E.
Horvath
K.
Gómez
V. M.
Rodellar
B. J.
2011
Simplified Modeling of a Laboratory Irrigation Canal for Control Purposes
.
IV Seminar for Advanced Industrial Control Applications
,
Barcelona
,
Spain
, pp.
83
88
.
Gill
P. E.
Murray
W.
Wright
M. H.
1982
Practical Optimization
.
Academic Press
,
London
,
UK
.
Gómez
M
.
1988
The Contribution of Case Study of Transient Flow in Sewer Systems
.
PhD Thesis
,
UPC
,
Barcelona
,
Spain
.
Goussard
J.
1993
Automation of Canal Irrigation Systems. International Commission on Irrigation and Drainage
,
ICID
,
New Delhi
,
India
.
Horváth
K.
2013
Model Predictive Control of Resonance Sensitive Irrigations Canals
.
PhD Thesis
,
Polytechnic University of Catalonia
,
Spain
.
Horváth
K.
Van Overloop
P. J.
Galvis
E.
Gómez
M.
Rodellar
J.
2012
Multivariable Model Predictive Control of Water Levels on a Laboratory Canal
.
SimHydro 2012: New frontiers of simulation
,
Sophia Antipolis, Nice
,
France
.
Institut Mediterranéen de l'Eau
2007
3 Éme atelier règional. In: La demande en eau en Méditerranée, Zaragoza, Spain, 19–21 March, p. 18
.
Liu
F.
Feyen
J.
Malaterre
P. O.
Baume
J. P.
Kosuth
P.
1998
Development and evaluation of canal automation algorithm CLIS
.
Journal of Irrigation and Drainage Engineering
124
(
1
),
40
46
.
Malaterre
P. O.
1994
Modelisation, Analysis and LQR Optimal Control of an Irrigation Canal
.
PhD Thesis
,
LAAS-CNRS-ENGREF-Cemagref, Etude EEE n14
,
France
.
Malaterre
P. O.
Rogers
D. C.
2008
Teaching canal hydraulics and control using a computer game or a scale model canal
. In:
13th IWRA World Water Congress
,
Montpellier
,
France
, p.
11
.
Mareels
I.
Weyer
E.
Ooi
S. K.
Cantoni
M.
Li
Y.
Nair
G.
2005
Systems engineering for irrigation systems: successes and challenges
.
Annual Reviews in Control
29
(
2
),
191
204
.
Plusquellec
H.
1988
Improving the Operation of Canal Irrigation Systems. A Seven-Part Slide Series
.
The Economic Development Institute of the World Bank
,
Washington, DC
,
USA
.
Press
H. W.
Teukolsky
S. A.
Vetterling
W. T.
Flannery
B. P.
1992
Numerical Recipes in FORTRAN 77
.
Cambridge University Press
,
Cambridge
,
UK
.
Sepúlveda
C.
2007
Instrumentation, Model Identification and Control of an Experimental Irrigation Canal
.
PhD Thesis
,
Polytechnic University of Catalonia
,
Spain
.
Soler
J.
2003
Study About Controllers in Irrigation Canals Using Nonlinear Numerical Methods
.
PhD Thesis
,
Polytechnic University of Catalonia
,
Spain
.
Soler
J.
Gómez
M.
Rodellar
J.
2013
Goroso: Feedforward control algorithm for irrigation canals based on sequential quadratic programming
.
Journal of Irrigation and Drainage Engineering
139
(
1
),
41
54
.
Strelkoff
T.
1970
Numerical solution of Saint-Venant equations
.
Journal of the Hydraulics Division, ASCE
96
,
223
252
.
Van Overloop
P. J.
Weijs
S.
Dijkstra
S.
2008
Multiple model predictive control on a drainage canal system
.
Control Engineering Practice
16
(
5
),
531
540
.
Walker
W. R.
Skogerboe
G. V.
1987
Surface Irrigation. Theory and Practice
.
Prentice-Hall
,
Englewood Cliffs, NJ
,
USA
.
Wylie
E. B.
1969
Control of transient free-surface flow
.
Journal of the Hydraulics Division ASCE
95
,
347
362
.