Modal decomposition of linearized open channel flow

Open channel flow is traditionally modeled as an hyperbolic system of conservation laws, 
which is an infinite dimensional system with complex dynamics. We consider in this paper 
an open channel represented by the Saint-Venant equations linearized around a non uniform 
steady flow regime. We use a frequency domain approach to fully characterize the open 
channel flow dynamics. The use of the Laplace transform enables us to derive the 
distributed transfer matrix, linking the boundary inputs to the state of the system. The 
poles of the system are then computed analytically, and each transfer function is 
decomposed in a series of eigenfunctions, where the influence of space and time variables 
can be decoupled. As a result, we can express the time-domain response of the whole canal 
pool to boundary inputs in terms of discharges. This study is first done in the uniform 
case, and finally extended to the non uniform case. The solution is studied and 
illustrated on two different canal pools.

Abstract. Open channel flow is traditionally modeled as an hyperbolic system of conservation laws, which is an infinite dimensional system with complex dynamics. We consider in this paper an open channel represented by the Saint-Venant equations linearized around a non uniform steady flow regime. We use a frequency domain approach to fully characterize the open channel flow dynamics. The use of the Laplace transform enables us to derive the distributed transfer matrix, linking the boundary inputs to the state of the system. The poles of the system are then computed analytically, and each transfer function is decomposed in a series of eigenfunctions, where the influence of space and time variables can be decoupled. As a result, we can express the time-domain response of the whole canal pool to boundary inputs in terms of discharges. This study is first done in the uniform case, and finally extended to the non uniform case. The solution is studied and illustrated on two different canal pools.
approach with time-domain simulations. The paper [17] and related studies provided a dimensionless expression of Saint-Venant equations, and time domain simulations were performed in order to evaluate the effect of downstream boundary condition on the flow routing. However, these studies did not end with general results for the characterization of open channel flow.
We propose to tackle these questions by using a frequency domain approach, which gives a better insight into dynamical systems' properties than the time domain approach. We follow in this way the approach of automatic control engineers, but we will also illustrate the analysis in the time domain, in order to show the link between the two methods. Both methods are strongly related, but this is unfortunately too rarely pointed out (see e.g. [7] for frequency domain analysis of Saint-Venant equations). Interesting results were obtained by Dooge et al. [6] concerning the dynamic analysis of linearized Saint-Venant equations for a semiinfinite channel around uniform flow, also using a frequency domain approach. Tsai [20] provided a detailed analysis of Saint-Venant equations and simplified models for flow routing accounting for backwater effects. However, the developed criteria uses a local characterization of the flow, and does not provide clues for developing an integrated criteria based on these results. As we will show in the sequel, the dynamics of a reach are strongly linked to the poles of the system's transfer matrix, which are themselves linked to the reach characteristics (length, geometry, etc.). Therefore, our results show that an integrated criteria is necessary to characterize the reach dynamics.
Moussa and Bocquillon [13] and [14] have used a frequency analysis of the linearized Saint-Venant equations to provide a criteria for the choice of a simplified model. They considered the relative weight of different terms in the solutions of the equations in order to justify the choice of simplified models. This approach does not characterize the flow dynamics per se, but is based on a chosen threshold to separate the different models.
Ridolfi et al. [15] computed the Green's function of the linearized Saint-Venant equations, which gives interesting results in terms of characterization of the flow. However, this does not provide a characterization of a given pool.
Thirriot and Benayada [19] used the Laplace transform to derive an approximate solution in the uniform case, but they could not compute the poles in an explicit way because of the boundary condition they used. With a downstream boundary condition in terms of discharge, we have shown that the poles can be computed explicitly (see [10]).
The approach developed in this paper is based on the poles of the system, which enable to fully characterize its dynamic behavior. The study is focused on the linearized equations, which capture the dynamics around stationary regimes. We provide a complete characterization in the uniform flow case, where closed-form solutions are available, and extend the results to the non uniform case using a numerical approach [10]. We show that the non uniform flow case is mathematically speaking very close to the uniform flow case, even if closed-form solutions are no longer available. This opens the way to a proper classification of open channel flow characteristics, which could be based on the modal decomposition proposed in this paper.
The outline of the paper is as follows: after this introduction, Section 2 presents the Saint-Venant equations of open channel flow and the linearized equations around non uniform flow. Section 3 introduces the frequency domain approach for the linearized model around a uniform steady flow, with a complete characterization of the transfer matrix in terms of poles, delays and modal factors. These results are then extended to the linearized equations around a non uniform steady flow in Section 4, using a numerical approach.

Modeling of Open Channel Flow.
2.1. Saint-Venant Equations. The Saint-Venant equations are two coupled partial derivative equations involving the discharge Q(x, t) (m 3 /s) and the water depth Y (x, t) (m). The first one is the mass conservation equation: and the second one is the momentum conservation equation: is the wetted area, and g (m/s 2 ) is the gravitational acceleration. The friction slope S f (x, t) (m/m) is modeled with the classical Manning formula [2]: with n the Manning coefficient (sm −1/3 ) and R(x, t) the hydraulic radius (m), defined by R(x, t) = A(x, t)/P (x, t), where P (x, t) is the wetted perimeter (m) (see Fig. 1). To complete the equations, we need to introduce initial and boundary conditions. The initial condition is given in terms of (Q(x, 0), Y (x, 0)), for all x ∈ [0, L], with L the length of the channel. (1) and (2). Then, denoting the variables corresponding to the equilibrium regime with a subscript zero (Q 0 (x), Y 0 (x), etc.), the Saint-Venant equations become: A0(x) the flow velocity (m/s). These two equations define an equilibrium regime given by Q 0 (x) = Q 0 and Y 0 (x) solution of the ordinary differential equation (4b), for a boundary condition in terms of downstream elevation.
When the right-hand side of equation (4b) is equal to zero, the water depth is constant along the channel. In this case, given Q 0 (x) = Q 0 , the equilibrium solution Y 0 (x) = Y n (also called normal depth) can be deduced by solving the following algebraic equation: This specific solution corresponds to the uniform flow regime. Equation (5) is usually solved numerically with a fixed-point or Newton-Raphson method.
In some cases, the uniform depth can be computed analytically. For large rectangular channels, the hydraulic radius R can be approximated by the water depth Y , and the Manning equation (3) reduces to: Combining equations (5) and (6) give the uniform depth Y n corresponding to a discharge Q 0 in the large rectangular case: The uniform flow regime will be specifically developed as an example, since this regime leads to closed-form analytical solutions for the transfer matrix. However, we will also show that this specific flow regime is not qualitatively different from realistic non uniform flow regimes.
The article will be illustrated on two trapezoidal prismatic channels, with different characteristics (see Fig. 2).

Figure 2. Section of a trapezoidal canal
The two example canals are presented in the following. Example 1. The channel characteristics are given in table 1, where L is the channel length (m), m the bank slope, B the bed width (m), S b the bed slope (m/m), n the Manning coefficient (m −1/3 s), Y n the normal depth (m) corresponding to the discharge Q 0 (m 3 s −1 ). Canal 1 is a short oscillating canal, and canal 2 is a long delayed canal.
Keeping a constant discharge, we compare the backwater curves obtained for different downstream boundary conditions Y 0 (L) = Y n × [0.8, 1, 1.2] (see figure 3). We have two types of backwater curves: if Y 0 (L) > Y n , then the flow is decelerating along x, this is a so-called 'M1' curve (see e.g. [18]), and if Y 0 (L) < Y n , then the   (1) and (2) and expanding in series.
After collecting terms in q, a, ∂q ∂x and ∂a ∂x , the linearized equations are obtained as follows: The solution in terms of variations of water level y(x, t) can readily be deduced by dividing by T 0 (x) the solution obtained in terms of variations of wetted area a(x, t) = T 0 (x)y(x, t).
To facilitate the mathematical analysis, we rewrite the linearized Saint-Venant equations as follows: The initial condition and the boundary conditions are given by: We also assume that the water level deviations are measured at each boundary: y(0, t) and y(L, t).

Uniform Flow Case.
In the uniform flow case, the parameters α, β, γ and δ are constant and given by: 3.1. Derivation of the Transfer Matrix. We will now derive the transfer matrix for the linearized Saint-Venant equations, taking q(0, t) and q(L, t) as boundary inputs, and the water depth deviations y(x, t) and discharge deviations q(x, t) as outputs.
We apply Laplace transform to the linear partial differential equations (9), using the classical relationd f dt = sf(s) − f (0) with s the Laplace variable, which yields: with A(s) = −A −1 (sI + B) and B = A −1 , i.e.: Let us diagonalize matrix A(s): and where λ 1 (s) and λ 2 (s) are the eigenvalues of A(s), solutions of the equation: This equation is a second order polynomial equation in λ: which has in general two solutions: 3.1.1. State-Transition Matrix. Once matrix A(s) is diagonalized, the differential equation (11) can be solved analytically: The elements of matrix Φ(x, s) are given by: Φ(x, s) is the state-transition matrix for the differential equation (11). For simplicity, we assume zero initial conditions for the following developments, thereforē ξ 0 (x, s) = 0.
3.1.2. Boundary Conditions. Specifying the boundary conditions leads to express the state as a function of the boundary discharges. To this end, we use Eq. (16) for x = 0, and get the equality: Once again, this is valid provided φ 21 (L, s) is not equal to zero. The values of s such that φ 21 (L, s) = 0 correspond to the poles of the Saint-Venant transfer matrix.
If s 0 is such that λ 1 (s 0 ) = λ 2 (s 0 ), we have: When λ 1 (s) = λ 2 (s), equation (20) is equivalent to There is a pole in zero p 0 = 0, which means that the canal pool acts as an integrator and the other poles verify the equation: with k ∈ N * . The case k = 0 leads to λ 1 (s) = λ 2 (s), therefore s = s 0 , which is not a pole. Finally, the poles are solutions of the second order equation: The poles (p ±k ) k∈N * are then given by: Then the poles obtained for 0 < k ≤ k m are negative real, and those obtained for k > k m are complex conjugate, with a constant real part equal to − (α−β)γ+2αβδ , which means that they are located on a vertical line in the left half plane. Canal pools with a dominant oscillating behavior correspond to Δ(1) < 0.
Example 2 (Poles of the example canals). The poles of the canals 1 and 2 are depicted in Fig. 4. Canal 1 has an oscillating behavior, all its poles being complex conjugate. Canal 2 has two negative real poles for k m = 1, and the other ones are complex conjugate. The complex poles of canal 2 have a larger negative real part than the ones of canal 1. It should therefore exhibit a damped behavior compared to the canal 1. In the Laplace domain, P (s) = (p ij (s)) denotes the 2 × 2 transfer matrix relating the inputs to the outputs. Using the algebraic relation (19) and the state-transition matrix (18), we obtain the input-output transfer matrix: This result is similar to the expressions obtained by several authors [3,4,16,1]. We depict in the following example the transfer functions p 21 (s) and p 22 (s) which are useful in a context of automatic control of open channels.
Example 3 (Bode plots of the input output transfer matrix). Bode plots of transfer functions p 21 (s) and p 22 (s) for example canals 1 and 2 are depicted in Figs. 5 and 6. We observe that the low frequency behavior of the system is dominated by the integrator, that there is a delay in transfer function p 21 (s). In high frequencies, the oscillating modes of canal 1 are clearly visible, while canal 2 is damped and tends towards a constant gain.  with The input-output transfer matrix P (s) can be obtained from the distributed transfer matrix G(x, s) as follows: p 11 (s) = g 11 (0, s), p 12 (s) = g 11 (L, s), p 21 (s) = g 12 (0, s) and p 22 (s) = g 12 (L, s). In the spatial Bode plot of canal 1, the integrator is clearly visible in low frequencies. For higher frequencies, the oscillating modes can be seen with the corresponding nodes and anti-nodes of oscillations along the channel. The canal 2 is much more damped than canal 1, but there are oscillating modes at the downstream end of the canal. This means that the water level response to the upstream discharge may exhibit local oscillations close to the downstream end of the canal. The surprising thing here is that these oscillations do not appear in the input-output transfer functions (see Fig. 6).

Properties of the Transfer Matrix.
3.2.1. Delays.  We finally get the following factorization of all elements of the transfer matrix: Let us note that when |s| → ∞, we have This explains why the terms e τ1s+λ1(s)L and e τ2s−λ2(s)L are delay free, since a delay can be pointed out in high frequencies by its effect on the phase of the transfer function. This result shows that the Saint-Venant transfer matrix includes delays only linked to the waves propagations, which correspond to the delays obtained by the characteristics.
Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi: 10.3934/nhm.2009.4.325 Distributed Transfer Matrix. The delays are obtained directly as e − x α s for the distributed transfer functions g 11 (x, s) and g 21 (x, s) and e − L−x β s for the distributed transfer functions g 12 (x, s) and g 22 (x, s). The transfer functions can therefore be factorized as where the delay-free partsg ij (x, s) are given by:

3.2.2.
Asymptotic Estimate of the Poles. The following proposition proved in [11] provides an asymptotic estimate of the poles for high frequencies.
This proposition shows that for high frequencies the poles of Saint-Venant transfer matrix are close to the ones of the following damped wave equation: with boundary conditions q(0, t) and q(L, t). Using Laplace transform, this equation reduces to an ODE in x, with eigenvalues equal to −r 1 − s α and r 2 + s β . The obtained transfer function has for denominator 1 − e −(r1+r2)L−(τ1+τ2)s , whose roots coincide with the poles approximation (26). This shows that the oscillating modes correspond to the interaction of two gravity waves, one traveling downstream at speed α = V 0 + C 0 with attenuation factor r 1 , and one traveling upstream at speed β = C 0 − V 0 with attenuation factor r 2 .

Modal Decompositions.
Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi: 10.3934/nhm.2009.4.325 Rational Series Expansion. For simplicity, we assume in the following that the poles have single multiplicity, i.e. that Δ(k) = 0, but the solution can easily be extended to the case where Δ(k) = 0. Then, the Cauchy residues theorem implies that each transfer function g ij (x, s) can be decomposed as an infinite sum: and The coefficient a ij (x) is the residue of transfer function g ij (x, s) at the pole p k .
Proof. The proof can be easily adapted from the one in Appendix A.
The residues can be computed analytically with Eq. (28), leading to, for k = 0: with ψ = γL αβ . The modal factors a (k) ij (x) are combination of sine and cosine functions of kπx L , modified by a multiplicative exponential term which depends on the difference α−β and the coefficient γ. Coefficients γ and δ are directly linked to the slope of the open channel. The larger the slope, the larger their influence on the modal factors. Rational Decomposition of the Delay-free Part. We have shown that the transfer functions g ij (x, s) could be factorized as a delay times a delay-free part. Therefore, one can find an approximation of g ij (x, s) by applying the Cauchy residues theorem on transfer functionsg ij (x, s), and then adding the delay. Finally, each transfer functiong ij (x, s) can be decomposed as: The coefficientã ij (x) is the residue of transfer functiong ij (x, s) at the pole p k .
Proof. The proof can be easily adapted from the one in Appendix A.
The residues can be computed analytically with Eq. (34), leading to, for k = 0 and, for k = 0:ã and theb ij (x) are given by: Example 5 (Modal factors of the example canals). The functions a (k) 11 (x) of the distributed modal decomposition of g 11 (x, s) of the example canals 1 and 2 are depicted in Fig. 9 and 10.
For canal 1, which is flat and has oscillating modes, the modes of oscillation are clearly visible. The integrator in this case is slightly modified by the slope and the friction. We clearly see the oscillating modes and the associated nodes and anti-nodes of oscillation. A node corresponds to a point where the traveling waves interfere negatively with each other. An anti-node corresponds to a positive interference: the perturbations are in phase and act additively.  For canal 2, which is steep, we see that the dominant term is the integrator, which is greatly modified by the slope. The poles p ±1 are negative real, therefore not oscillating, and the ones obtained for |k| > 1 are oscillating. The corresponding modal factors are strongly modified by the slope, compared to the ones of the example canal 1. Figure 11 depicts the value of coefficients b 11 (x) and b 12 (x) of the distributed modal decomposition of the transfer functions g 11 (x, s) and g 12 (x, s) for canals 1 and 2.
The term b ij (x) corresponds to a direct term in transfer function g ij (x, s). We see that this term varies a lot with x for both canals. It even changes sign, which is a surprising result: this means that the water level change resulting from a step Due to structural properties of the considered differential equation, the classical compromise between precision of the approximate solution and integration step size necessarily leads to a large computational time for each value of s. The frequency domain of interest goes from the low frequency behavior of the canal (typically ω r /100 with ω r the resonant frequency ω r = 2π/(τ 1 + τ 2 ) with τ 1 = L/(V 0 + C 0 ) and τ 2 = L/(C 0 − V 0 ) respectively the downstream and upstream propagation timedelays) to the high frequency behavior (typically 10 times the resonant frequency). Since it is important to have a good representation in between these extreme frequencies (and especially to correctly reproduce the resonant modes), it is necessary to compute a large number of frequency points (typically 500 points) [10]. To this end, it is therefore essential to have an efficient numerical method to compute the frequency response.

Computation of the Transfer Matrix.
In the non uniform flow case, it is no longer possible to diagonalize the system. Indeed, applying Laplace transform and after elementary manipulations, we obtain the following differential equation: with ξ(x, t) = (T 0 (x)y(x, t), q(x, t)) T and where A(x, s) and B(x) depend on x: To simplify the exposition, we assume in the following that ξ(x, 0) = 0.
One may show that the change of variable that diagonalizes A(x, s) introduces a supplementary term due to the derivative of matrix X (x, s) with respect to x. Indeed, let us introduceζ(x, s) = X (x, s)ξ(x, s), where X (x, s) is such that A(x, s) = X (x, s) −1 D(x, s)X (x, s). Plugging this variable into Eq. (39) leads to: The matrix D(x, s) X (x, s) −1 . Therefore we can no longer provide an analytical solution to this equation.
However, since the differential equation (39) is linear, we know that its general solution always exists, is unique and is given by [9, p. 598]: where Φ(x, s) is the state-transition matrix associated to the differential equation (39). As a matter of fact, one may think that, as in the scalar case, the state-transition matrix could be written as Φ(x, s) = e and therefore: The only case where this holds is when A(x, s) and x 0 A(v, s)dv commute. This is true e.g. when A(x, s) is constant with respect to x, i.e. in the uniform flow case. In the non uniform flow case, this is generally not true, this is why we need to find a numerical way to compute the state-transition matrix.
Once the state-transition matrix is obtained, the transfer matrix corresponding to the original differential equation (39) is then given by: The problem is therefore to solve the differential equation (39) and compute the state-transition matrix Φ(x, s). This problem is not easy to solve, due to its oscillatory and unstable nature.
We proposed an efficient numerical solution to this problem in [10], leading to a general method to compute the Saint-Venant transfer matrix for any kind of non uniform flow. 4.1.1. Poles Computation. As we have seen, the poles correspond to the values of s such that φ 21 (L, s) = 0. In the non uniform case, the poles are obtained numerically, by finding the zeros of equation (20). Since we can compute numerically the transition matrix, we know how to compute φ 21 (L, s) for any s ∈ C, and we can look for values of s such that φ 21 (L, s) = 0.
Practically, we seek the value of s that minimizes the modulus of φ 21 (L, s). This search is done starting from the values obtained analytically for the uniform regime.
We used this numerical method to find the values of s that minimize the modulus of φ 21 (L, s). This has been done for the first 7 pairs of poles for the canals in example 1.
The results are depicted in Figs. 12 and 13 for canals 1 and 2, respectively. In both figures, we see that an accelerating flow tends to dampen the oscillating poles, while a decelerating flow tends to reduce their damping. However, the first negative real pole of canal 2 reacts differently: its damping increases for decelerating flow.
This confirms our results from the asymptotic study in high frequencies, but this does not necessarily apply for low frequencies.  with Example 6 (Bode plots for non uniform flow). We compute the Bode plots for the three situations: uniform flow, decelerating flow (M1 backwater curve) and accelerating flow (M2 backwater curve). The results are depicted in Figs. 14 and 15. The integrator gain increases for accelerating flow, while it decreases for decelerating flow. This also is consistent with the physical intuition: in a first approximation, the integrator gain is inversely proportional to the area of the pool. In the case of a trapezoidal geometry, the area increases when the water level increases. For a constant discharge, this corresponds to a decelerating flow. This explains why the integrator gain tends to decrease for a decelerating flow. The reverse occurs for accelerating flow.
For the oscillating modes, we see that their damping decreases for decelerating flow, while it increases for accelerating flow. This will be confirmed by the direct study of the poles. x ∈ [0, L]. Then, the distributed transfer matrix G(x, s) is directly obtained as follows: ŷ(x, s) q(x, s) = G(x, s) q(0, s) q(L, s) with These expressions generalize the results already obtained in Section 3.1.5 to the non uniform flow case.

4.2.1.
Delays. In the non uniform flow case, the delays of the Saint-Venant transfer functions can be computed by integrating the characteristics lines along the channel. Indeed, we can show that the delays for the transfer functions p 21 (s) and p 12 (s) in the non uniform case are given by: for p 21 (s) and This result is consistent with our physical intuition, and confirms that an open channel in non uniform flow behaves very similarly to an open channel in uniform flow: the main difference lies in the fact that the results can no longer be expressed in a closed form, but the proposed numerical scheme enables to compute the transfer function for any realistic flow configuration.
Example 7 (Delays for non uniform flow). Let us define the dimensionless time t * = αnt L , where α n stands for the uniform value of α = V 0 + C 0 . This enables us to compare the delays τ 1 and τ 2 for various non uniform flow conditions: for the uniform flow, we have τ * 1 = 1, and τ * 2 = αn βn , where the subscript n denotes the uniform flow values. Fig. 16 depicts the way the dimensionless delays vary with the flow in both canals. The two canals have a very different behavior: the delay τ 1 varies much more with the downstream water level in the case of canal 1 than in the case of canal 2. Indeed, the backwater curve affects only a small portion of the canal 2, while it affects the whole canal 1. The delay τ 2 is much larger than τ 1 in the case of canal 2, because β is relatively much lower than α for canal 2 than for canal 1.  Figure 16. Variation of the dimensionless delays τ * 1 and τ * 2 as function of Y 0 (L)/Y n 4.2.2. Asymptotic Estimates of the Poles. It is also possible to derive a closed-form expression for the high frequency behavior of the system.
We start from the characteristics form: where χ(x, t) = X(x)ξ(x, t) with ξ(x, t) = (T 0 (x)y(x, t), q(x, t)) T , and . The elements of matrix E(x) are given by: where we have dropped the argument x for readability, and α and β denote the derivatives of α and β with respect to x.
For high frequencies, the diagonal terms in C(x, s) dominate the anti-diagonal terms. Then, the high frequency approximate solution of Eq. (50) is given by: Therefore, using the previous developments, a high frequency approximation of the poles is obtained by solving: which leads to: . We recover an expression similar to the one already obtained in the uniform case. The high frequency poles have an imaginary part which is linked to the forward and backward delays, and a real part which is linked to the forward and backward damping.
Let us note that: where the argument x has been dropped for readability.
Therefore, the high frequency poles damping is the sum of a positive term plus a term whose sign depends on the derivative of α(x) − β(x) = 2V 0 (x) with respect to x. This means that the system will tend to be more damped if the velocity tends to increase along the canal. On the contrary, it shows that the system will tend to be more oscillating if the velocity tends to decrease along the canal. This is consistent with our intuition: we expect more damping in an accelerating flow, while a decelerating flow will be more prone to possible oscillations.
As in the uniform case, we may show that for high frequencies, the Saint-Venant equations are close to the following damped wave equation: with boundary conditions q(0, t) and q(L, t).
The high frequency modes are therefore the result of the interaction of two gravity waves, one traveling downstream at speed α(x) = V 0 (x)+C 0 (x) with attenuation factor r 1 (x), and one traveling upstream at speed β(x) = C 0 (x) − V 0 (x) with attenuation factor r 2 (x). The high frequency estimate of the poles is the direct extension of the result obtained in the uniform case. We know using Eq. (51) that the high frequency poles are complex conjugate, and close to a vertical line in the left half plane. Using similar arguments as in the uniform case, the Cauchy residues theorem implies that each transfer function g ij (x, s) can be decomposed as an infinite sum: where the coefficient a ij (x) = lim s→p k (s − p k )g ij (x, s) is the residue of transfer function g ij (x, s) at the pole p k , and b ij (x) = ∂ ∂s [sg ij (x, s)]| s=0 .
Contrarily to the uniform case, we have no analytical expression for the residues a (k) ij (x). However, the proposed numerical method also enables us to compute the terms numerically.

Proposition 2.
In the non uniform case, the coefficient of the distributed modal decomposition can be computed as follows: • The coefficients a ij (x) are given by: whereφ 21 (L, 0) is defined by: • The modal coefficients a (k) ij (x) are given by: where the N ij (x, p k ) are given by: and φ 21 (L, p k ) is the term of the matrix Φ (x, s) given by: • The coefficients b ij (x) are given by: where the N ij (x, 0) are given by: φ 21 (L, 0) is given by Eq. (54), andφ 21 (L, 0) is given by: Proof. See Appendix B.  11 (x) of the distributed modal decomposition of g 11 (x, s) of the example canals 1 and 2 are depicted in Fig. 17 and 18. We see that in both canals the modes in non uniform flow are similar to the ones in uniform flow. For canal 1, the spatial functions a  For canal 2, the non uniformity has an important effect on the modal decomposition. The nodes are modified following the same pattern observed for canal 1, but the change in amplitude is much more important, especially in the case of accelerating flow.  Figure 19 depicts the value of coefficients b 11 (x) and b 12 (x) of the distributed modal decomposition of the transfer function g 11 (x, s) and g 12 (x, s) for canals 1 and 2. We observe the same trends already discussed for coefficients a

Example 8 (Coefficients a
where the coefficientã ij (x) = lim s→p k (s − p k )g ij (x, s) is the residue of transfer functiong ij (x, s) at the pole p k .
For the delay-free decomposition, we have:ã

Conclusion.
We have derived the Saint-Venant transfer matrix, first for the uniform flow case, and then for the general non uniform flow case. We have analyzed the transfer matrix for boundary conditions in terms of discharges, and have computed its poles, first analytically for the uniform case, then numerically for the backwater case. We also characterized the model in terms of delays, which gives an a priori evaluation of the available bandwidth for the controlled system. These results are important for the control of open channel flow. As shown in [12], the modal decomposition can be used to show that the transfer matrix belongs to the Callier-Desoer class of transfer matrices. This fact allows us to use classical frequency domain tools such as the Nyquist criteria for closed-loop stability.
Finally, we have shown that the dynamic behavior of an open channel in non uniform flow is qualitatively very similar to its behavior in uniform flow. The only difference is that in non uniform flow, the nice analytical expressions can no longer be used, but efficient numerical methods are available to compute poles, delays and modal factors. This transfer function is well defined in s = 0 since 0 is a root of φ 21 (x, s). Using Eq.(66c),φ 21 (x, s) verifies the following differential equation: with an initial condition given byφ 21 (0, s) = 0. Therefore the value of Φ (x, s) at x = 0 is 0 and we finally get Eq. (57). This equation provides a way to compute Φ (x, s) for any s ∈ C.
Finally, in order to compute the coefficients a (k) ij (x), we need to compute the N ij (x, p k ). For any s ∈ C, the N ij (x, s) are given by: The N ij (x, 0) can be computed using Eqs. (70) and the solution Φ (x, 0) given by Eq. (57) for s = 0. In that case, the expressions simplify, and lead to Eqs. (59).
To computeφ 21 (L, 0), we differentiate Eq. (68) with respect to s, which gives: with an initial condition given byφ 21 (0, s) = 0. Therefore, once φ 11 (x, 0) is computed,φ 21 (L, 0) is directly obtained by integration, as in Eq. (60). Therefore, we have derived a numerical means to compute all the coefficients of the modal decomposition in the non uniform case, which ends the proof.