Wael Suleiman
email: suleiman@laas.fr

André Monin
email: monin@laas.fr

M Wael

Suleiman

New Method for Identifying Finite Degree Volterra Series ⋆

Keywords: Identification, nonlinear systems, Volterra series, optimization

In this paper, the identification of a class of nonlinear systems which admits input-output maps described by a finite degree Volterra series is considered. In actual fact, it appears that this class can model many important nonlinear multivariable processes not only in engineering, but also in biology, socio-economics, and ecology.

To solve this identification problem, we propose a method based on a local gradient search in a local parameterization of the state space realization of finite degree Volterra series with infinite horizon. Using the local parameterization not only reduces the amount of the gradient calculations to the minimal value, but also overcomes the nonuniqueness problem of the optimal solution. Moreover, we propose a sequential projection method to provide an initial estimation of the parameters of finite degree Volterra series realization. This estimation is used to initialize the gradient search method.

Introduction

During the last thirty years, the development of identification theory of dynamic systems has been a subject of active research, such a development is the result of the necessity of reliable models in process control, aerospace applications, biomedical systems, ecology, physiology, biology, and socio-economics [START_REF] Lakshmikantham | Stability Analysis of Nonlinear Systems[END_REF][START_REF] Nijmeijer | Nonlinear Dynamical Control Systems[END_REF][START_REF] Sastry | Nonlinear Systems: Analysis, Stability and Control[END_REF][START_REF] Schetzen | Nonlinear system modeling based on the wiener theory[END_REF][START_REF] Korenberg | The identification of nonlinear biological systems: Wiener kernel approaches[END_REF][START_REF] Schuppen | System theory for system identification[END_REF].

The developments of linear subspace identification methods have recently offered a significantly practical tool to deal with the identification of multivariable linear or pseudo linear system [START_REF] Söderström | System identification[END_REF][START_REF] Moonen | On-and off-line identification of linear statespace models[END_REF][START_REF] Overschee | N4SID: subspace algorithms for the identification of combined deterministicstochastic systems[END_REF][START_REF] Viberg | Subspace-based methods for the identification of linear time-invariant systems[END_REF][START_REF] Ljung | System Identification: Theory for the User[END_REF]. However, in practice, the identification of nonlinear multivariable systems is becoming crucial since many systems in nature are nonlinear. In the theory of nonlinear systems, the term "nonlinear" defines a class of systems for which the linear approximation fails to be an efficient model able to capture the dynamic of the system. The class of nonlinear systems is a complex class, and defining an universal represented model of this class is a very complicated task. Therefore, anyone interested in nonlinear modeling is compelled to focus on specific model classes [START_REF] Billings | Identification of nonlinear systems: a suvey[END_REF][START_REF] Haber | Structure identification of nonlinear dynamic systems -a survey on input/output approaches[END_REF][START_REF] Lee | Modeling and identification for nonlinear predictive control: Requirements, current status and future research needs[END_REF][START_REF] Mathews | Polynomial Signal Processing[END_REF][START_REF] Pearson | Nonlinear input/output modeling[END_REF][START_REF] Pearson | Discrete-Time Dynamic Models[END_REF][START_REF] Doyle | Identification and Control Using Volterra Models[END_REF].

In this paper, we focus on the representation of nonlinear system by Volterra series with infinite horizon [START_REF] Volterra | Theory of Functionals and of Integral and Integro-Differential Equations[END_REF][START_REF] Brockett | System theory on group manifolds and coset spaces[END_REF][START_REF] Brockett | Volterra series and geometric control theory[END_REF]. The expansion of a finite degree Volterra series has the following form

y t = l ∑ n=1 ∞ ∑ τ 1 =0 τ 1 ∑ τ 2 =0 . . . τ n-1 ∑ τ n =0 w n (τ 1 , τ 2 , . . . , τ n ) ×u t-τ 1 ⊗ u t-τ 2 ⊗ . . . ⊗ u t-τ n + v t (1)
where u t ∈ R m is the input signals, y t ∈ R p is the output signals. The measurement noise v t is assumed to be a whitenoise that is independent of the input signal and with zero mean, and ⊗ denotes the Kronecker product. The functions w n (τ 1 , τ 2 , . . . , τ n ) ∈ R p×m n denote the Volterra kernels. As the last term in the series involves the l th kernel, they will be called Volterra series of degree l with infinite horizon. This class of Volterra series is dense, in the L 2 sense, in nonlinear analytic input-output systems [START_REF] Brockett | Volterra series and geometric control theory[END_REF]. Approximating nonlinear systems by Volterra series has been used in many research areas, such as filtering [START_REF] Monin | Volterra filtering with application to bilinear systems[END_REF], predictive control [START_REF] Allgöwer | Nonlinear Model Predictive Control[END_REF], ecology [START_REF] Takeuchi | Global Dynamical Properties of Lotka-Volterra Systems[END_REF], biology [START_REF] Thieme | Mathematics in Population Biology[END_REF], and system control [START_REF] Nijmeijer | Nonlinear Dynamical Control Systems[END_REF][START_REF] Doyle | Identification and Control Using Volterra Models[END_REF]. In reality, the direct calculation of the Volterra kernels from Equation ( 1) is not feasible in practice. This because it involves the integral past of input signal and the size of the matrices which should be manipulated increases exponentially. For that, the horizon of finite degree Volterra series is fixed [START_REF] Rugh | Nonlinear System Theory: The Volterra/Wiener Approach[END_REF][START_REF] Schetzen | The Volterra and Wiener theories of nonlinear systems[END_REF]. The obtained class is called truncated finite degree Volterra series of degree l. It has the following form

y t = l ∑ n=1 T ∑ τ 1 =0 τ 1 ∑ τ 2 =0 . . . τ n-1 ∑ τ n =0 w n (τ 1 , τ 2 , . . . , τ n ) ×u t-τ 1 ⊗ u t-τ 2 ⊗ . . . ⊗ u t-τ n + v t (2)
where T is the length of the considered horizon. In this case, the estimation of Volterra kernels becomes trivial and can be obtained by least-square methods. However, if we consider that the Volterra kernels are fully parameterized, then the number of parameters which should be estimated is the following 1

C (l, T ) = p l ∑ n=1 (T + n)! n! T ! m n (3)
recall that p and m are the dimensions of the output and input signals respectively. It is clear that C (l, T ) increases exponentially with respect to l and T . For that, in the case of multivariable systems, we are obligated to operate high dimensional matrices, which are often ill-conditioned [START_REF] Nowak | Efficient methods for identification of Volterra filter models[END_REF].

The main contribution of this paper is defining a new method to identify a state-space realization of finite degree Volterra series with infinite horizon (1). The initial estimation of the realization's parameters is obtained by a sequential projection method that we have developed thanks to the recursive property of the realization structure, then the realization parameters are optimized using a local gradient search method.

The paper is organized as follows. In Section 2, a realization of finite degree Volterra series in finite dimension state space representation is defined. In Section 3, the output error identification problem is formulated and the structure's parameters are chosen. A local parameterization of the realization of finite degree Volterra series is developed in Section 4.

In Section 5, we propose a sequential projection method to calculate an initial estimation of the structure's parameters.

In Section 6, we summarize the algorithm of identification.

Finally, Section 7 presents some illustrative examples and a comparison with some system identification methods for nonlinear systems.

1 See Appendix.A for the details of calculation.

General notations

As a general rule in this paper, selecting elements of matrices is done using MATLAB TM standard matrix operations, e.g. M(:, i : j) stands for the sub-matrix of the matrix M which contains the columns from the i th to j th columns.

2 Realization of a finite degree Volterra series [START_REF] Brockett | Volterra series and geometric control theory[END_REF] has proved that any observable realization of a finite Volterra series of degree l can be approximated by an recursive realization of the form

Z 1 t = A 1 Z 1 t-1 + B 1 u t , Z 1 0 = 0 Z 2 t = A 2 Z 2 t-1 + B 2 (u t ⊗ Z 1 t ) , Z 2 0 = 0 . . . Z l t = A l Z l t-1 + B l (u t ⊗ Z l-1 t ) , Z l 0 = 0 y t = l ∑ i=1 C i Z i t + v t (4)
The direct relation between y t and u t , can be obtained by using the property FG ⊗ HJ = (F ⊗ H) (G ⊗ J) and a simple development of (4) leads to

y t = l ∑ i=1 t-1 ∑ τ 1 =0 τ 1 ∑ τ 2 =0
. . .

τ i-1 ∑ τ i =0 C i Φ i (τ 1 , . . . , τ i ) ×u t-τ 1 ⊗ u t-τ 2 . . . ⊗ u t-τ i + v t (5)
where

Φ 1 (τ 1 ) = (A 1 ) τ 1 B 1 and for i ≥ 2: Φ i (τ 1 , τ 2 , . . . , τ i ) = (A i ) τ 1 B i I m ⊗ Φ i-1 (τ 2 , . . . , τ i ) (6)
By comparing (1) and ( 5), we observe that Realization (4) leads to Volterra series similar to (1), in which we suppose that u t = 0 for t ≤ 0, and the Volterra kernels are approximated by a sum of matrix products. In fact, this approximation depends on the dimensions of the states {Z i t : i = 1, 2, • • • , l}, and it is clear that it becomes better if we increase the dimensions of states. On the other hand, for practical purposes these dimensions should be chosen as lower as possible in order to find a compromise between the approximation and numerical complexity issues. Obviously, Realization (4) is not unique, assume that the states of Structure (4) transform into X i t = T i -1 Z i t where T i ∈ R n i ×n i nonsingular matrices, and n i is the dimension of the state Z i t . The structure becomes

X 1 t = Ā1 X 1 t-1 + B1 u t X 2 t = Ā2 X 2 t-1 + B2 (u t ⊗ X 1 t ) . . . X l t = Āl X l t-1 + Bl (u t ⊗ X l-1 t ) y t = l ∑ i=1 Ci X i t + v t where Āi Bi Ci 0 = T i -1 A i T i T i -1 B i T i-1 m C i T i 0 (7)
and {T j m : j = 0, 1, 2, • • • , l -1} are defined as follows

T j m = I m ⊗ T j , T 0 m = I m : I m is the identity matrix Lemma 1
The realization of finite degree Volterra series ( 4) is asymptotically stable if and only if:

For i = 1, 2, • • • , l : ρ(A i ) < 1 (8)
where ρ(A i ) stands for the spectral radius of A i .

PROOF. Each subsystem may be viewed as a linear system with input (u t ⊗ Z i-1 t ). So, Z i t will be asymptotically stable if and only if the linear system is stable ρ(A i ) < 1 and

u t ⊗ Z i-1 t is bounded. Assuming that Z 1 t is stable ρ(A 1 ) < 1 , so u t ⊗ Z 1
t is bounded and Z 2 t is stable if and only if ρ(A 2 ) < 1 . And so on.

Output error identification

Our goal is to determine the coefficient matrices of Structure (4). Assume that all matrices are fully parameterized, so the structure's parameters can be given by

θ =                vec(A 1 ) vec(B 1 ) vec(C 1 ) . . . vec(A l ) vec(B l ) vec(C l )                (9)
where vec(.) denotes the vectorization operator defined as follows

vec : M ∈ R m×n → R m•n vec(M) = vec m 1 m 2 • • • m n = m T 1 m T 2 • • • m T n T
Given the input u t and output y t of the real system, the model corresponding to θ can be given as follows

Ẑ1 t = A 1 (θ ) Ẑ1 t-1 + B 1 (θ )u t Ẑ2 t = A 2 (θ ) Ẑ2 t-1 + B 2 (θ )(u t ⊗ Ẑ1 t ) . . . Ẑl t = A l (θ ) Ẑl t-1 + B l (θ )(u t ⊗ Ẑl-1 t ) ŷt (θ ) = l ∑ i=1 C i (θ ) Ẑi t ( 10
)
Note that as Z i t depends on θ , the mapping ŷt (θ ) is nonlinear with θ . Our goal is achieved if the output ŷt (θ ) approximates the output of the real system accurately. This criterion can be transformed into the minimization of the output error with respect to the parameters θ . Considering a data of length N, the output-error cost function is given by

J N (θ ) = 1 N N ∑ k=1 y k -ŷk (θ ) 2 2 = 1 N E N (θ ) T E N (θ ) (11)
where

E N (θ ) = e(1) T e(2) T • • • e(N) T T
(12) is the error vector in which e(k) = y kŷk (θ ). The minimization of ( 11) is clearly a nonlinear, nonconvex optimization problem. The numerical solution of this problem can be calculated by different algorithms, e.g. gradient search method (Levenberg-Marquard method) is a popular one. This iterative method is based on the updating of the system parameters θ as follows

θ i+1 = θ i -(ψ T N (θ i )ψ N (θ i ) + λ i+1 I) -1 ψ T N (θ i )E N (θ i ) (13)
Where λ i is the regularization parameter and

ψ N (θ ) ∂ E N (θ ) ∂ θ T (14)
is the jacobian of the error vector E N (θ ). As we mentioned in Section (2), the structure (4) is not unique. As a consequence, the minimization of J N (θ ) does not have a unique solution.

The nonuniqueness solution of θ is the consequence of the full parameterization of the matrices of the realization. However, the optimal solution can be made unique by choosing a suitable canonical parameterization. As each subsystem of the realization can be considered as a linear system, one could use a classical parameterization of the various parameterization proposed for the linear systems. Unfortunately, these parameterizations are not numerically robust [START_REF] Mckelvey | System identification using an over-parametrized model class-improving the optimization algorithm[END_REF].

To overcome the nonuniqueness problem of the optimal θ and keep the full parameterization of the matrices, Ribarits et al [START_REF] Ribarits | An analysis of the parametrization by data driven local coordinates for multivariable linear systems[END_REF] have proposed a method for the linear systems, in which the directions that do not change the cost of output error function are identified and projected out at each iteration.

For that only the active parameters are updated. In analogues way, we will define a local parameterization of the realization of finite degree Volterra series (4) in order to define the directions in which the cost function J N (θ ) does not change.

Local parameterization

Recall that two realizations of a finite degree Volterra series ( 4) are similar if their coefficient matrices are related by Equation ( 7), where the transformation matrices

{T i : i = 1, 2, • • • , l} parameterize the subset of equivalent models.
The obtained subset defines a manifold. In order to identify the tangent plane of this manifold at

(A i , B i ,C i : i = 1, 2, • • • , l),
we linearize Relation (7) around the identity matrices. Considering a small perturbation T i = I n i + ∆T i , we suppose that ρ(∆T i ) ≪ 1 then by using the approximation

I n i + ∆T i -1 ≃ I n i -∆T i and neglecting all second order terms, we obtain Āi Bi Ci 0 = A i B i C i 0 + -∆T i A i + A i ∆T i -∆T i B i C i ∆T i 0 + 0 B i ∆T i-1 m 0 0 (15)
By defining

Λ j i 0 n j ×(i-1)n j I n j 0 n j ×(m-i)n j
it is possible to write ∆T j m as follows

∆T j m = m ∑ i=1 Λ j i T ∆T j Λ j i
If we consider the following vectors of parameters (16) the relation between θ and θ can be obtained using the property vec(ABC) = C T ⊗ A vec(B)

θ =                vec(A 1 ) vec(B 1 ) vec(C 1 ) . . . vec(A l ) vec(B l ) vec(C l )                and θ =                vec( Ā1 ) vec( B1 ) vec( C1 ) . . . vec( Āl ) vec( Bl ) vec( Cl )               
θ = θ + M θ     vec(∆T 1 ) . . . vec(∆T l )     (17)
where for 1 ≤ j ≤ l -1

M θ :, 1 + j-1 ∑ k=1 n 2 k : j ∑ k=1 n 2 k =                  0 α j ×n 2 j -A j T ⊗ I n j + I n j ⊗ A j -B j T ⊗ I n j I n j ⊗C j 0 n 2 j+1 ×n 2 j m ∑ i=1 Λ j i T ⊗ B j+1 Λ j i T 0 γ j ×n 2 j                 
and

M θ :, 1 + l-1 ∑ k=1 n 2 k : l ∑ k=1 n 2 k =        0 α l ×n 2 l -A l T ⊗ I n l + I n l ⊗ A l -B l T ⊗ I n l I n l ⊗C l        (18) with α j = ∑ j-1 i=1 n i (n i + m n i-1 + p) γ j = p n j+1 + ∑ l i= j+2 n i (n i + m n i-1 + p) Lemma 2
The left null space of M θ (18) defines a basis of the directions in which the parameters should be modified to lead a change in the value of cost function J N (θ ).

PROOF. Equation ( 17) shows that, the tangent space of the manifold of equivalent realizations at

(A i , B i ,C i : i = 1, 2, • • • , l
) is equal to the column space of the matrix M θ (18). Since the left null space of the matrix M θ is orthogonal complement to the column space, the directions in which the value of the cost function J N (θ ) changes are those related to left null space of M θ .

Let the QR decomposition of M θ be given by

M θ = Q 1 Q 2 R 0 (19)
then a basis of the null space of M θ is Q 2 . Thus, the update rule should be modified such that we project out the directions in which the cost function does not change.

The new update rule becomes

θ i = θ i-1 -Q 2 Q T 2 ψ T N ψ N Q 2 + λ i I -1 Q T 2 ψ T N E N ( 20
)
where Q 2 , ψ N and E N depend on θ i-1 . Note that since Q 2 depends on the past parameter θ i-1 , the QR decomposition (19) must be computed at each iteration.

Computing the iterative parameter update

In order to compute the update rule (20), the quantities E N (θ ) and ψ N (θ ) should be computed. Computing the vector E N (θ ) can be done by simulating the system (10) with θ = θ i-1 . At the same time, this simulation brings out the states { Ẑi

t : i = 1, 2, • • • , l}.
In order to simulate the gradient ψ N (θ i-1 ), we should compute the derivative of ŷt with respect to θ i-1 . Let us define

ζ j t,k = ∂ Ẑ j t ∂ θ k (21) then the computation of ∂ ŷt ∂ θ T = ∂ ŷt ∂ θ 1 • • • ∂ ŷt ∂ θ q
, where q is the number of parameters in θ , can be made as follows

ζ 1 t,k =A 1 ζ 1 t-1,k + ∂ A 1 ∂ θ k Ẑ1 t-1 + ∂ B 1 ∂ θ k u t ζ 2 t,k =A 2 ζ 2 t-1,k + ∂ A 2 ∂ θ k Ẑ2 t-1 + ∂ B 2 ∂ θ k (u t ⊗ Ẑ1 t ) + B 2 (u t ⊗ ζ 1 t,k ) . . . ζ l t,k =A l ζ l t-1,k + ∂ A l ∂ θ k Ẑl t-1 + ∂ B l ∂ θ k (u t ⊗ Ẑl-1 t ) + B l (u t ⊗ ζ l-1 t,k ) ∂ ŷt ∂ θ k = l ∑ j=1 C j ζ j t,k + ∂C j ∂ θ k Ẑ j t (22)
In fact, if we consider the local parameterization of finite degree Volterra realization, then computing the update rule (20) can be done without calculate ψ N first. Let us define

Ω N ψ N Q 2 = ∂ E N ∂ θ T Q 2 (23)
So, the update rule can be expressed as follows

θ i = θ i-1 -Q 2 Ω T N Ω N + λ i I -1 Ω T N E N (24)
Consider the t p elements of the s th column of

Ω N Ω N ((t -1) p + 1 : t p, s) = q ∑ k=1 ∂ ŷt ∂ θ k Q 2 (k, s) (25)
In order to calculate this sum, let us define

ζ j t = q ∑ k=1 ∂ Ẑ j t ∂ θ k Q 2 (k, s) (26)
Using ( 22), the computation of ( 25) can be done as follows

ζ 1 t =A 1 ζ 1 t-1 + ∆A 1 Ẑ1 t-1 + ∆B 1 u t ζ 2 t =A 2 ζ 2 t-1 + ∆A 2 Ẑ2 t-1 + ∆B 2 (u t ⊗ Ẑ1 t ) + B 2 (u t ⊗ ζ 1 t ) . . . ζ l t =A l ζ l t-1 + ∆A l Ẑl t-1 + ∆B l (u t ⊗ Ẑl-1 t ) + B l (u t ⊗ ζ l-1 t ) q ∑ k=1 ∂ ŷt ∂ θ k Q 2 (k, s) = l ∑ j=1 C j ζ j t + ∆C j Ẑ j t ( 27
) where

∆A j = q ∑ k=1 ∂ A j ∂ θ k Q 2 (k, s) : j = 1, 2, • • • , l (28)
and ∆B j , ∆C j are defined analogously. These matrices can be obtained from

               vec(∆A 1 ) vec(∆B 1 ) vec(∆C 1 ) . . . vec(∆A l ) vec(∆B l ) vec(∆C l )                = q ∑ k=1 ∂ θ ∂ θ k Q 2 (k, s) = q ∑ k=1 e k Q 2 (k, s) (29)
where

e k = [0 . . . 0 1 ↑ k 0 . . . 0] T
Note that the number of columns of Ω N is smaller than that of ψ N because the first one is the number of active parameters and the second one is the total number of parameters.

Computing an initial estimation

The reason why a good initial estimation is crucial is that the local gradient search method converges to the nearest local optimum in the neighborhood of initial guess of solution. Moreover, the initial estimation should also provide a stable realization that verifies the constrain (8).

In this section, we propose a sequential projection method in order to provide an initial estimation of the realization's vector of parameters (θ ). The basic idea is to use the consecutive property of finite degree Volterra series realization to define a sequential projection procedure. Recall that, the error introduced in the projection process is due to the non orthogonality of the subprocesses (linear, quadratic, cubic, • • • ), and as well the effect of the noise on the observed output signal. The algorithm of projection can be resumed as follows

(1) Estimate the best linear stable approximation (the linear subsystem) of the nonlinear system. (2) From the simulation of the linear subsystem, we project the residual on the class of 2-degree subsystems.

(3) From the simulation of linear + 2-degree subsystems, we project the residual on the class of 3-degree subsystem. (4) Etc.

Moreover, this projection procedure yields an estimation of the dimensions of the states {Z i t : i = 1, 2, • • • , l}, which means the set {n i : i = 1, 2, • • • , l}. Furthermore, it yields the degree of the realization (l) by defining a precision criterion as it will appear in the sequel.

Identification of linear subsystem

Estimating the best linear stable approximation of nonlinear system can be formulated as a minimization problem

min L 1 N = 1 N N ∑ t=1 y t -y 1 t 2 2
where

Z 1 t = A 1 Z 1 t-1 + B 1 u t y 1 t = C 1 Z 1 t with the constraint ρ A 1 < 1 (30)
The parameters of the optimization problem, which should be estimated, are the triple (A 1 , B 1 ,C 1 ). The above minimization problem can be transformed into the following matricial form [START_REF] Gopinath | On the identification of linear timeinvariant systems from input-ouput data[END_REF][START_REF] Demoor | A geometrical strategy for the identication of state space models of linear multivariable systems a geometrical strategy for the identication of state space models of linear multivariable system with singular value decomposition[END_REF])

min Y α,N -Γ 1 α Z 1 0,N-α -Φ 1 α U α,N 2 F subject to ρ A 1 < 1 (31)
where . F denotes Frobenius norm. The input and output are stocked in Hankel matrices form

U α,N        u 1 u 2 • • • u N-α+1 u 2 u 3 • • • u N-α+2 . . . . . . . . . . . . u α u α+1 • • • u N       
where α and N refer to the number of rows in the matrix and data length respectively. The output Hankel matrix Y α,N is defined analogously to U α,N . The matrices Γ 1 α , Φ 1 α and Z 1 0,N-α are defined as follows

Γ 1 α =        C 1 C 1 A 1 . . . C 1 (A 1 ) (α-1)        Φ 1 α =        C 1 B 1 0 0 • • • 0 C 1 A 1 B 1 C 1 B 1 0 • • • 0 . . . . . . . . . . . . . . . C 1 (A 1 ) (α-1) B 1 • • • • • • • • • C 1 B 1        Z 1 0,N-α = Z 1 0 Z 1 1 • • • Z 1 N-α .
Note that the number of rows (α) should be roughly chosen to be greater than the expected linear system order n 1 . This condition guarantees that the extended matrix of observability Γ 1 α has a full rank. Our objective is to estimate the matrices Γ 1 α and Φ 1 α . Since only the data matrices U α,N and Y α,N are known, instead of solving the minimization problem (30) with respect to Γ 1 α and Φ 1 α , we can transform it into an equivalent one involving only the matrix Γ 1 α as follows

min Y α,N Π ⊥ U T α,N -Γ 1 α Z 1 0,N-α Π ⊥ U T α,N 2 F subject to ρ A 1 < 1 (32)
where

Π ⊥ U T α,N
is the orthogonal projection onto the nullspace

of U α,N Π ⊥ U T α,N = I -U T α,N U α,N U T α,N -1 U α,N (33)
such that U α,N Π ⊥ U T α,N = 0 (34)
The inverse of the matrix U α,N U T α,N exists if the input is persistently exciting and N > mα. In fact, the origin of the idea of subtracting the term that involves the input signal belongs to the direct 4SID method, which is a classical subspace method [START_REF] Demoor | A geometrical strategy for the identication of state space models of linear multivariable systems a geometrical strategy for the identication of state space models of linear multivariable system with singular value decomposition[END_REF]. To solve the problem (32), one could use the Singular Value Decomposition (SVD). Recall that, the extended observability matrix (Γ 1 α ) is a full rank matrix. Consider the following SVD

Y α,N Π ⊥ U T α,N = Q s Q n S s 0 0 S n V T s V T n ( 35
)
where the matrix S s contains the principals singular values (further than threshold).

The dimension of this matrix yields n 1 (the dimension of Z 1 t ). The estimation of the matrix Γ 1 α can be done by solving the following problem of minimization

min Γ 1 α -Q s 2 F subject to ρ A 1 < 1 (36)
By using the property

Γ 1 α (1 : (α -1)p, :)A 1 =Γ 1 α (p + 1 : α p, :) (37)
and replacing the stability condition ρ(A 1 ) < 1 by equivalent Lyapunov inequalities

ρ(A 1 ) < 1 ⇐⇒ ∃P ≥ δ I n 1 : P -A 1 PA 1 T ≥ δ I n 1
where δ > 0, we can transform the minimization problem (36) into the following one

min J A 1 L A 1 Q † s -Q ‡ s A 1 R A 1 2 F subject to P -δ I n 1 X T X P ≥ 0 2n 1 (38) where Q † s Q s (p + 1 : α p, :) Q ‡ s Q s (1 : (α -1) p, :) X A 1 P (39)
If we let R A 1 P and L A 1 be equal to the identity matrix, the problem is converted to an optimization problem that involves minimizing a linear function over symmetric cones and can be solved using SeDuMi MATLAB TM package [START_REF] Sturm | Using SeDuMi 1.0x, a MATLAB Toolbox for Optimization Over Symmetric Cones[END_REF] (see [START_REF] Lacy | Subspace identification with guaranteed stability using constrained optimization[END_REF] for more details).

Solving the minimization problem ( 38) provides an estimation of a stable matrix A 1 . The matrix C 1 can be estimated using the following formula

C 1 = α-1 ∑ k=0 Q s (kp + 1 : (k + 1)p, :) Λ T A 1 Λ A 1 Λ T A 1 -1 (40) where Λ A 1 I n 1 + ∑ α-1
i=1 A 1 i . Note that as A 1 verifies that ρ(A) < 1, Λ A 1 is a full rank matrix and the matrix inverse exists. The matrix Φ 1 α can be estimated by considering the least-squares solution to the overdetermined system of equations

Q T n Y α,N U T α,N U α,N U T α,N -1 = Q T n Φ 1 α (41)
Finally, the matrix B 1 can be easily calculated from the estimation of Φ 1 α .

Identification of higher order subsystems

After using the procedure described in the previous section, an estimation of (A 1 , B 1 ,C 1 ) is available. Calculating the state Z 1 t and the output y 1 t for t = 1, .., N of the linear subsystem can be done by evaluating the following model

Z 1 t = A 1 Z 1 t-1 + B 1 u t y 1 t = C 1 Z 1 t
The next task is to project the residual on the class of 2degree subsystem. We define the residual as follows ỹt = y ty 1 t (42)

The estimation of the best 2-degree subsystem can be done by solving the following optimization problem

min L 2 N = 1 N N ∑ t=1 ỹt -y 2 t 2 2
where

Z 2 t = A 2 Z 2 t-1 + B 2 u t ⊗ Z 1 t y 2 t = C 2 Z 2 t with the constraint ρ A 2 < 1 (43)
As u t ⊗ Z 1 t is computed thanks to the previous step, this minimization problem is similar to (30). So, by using an analogous logic, we obtain an estimation of the matrices

A 2 , B 2 ,C 2 .
Then, this procedure is reiterated until we get an appreciated precision. Defining this precision criterion is the objective of the following section. It is clear from Fig. 1 that using local parameterization reduces the time of computation. This reducing becomes advantageous when the realization degree increases. In order to corroborate that the proposed method is able to identify the finite Volterra series efficiently, we calculate the average and the standard deviation of the accuracies of the obtained models. The accuracies are reported in Table 1 where the initial estimated model is the model obtained by the sequential projection explained in Section 5 and the optimized model is that obtained by local gradient search in local parameterization space. These accuracies demonstrate that the proposed method has efficiently identified the finite degree Volterra realization and the sequential projection method provides a good initial guess for the optimization process.

First example of nonlinear system

Consider the following example

x 1 t =0.7 x 1 t-1 + u t , x 1 0 = 0 x 2 t =0.7 x 2 t-1 + x 1 t 2 + x 1 t + u t , x 2 0 = 0 y t = x 1 t 2 + x 2 t + v t (47)
The input signal u t ∈ R is chosen to be uniform white noise with standard deviation is equal to 0.5 and zero mean, its length is equal to 1000 samples. The measurement noise v t is a Gaussian white noise scaled such that, the signal to noise ratio SNR = 10 dB. We compare the performance of finite Volterra series realization with neural network based nonlinear system identification, which is a popular approach [START_REF] Billings | Properties of neural networks with applications to modelling nonlinear dynamical systems[END_REF][START_REF] Sjoberg | Neural networks in system identification[END_REF][START_REF] Nørgaard | Neural networks for Modelling and Control of Dynamic Systems[END_REF]. For comparison, we have chosen two nonlinear model structures

(1) Neural Network Output Error (NNOE) model with 10 hidden hyperbolic tangent units in the hidden layer.

The output function generated by the neural network can be calculated as follows

ϕ t = ŷt-1 (θ ) ŷt-2 (θ ) ŷt-3 (θ ) u t u t-1 u t-2 T ŷt (θ ) = g (ϕ t , θ ) (48
) where ϕ t is a vector containing the regressors, θ is a vector containing the weights and g is the function realized by the neural network.

(2) Neural Network State Space Innovation Form (NNSSIF) model which determines a nonlinear state space model of the dynamic system

xt (θ ) = g ( xt-1 (θ ), u t ) ŷt = C xt (θ ) (49)
where xt (θ ) is the state vector and its dimension is n x , g is the function realized by the neural network and C is a constant matrix.

In this example, the best performance was obtained by n x = 6 and neural network with 10 hidden hyperbolic tangent units in the hidden layer. The neural network is trained with the Levenberg-Marquardt method.

The implementation of two models (NNOE and NNSSIF) is done using Neural Network Based System Identification toolbox [START_REF] Nørgaard | Neural Network Based System Identification Toolbox[END_REF].

Note that the dimensions of the states (Z i t : i = 1, 2) of the quadratic system (Volterra series realization of degree two) are 2 and 3 respectively. Recall that these dimensions are obtained using the sequential projection method explained in Section 5. The accuracies reported in Table 2 show the outperformance of a quadratic system in comparison with NNOE and NNSSIF models.

To deal with a clear graphical representation, we represent the results in the validation interval over a window of 150 samples in Fig. 3.

Second example of nonlinear system

Consider the following example of nonlinear system

x 1 t =0.8x 1 t-1 + 0.8u t (1) + 0.6u t (2) + 0.8u t (3) x 2 t =0.5 x 1 t-1 2 + 0.5x 2 t-1 + 0.9 (u t (2)) 3 + u t (1)u t (2)x 1 t-1 x 3 t =0.7x 2 t-1 + 0.7x 3 t-1 + 0.8u t (2)u t (3)x 1 t-1 + 0.9 (u t (1)) 3 • • • + 0.9 (u t (2)) 3 y t =    
0.6 0.3 0.6 0.2 0.7 0.6 0.3 0.0 0.5

        x 1 t x 2 t x 3 t     + v t (50)
where the non observed nonlinear state is

x t = x 1 t , x 2 t , x 3 t T ,
and u t (i) denotes the i th element of the input signal u t ∈ R 3 .

The system is simulated with three dimensional uniform white noise with standard deviation is equal to 1 as input and with the initial condition x 1 0 , x 2 0 , x 3 0 T = 0, 0, 0 T . The length of the input data is equal to 2000 samples. The measurement noise v t is a Gaussian white noise is scaled such that the SNR = 15 dB. We compare the performance of finite degree Volterra series realization with NNSSIF model which can handle Multiple Inputs -Multiple Outputs (MIMO) systems. The best performance was obtained by n x = 10 and 15 hidden hyperbolic tangent units in the hidden layer.

Note that the dimensions of the states (Z i t : i = 1, 2, 3) of the cubic system (Volterra series realization of degree three) are 6, 9 and 13 respectively. The accuracies reported in Table 3 show the outperformance of a cubic system in comparison with NNSSIF model. To deal with a clear graphical representation, we represent the results in the validation interval over a window of 150 samples in Fig. 4. It well be of interest to compare the performance of finite Volterra series realization with truncated Volterra series which has the same degree. Therefore, the accuracies and the number of parameters of truncated Volterra series of degree three with various horizon lengths T are calculated and reported in Table 4. Note that the accuracy reported in Table 4 is the total accuracy of the output signal y t ∈ R 3 . From Table 4, three remarks can be concluded

(1) The number of parameters of cubic system is much smaller than that of truncated Volterra series of degree three .

(2) The validation accuracy of cubic system is much higher than that of truncated Volterra series of degree three. This result is particularly interesting. The reason why the prediction of cubic system is more accurate than truncated Volterra series is that the cubic system consider the integral past of the input signal and its parameters are optimized such that the output of the model accurately approximates the real output of system. On the contrary, the truncated Volterra series consider a fixed horizon and calculate the Volterra kernels which correspond to this horizon.

(3) On the contrary of theoretical expectations, the validation accuracy of truncated Volterra series decreases with increasing the horizon T . This because the number of parameters increases and as well the size of matrices which become ill-conditioned.

Conclusion

In this paper, we have presented a new method to identify a recursive state-space realization of finite degree Volterra series with infinite horizon. The method is based on a local parameterization of the state-space representation of the realization and subsequent gradient search in the resulting local parameter space. Furthermore, we have proposed a sequential projection procedure to calculate an initial estimation of the realization's parameters.

The method has successfully applied to identify various illustrative examples, and a comparison with some methods from the state-of-the-art nonlinear system identification methods have pointed out the outperformance of finite Volterra series realization. Moreover, a comparison with truncated Volterra series of the same degree has borne out that not only the number of parameters of finite Volterra series realization is much smaller than that of truncated Volterra series, but also the prediction accuracy of finite Volterra series realization is superior to that of truncated Volterra series.

Our future work will focus on the improvement of the optimization algorithm in order to guarantee the convergence to the global optimum. Moreover, we will work on the development of control techniques for this realization in order to control real world nonlinear systems.

A Appendix

Consider that the total number of parameters of a truncated Volterra series of degree l and horizon length T is C (l, T ).

Let us define C n (T ) as the total number of parameters of the Volterra kernels of degree n ( the number of elements of the matrices w n (τ 1 , τ 2 , . . . ,

τ n ), 0 ≤ τ n ≤ τ n-1 ≤ • • • ≤ τ 1 ≤ T ).
It is obvious that

C (l, T ) = l ∑ n=1 C n (T ) (A.1)
Suppose S is any finite set of elements and denote the number of these elements ♯S, then

C n (T ) = ♯ {0 ≤ τ n ≤ τ n-1 ≤ • • • ≤ τ 1 ≤ T } p m n = T ∑ k=0 ♯ {0 ≤ τ n ≤ τ n-1 ≤ • • • ≤ τ 1 ≤ k} p m n = T ∑ k=0 k + n -1 n -1 p m n

Fig. 1 .Fig. 1 .

 11 Fig.1. The implementation of two approaches is done using MATLAB TM programming environment.

Fig. 2 .

 2 Fig. 2. Computation time of the sequential projection method.

  The real output without the measurement noise (solid line) with the outputs of quadratic system (dashed line) and NNOE model (dashdotted line), and their errors.

  The real output without the measurement noise (solid line) with the outputs of quadratic system (dashed line) and NNSSIF model (dashdotted line), and their errors.

Fig. 3 .

 3 Fig. 3. Comparison of the results of quadratic system, NNOE model, and NNSSIF model.

Fig. 4 .

 4 Fig.4. The real outputs without the measurement noise (solid line) with the outputs of cubic system (dashed line) and NNSSIF model (dash-dotted line), and their errors are superimposed.

Table 1

 1 Accuracy of the identified models

	Identification	Initial estimated Optimized
		model	model
	Accuracy (%VAF)	82.3 ± 2.7	89.1± 1.8
	Validation
	Accuracy (%VAF)	86.9 ± 1.3	98.6 ± 0.7

Table 2

 2 Accuracy of the identified models

	Identification	Quadratic NNOE NNSSIF
		system	model	model
	Accuracy (%VAF)	85.7	76.5	63.9
	Validation
	Accuracy (%VAF)	98.0	81.9	73.4

Table 3

 3 Accuracy of the identified models

	Identification	Accuracy (%VAF)
		Cubic system NNSSIF model
	First output	96.8	79.2
	Second output	96.9	82.4
	Third output	96.0	86.5
	Validation
	First output	98.0	82.3
	Second output	97.6	83.8
	Third output	97.3	88.5

Determining the realization's degree

Consider that, the subsystems of order inferior to j are available. In order to verify that the nonlinear system is accurately fitted by the finite degree Volterra series realization, we evaluate the following criterion ŷt

where η is a user-defined constant and ȳ denotes the mean of the output signals defined as follows

If the criterion ( 44) is true, the algorithm stops and the realization's degree is j. Otherwise we estimate the subsystem of order j + 1. Note that the measurement noise v t should be considered when we define η as follows

× 100 (46)

Identification algorithm

We can resume the algorithm of identification as follows

(1) Calculate an initial estimation of the realization's parameters using the sequential projection described in Section 5. Recall that, this procedure yields the realization degree (l), the dimensions of the states

The estimated vector of parameters θ 0 is used as an initial guess for the optimization process and k ← 0.

(2) Calculate the states {Z i t : i = 1, 2, • • • , l} and ŷt by simulating the system (10) with θ = θ k .

(3) Compute E N (θ ) using ( 12). (4) Calculate the matrix M θ using (18).

(5) Calculate the QR decomposition of M θ (19), from which we obtain Q 2 . (6) Calculate ∆A j , ∆B j and ∆C j : j = 1, 2, • • • , l using (29). ( 7) Calculate the matrix Ω N using ( 25) and ( 27), we suppose that {ζ j 0 = 0 : j = 1, 2, • • • , l}. (8) Calculate the update rule of the gradient search algorithm using (24) and k ← k + 1.

(9) Perform the termination test for minimization, If true, the algorithm stops. Otherwise, return to step (2), i.e. compute the values J N θ k-1 and J N θ k using (11) and test if J N θ k -J N θ k-1 2 is small enough.

Illustrative examples

In this section, first we compare the computational complexity of the direct gradient search and gradient search in the local parameterization of finite degree Volterra realization as a function of realization degree. At the same time, we corroborate that the method is able to identify the realization of finite degree Volterra series. Second, we consider two examples of nonlinear system, and we compare the capability of finite degree Volterra realization to approximate the considered examples with some of state-of-the-art system identification methods for nonlinear systems.

In order to validate the obtained models, for a data of N samples, we have used T id = N 2 samples for the identification purpose, and the rest of them T val = T -T id = N 2 samples for the validation purpose. The validation step, can be viewed as an evaluation of prediction accuracy of the models, and in order to verify that the models are capable to extract the real output signal from the measured noised one, we consider the output signal without the measurement noise in the validation step. The model accuracy is defined as the Percent Variance Accounted For (%VAF)

where ŷt denotes the estimated output signal and ȳ is the mean of the output signals.

Computational complexity

In this section, we figure out that using a local parameterization accelerates the convergence of gradient search method.

For that, we consider a finite degree Volterra series realization (4), where the states {Z i t : i = 1, 2, • • • , l} have the same dimension which is equal to n = 10. The identification experiment for each value of the realization's degree is repeated 10 times. Each one has a length of 4000 samples. The average of optimization time is then computed. The optimization time is the required time to reach a specific precision (local optimum). Each experiment has three inputs ( u t ∈ R 3 ) which have been chosen to be uniform white noise and three outputs ( y t ∈ R 3 ) . The measurement noise v t is a Gaussian white noise scaled such that the signal to noise ratio SNR = 10 dB. The initial estimation of the realization's vector of parameters (θ ) is computed using the sequential projection method explained in Section 5. The computation time of sequential projection method is calculated and given in Fig. 2. The optimization time as a function of realization degree (l) for the direct gradient search method and the gradient search in local parameterization space is given in


Mouhamed Ndong
email: mouhamed.ndong@polymtl.ca

David Bird

Tri Nguyen-Quang

Marie-Laure De Boutray

Arash Zamyadi

Brigitte Vinçon-Leite

Bruno J Lemaire

Michèle Prévost

Sarah Dorner

Estimating the risk of cyanobacterial occurrence using an index integrating meteorological factors: application to drinking water production

Keywords: cyanobacteria, meteorological index, drinking water intake, early warning systems, wind, nutrient mass ratios

The sudden appearance of toxic cyanobacteria (CB) blooms is still largely unpredictable in waters worldwide. Many post-hoc explanations for CB bloom occurrence relating to physical and biochemical conditions in lakes have been developed. As potentially toxic CB can accumulate in drinking water treatment plants and disrupt water treatment, there is a need for water treatment operators to determine whether conditions are favourable for the proliferation and accumulation of CB in source waters in order to adjust drinking water treatment accordingly. Thus, a new methodology with locally adaptable variables is proposed in order to have a single index, f(p), related to various environmental factors such as temperature, wind speed and direction. The index is used in conjunction with real time monitoring data to determine the probability of CB occurrence in relation to meteorological factors, and was tested at a drinking water intake in Missisquoi Bay, a shallow transboundary bay in Lake Champlain, Québec, Canada. These environmental factors alone were able to explain a maximum probability of 68% that a CB bloom would occur at the drinking water treatment plant. Nutrient limitation also influences CB blooms and intense blooms only occurred when the dissolved inorganic nitrogen (DIN) to total phosphorus (TP) mass ratio was below 3. Additional monitoring of DIN and TP could be considered for these source waters prone to cyanobacterial blooms to determine periods of favourable growth. Real time monitoring and the use of the index could permit an adequate and timely response to CB blooms in drinking water sources.

Introduction

Harmful cyanobacteria (CB) blooms have been noted worldwide. CB blooms may release toxins, which can be harmful for aquatic communities [START_REF] Bartram | Toxic cyanobacteria in Water: A guide to their public health consequences, monitoring and management[END_REF], limit recreational and economic activities [START_REF] Steffensen | Cyanobacterial Harmful Algal Blooms: State of the Science and Research Needs[END_REF] and constitute a threat to drinking water sources (Zamyadi et al., 2012a). Long-term solutions for restoring water bodies require an understanding of the factors that are contributing to the increase of CB blooms [START_REF] Wongsai | Understanding environmental factors associated with cyanobacterial bloom[END_REF]. Thus, considerable efforts have been made globally in order to understand the environmental factors related to these phenomena for better management of the affected water bodies [START_REF] Trojanowska | Phosphorus fractions transformation in sediments before and after cyanobacterial bloom: implications for reduction of eutrophication symptoms in dam reservoir[END_REF]Izydorczyk, 2010, Le Vu et al., 2010;[START_REF] Falconer | Cyanobacterial toxins of drinking water supplies-Cylindrospermopsins and microcystins: Boca Raton, Florida[END_REF]. Both endogenous factors (e.g. nutrients, water temperature) and exogenous factors (air temperature, solar radiation, wind speed and wind direction) can be considered to play a role in CB bloom development [START_REF] Howard | Problem cyanobacteria blooms: Explanation and simulation modeling[END_REF]. Short term fluctuations of meteorological variables play an important role with regards to CB occurrence [START_REF] Wu | Dynamics of cyanobacterial bloom formation during short-term hydrodynamic fluctuation in a large shallow, eutrophic, and wind-exposed Lake Taihu, China[END_REF]. Generally, the principal focus of past studies has been on the effects of nutrients. Thus, the mass ratio (TN/TP) of total nitrogen (TN) to total phosphorus (TP) in water and in biomass has been studied as an index for its relationship to CB blooms [START_REF] Lilover | The variability of parameters controlling the cyanobacteria bloom biomass in the baltic Sea[END_REF][START_REF] Havens | N:P ratios, light limitation, and cyanobacterial dominance in a subtropical lake impacted by non-point source nutrient pollution[END_REF][START_REF] Smith | Light and nutrient effects on the relative biomass of blue green algae in lake phytoplankton[END_REF][START_REF] Schindler | Evolution of phosphorus limitation in lakes[END_REF]. However, others have had concerns [START_REF] Ferber | Do cyanobacteria dominate in eutrophic lakes because they fix atmospheric nitrogen?[END_REF] or do not support the use of TN/TP mass ratio [START_REF] Jayatissa | Occurrence of toxigenic cyanobacterial blooms in freshwaters of Sri Lanka[END_REF]. It has been suggested that eutrophication generated by phosphorus is the main cause of CB blooms [START_REF] Christian | Planktonic bloom-forming cyanobacteria and the eutrophication of lakes and rivers[END_REF], but this view has been qualified of overly simplistic and is not applicable to all water bodies (e.g. [START_REF] Vincent | Rapid physiological assays for nutrient demand by the plankton. I. Nitrogen[END_REF]. Thus, the TN/TP mass ratio is not the only indicator to predict CB blooms and a consensus does not exist with regards to a quantitative value. Another approach related to the mass ratio of Dissolved Inorganic Nitrogen (DIN) to total phosphorus (TP) has also been proposed and was shown to be a better predictor of CB biomass than the TN/TP mass ratio with CB bloom occurrences observed when the DIN/TP mass ratio was below 2 [START_REF] Bergstrom | The use of TN:TP and DIN:TP ratios as indicators for phytoplankton nutrient limitation in oligotrophic lakes affected by N deposition[END_REF][START_REF] Ptacnik | Performance of the Redfield Ratio and a Family of Nutrient Limitation Indicators as Thresholds for Phytoplankton N vs. P Limitation[END_REF].

In addition to nutrients, meteorological conditions are required to explain CB bloom occurrences. Meteorological variables are known to have a significant influence on CB development [START_REF] Elliott | Is the future blue-green? A review of the current model predictions of how climate change could affect pelagic freshwater cyanobacteria[END_REF][START_REF] Reichwaldt | Effects of rainfall patterns on toxic cyanobacterial blooms in a changing climate: Between simplistic scenarios and complex dynamics[END_REF]Liu et al., 2012a). Many authors have explored the relationships among meteorological variables and CB (e.g. [START_REF] Lilover | The variability of parameters controlling the cyanobacteria bloom biomass in the baltic Sea[END_REF][START_REF] Bormans | Spatial and temporal variability in cyanobacterial populations controlled by physical processes[END_REF][START_REF] Kanoshina | The influence of weather conditions (temperature and wind) on cyanobacterial bloom development in the Gulf of Finland (Baltic Sea)[END_REF]Liu et al., 2012a;[START_REF] Liu | Modelling the Effect of Weather Conditions on Cyanobacterial Bloom Outbreaks in Lake Dianchi: a Rough Decision-Adjusted Logistic Regression Model[END_REF]. Various weather variables have been related to CB blooms such as air temperature [START_REF] Elliott | Is the future blue-green? A review of the current model predictions of how climate change could affect pelagic freshwater cyanobacteria[END_REF][START_REF] Reichwaldt | Effects of rainfall patterns on toxic cyanobacterial blooms in a changing climate: Between simplistic scenarios and complex dynamics[END_REF][START_REF] Zhang | Contributions of meteorology to the phenology of cyanobacterial blooms: Implications for future climate change[END_REF], rainfall [START_REF] Reichwaldt | Effects of rainfall patterns on toxic cyanobacterial blooms in a changing climate: Between simplistic scenarios and complex dynamics[END_REF], wind (Liu et al., 2012a;[START_REF] Liu | Modelling the Effect of Weather Conditions on Cyanobacterial Bloom Outbreaks in Lake Dianchi: a Rough Decision-Adjusted Logistic Regression Model[END_REF][START_REF] Zhang | Contributions of meteorology to the phenology of cyanobacterial blooms: Implications for future climate change[END_REF] and sunshine hours [START_REF] Zhang | Contributions of meteorology to the phenology of cyanobacterial blooms: Implications for future climate change[END_REF]. Wind in particular may play both direct and indirect roles in determining cyanobacterial densities [START_REF] Falconer | Safe levels and safe practices[END_REF]. Wind can generate complex physical processes, such as internal waves or seiches, which produce periodic motion of the water layers in stratified lakes (e.g. [START_REF] Lemmin | Test of an extension to internal seiches of Defant's procedure for determination of surface seiche characteristics in real lakes[END_REF]. These internal waves affect the spatial distribution of phytoplankton in the epilimnetic layer [START_REF] Marcé | Interaction between windinduced seiches and convective cooling governs algal distribution in a canyon-shaped reservoir[END_REF]. In deep lakes, they have a great impact on the distribution of cyanobacteria proliferating in the metalimnion and on their growth though a direct influence on light availability [START_REF] Cuypers | Impact of internal waves on the spatial distribution of Planktothrix rubescens (cyanobacteria) in an alpine lake[END_REF]. As CB are buoyant, wind can result in accumulation of cyanobacteria at specific locations during a bloom [START_REF] Cao | Effects of Wind and Wind-Induced Waves on Vertical Phytoplankton Distribution and Surface Blooms of Microcystis aeruginosa in Lake Taihu[END_REF][START_REF] Ostos | The influence of wind-induced mixing on the vertical distribution of buoyant and sinking phytoplankton species[END_REF][START_REF] Kanoshina | The influence of weather conditions (temperature and wind) on cyanobacterial bloom development in the Gulf of Finland (Baltic Sea)[END_REF], but it may also indirectly influence nutrient availability through mixing processes of the water column (e.g. [START_REF] Macintyre | Boundary Mixing and Nutrient Fluxes in Mono Lake, California[END_REF][START_REF] Livingstone | Impact of secular climate change on the thermal structure of a large temperate central European lake[END_REF][START_REF] Schmittner | Decline of the marine ecosystem caused by a reduction in the Atlantic overturning circulation[END_REF].

Other meteorological variables such as the relative humidity, minimum daily air temperature and the amount of solar radiation (sunspots) have also been related to CB density [START_REF] Hu | Weather Variability, Sunspots, and the Blooms of Cyanobacteria[END_REF]. It has been hypothesized that the increase in the number of CB blooms is at least partly related to global climate change and increasing air temperatures [START_REF] Zhang | Contributions of meteorology to the phenology of cyanobacterial blooms: Implications for future climate change[END_REF][START_REF] Wagner | Cyanobacteria dominance: quantifying the effects of climate change[END_REF][START_REF] Mooij | Predicting the effect of climate change on temperate shallow lakes with the ecosystem model PCLake[END_REF]. Although the effect of water temperature on growth is species-specific [START_REF] Mehnert | Competitiveness of invasive and native cyanobacteria from temperate freshwaters under various light and temperature conditions[END_REF]) some generalizations can be made for the CB native to temperate environments. [START_REF] Seaburg | Temperature-growth reponses of algal isolates from Antartic oaes[END_REF] found that the majority of CB have an optimal growth temperature of 25°C. However, many bloom forming CB show optimal growth rates at temperatures of 25°C or higher [START_REF] Robarts | Temperature effects on photosynthetic capacity, respiration, and growth rates of bloom-forming cyanobacteria[END_REF]. [START_REF] Tang | Cyanobacterial dominance of polar freshwater ecosystems: are high-latitude mat-formers adapted to low temperature[END_REF], in studying the predominance of CB in freshwater in polar environments, showed that the growth was undetectable when the temperature of the aquatic environment in which the CB live was less than or equal to 5°C or greater than 35°C.

Consideration of meteorological variables is important when implementing an early warning system for the proliferation of CB in drinking water sources [START_REF] Hu | Weather Variability, Sunspots, and the Blooms of Cyanobacteria[END_REF]. In the United States, NOAA issues weekly harmful algal bloom reports based upon satellite imagery and meteorological conditions [START_REF] Wynne | Evolution of a Cyanobacterial Bloom Forecast System in Western Lake Erie: Development and Initial Evaluation[END_REF]. Early warning systems are needed for drinking water treatment plant operators who must adjust treatment according to the biovolume and species of CB present, including such measures as changing coagulant, or adjusting the doses of treatment chemicals to be added [START_REF] Zamyadi | Speciesdependence of cyanobacteria removal efficiency by different drinking water treatment processes[END_REF]. For drinking water treatment, in sources highly impacted by cyanobacteria, operational decisions must be made in real time. Thus, weekly forecasts can provide useful information, but higher resolution forecasts are required (i.e. daily or sub-daily).

For the development of a reliable index which can serve as a 'universal threshold' for the occurrence of CB blooms in the aquatic environment, the effects of endogenous factors and exogenous parameters must be investigated. It is anticipated that the threshold for CB bloom occurrence will be based on the coupled effects of these external and internal factors. Thus, the objectives of this study are to: (1) through a new approach, propose an index in the form of a unique factor that integrates the effects of exogenous factors such as meteorological variables and water temperature on CB occurrence at a drinking water treatment plant intake, (2) evaluate this index to determine the probability of CB occurrence as an aid to operational decision making for water treatment plants, and (3) assess the role of nutrient limitation and the DIN/TP mass ratio as additional variables to be monitored along with phycocyanin at drinking water intakes. This paper proposes a novel index approach integrating meteorological factors and continuous phycocyanin probe measurements to determine drinking water intake risk as a result of intense cyanobacterial blooms.

Methodology 2.1 Site description

The research was performed at Missisquoi Bay, a large bay (77.5 km 2 at the latitude and longitude coordinates of of Lake Champlain, which straddles the Canada -United States border (Supplemental Information, Fig. SD-I1). The characteristics of tributaries are as described by [START_REF] Medalie | Use of flow-normalization to evaluate nutrient concentration and flux changes in Lake Champlain tributaries, 1990-2009[END_REF] and [START_REF] Adhikari | Temporal variability of phosphorus flux from Pike River watershed to the Missisquoi Bay of Quebec[END_REF]. Missisquoi Bay is a shallow water body with a mean depth of 2.8 m [START_REF] Galvez | Canada-USA Agreement on the Eutrophication Control and Protection of Lake Champlain Quebec[END_REF]. [START_REF] Myer | Limnology of Lake Champlain[END_REF] reported that in areas without rooted aquatic vegetation, the wave base can be deep enough during storms to cause mixing of bottom sediments. It is also the principal source of drinking water for a population of approximately 4100 residents (Statistics Canada, 2006). From 1992 to 2009, total phosphorus seasonal mean concentrations (sampling from June to September) increased in Missisquoi Bay by 72% (representing an increase of 20µg/L from an initial seasonal mean concentration of approximately 28µg/L). In contrast, total nitrogen seasonal mean concentrations decreased by approximately 25% (representing a decrease of approximately 0.15mg/L from an initial seasonal mean concentration of approximately 0.65mg/L) over the same period, possibly related to a reduction in atmospheric nitrogen deposition [START_REF] Smeltzer | Environmental change in Lake Champlain revealed by long-term monitoring[END_REF]. The reduction of the phosphorus is critical for the health of Missisquoi Bay and has led to a bilateral agreement between the province of Quebec (Canada) and state of Vermont (United States) to take action on reducing phosphorus to an acceptable limit of 25 µg/L (Comité inter-ministériel de concertation sur la Baie Missisquoi., 2003) and to guide related research activities [START_REF] Adhikari | Temporal variability of phosphorus flux from Pike River watershed to the Missisquoi Bay of Quebec[END_REF].

Field sampling and laboratory analyses

In 2007, a real-time monitoring program was developed and implemented in Missisquoi Bay using an online multi-probe system from YSI (YSI, Yellow Springs, Ohio) as described by [START_REF] Mcquaid | Use of in vivo phycocyanin fluorescence to monitor potential microcystin producing cyanobacterial biovolume in a drinking water source[END_REF]. The multi-probe system consists of eight sensors: temperature, depth, conductivity, turbidity, pH, dissolved oxygen, chlorophyll fluorescence (excitation: 470 nm; emission: 680 nm) and phycocyanin fluorescence (excitation: 590 nm; emission: 680 nm). It was installed to monitor the raw water intake of the drinking water treatment plant on the eastern shore of Missisquoi Bay every hour. Information on additional sampling and laboratory analyses of phycocyanin, chlorophyll, taxonomic counts and nutrients are provided in supplemental information (Section S1).

Meteorological data collection and treatment

Meteorological data were obtained from the weather station in Frelighsburg situated close to Missisquoi Bay. Frelighsburg data (Station name: FRELIGHSBURG; Province: QUEBEC; Latitude: 45.05°North; Longitude: -72.86°West; Altitude: 152.4m; climate identification: 7022579), provided by the National Climate Data and Information Archive (2013), were used to develop the model relating meteorological variables with cyanobacterial biovolume (CBV). Hourly data were available, thus the daily means of wind speed and direction were calculated using the method proposed in [START_REF] Roland | Meteorology for Scientists and Engineers[END_REF]. Statistica 9.0 (StatSoft.Inc, Tulsa, OK USA), Matlab (version 7.1.0.246 (R14) Service Pack 3, MathWorks, Inc.) and MS Excel (Microsoft, Redmond, WA, USA) were used for statistical analyses.

Development of a Meteorological Index

2.4.1 Description of the approach In order to evaluate the effects of meteorological factors on CB biomass, the development of this simple index is based on similar approach from [START_REF] Yu | Simulation of rice biomass accumulation by an extended logistic model including influence of meteorological factors[END_REF] for integrating the influence of climate components on rice biomass accumulation. [START_REF] Yu | Simulation of rice biomass accumulation by an extended logistic model including influence of meteorological factors[END_REF] proposed using a classical logistic growth model in which the maximum relative growth rate (µ) is defined as shown in equation (1):

) ( 0 p f × = µ µ (1)
where µ0 is the growth rate when meteorological factors are most favourable for population growth. The function f(p) reflects the meteorological factors being considered and can take on values between 0 and 1 with f(p) = 1 signifying growth without the constraint of meteorological factors and f(p) = 0 signalling the end of growth because of completely unfavourable conditions. However, different from the approach by [START_REF] Yu | Simulation of rice biomass accumulation by an extended logistic model including influence of meteorological factors[END_REF], in this paper, our focus is on the parameterization of f(p) and the relationship between this index and observed biomass accumulation in Missisquoi Bay to provide a functional index at the drinking water treatment plant and to assess the importance of meteorological factors in driving CB blooms. [START_REF] Yu | Simulation of rice biomass accumulation by an extended logistic model including influence of meteorological factors[END_REF] integrate f(p) directly into their logistic growth model and use multiple linear regression to fit model parameters for the prediction of biomass growth. Our objective in this study was not to predict the amount of cyanobacterial growth, which would require continuous data on nutrient limitation which are not available and are not typically monitored at high frequency at drinking water treatment plant, but rather to estimate the probability of CB blooms occurrence that can be explained by a series of easily measured meteorological parameters. The meteorological index does not distinguish between the direct and indirect influences that meteorological conditions play on biomass accumulation and it is not possible to directly distinguish accumulation from growth. However, the primary objective is to determine the probability of bloom occurrence given meteorological conditions influencing their presence at a drinking water treatment plant intake. In this context, our expectation is to develop an approach that can easily integrate meteorological factors into an index to adequately describe both CB growth and accumulation under various environmental conditions.

The meteorological index was developed in three stages: Stage 1) Meteorological index development: Data collected during the years 2007 and 2008 were used to define functions and calibrate all the parameters. The first two years of the 5 years of data were selected for model calibration because they represented a broad range of bloom conditions. In 2007, there was little accumulation of cyanobacteria at the drinking water intake in contrast to 2008 when intense blooms were observed. Additional data were not used in calibration in order to avoid -over-fitting parameters. The definition of this function f(p) takes into account the physical position of the drinking water intake relative to the wind direction.

Stage 2) Model validation: Data recorded from 2009 -2011 (YSI probe and FRELIGHSBURG weather station) were used to confirm the link between the function of factor f(p) on cyanobacteria biovolume (CBV). Similar model parameters found for 2007 and 2008 were used for the validation of the following years (2009 -2011).

Stage 3) Logistic regression model development: the logistic regression was developed with five years of data (2007)(2008)(2009)(2010)(2011) in order to estimate the probability of CB blooms occurrence given the meteorological conditions included in the f(p) function.

Selection of Meteorological Parameters for f(p)

The selection of meteorological variables was based upon local considerations (position of drinking water intake) and literature review. Wind direction measured in degrees is not monotonically related to CBV because only a specific range of wind directions can physically have an impact on the drinking water treatment plant based on its location. Even with strong winds, if the wind is moving in the direction away from the plant, CB accumulation at the intake is unlikely. There exists an optimal wind direction that could transport CB towards the plant. Thus, wind direction was included in the model based on physical rather than statistical considerations. The model used water temperature (as opposed to air temperature) because it more directly represents the environmental conditions experienced by the cyanobacteria. Similar variables (wind speed and direction, and water temperature) were used by [START_REF] Liu | Modelling the Effect of Weather Conditions on Cyanobacterial Bloom Outbreaks in Lake Dianchi: a Rough Decision-Adjusted Logistic Regression Model[END_REF] for the development of a model to estimate the probability of bloom occurrence. The measurements of CBV at the drinking water intake represent CB that are at the bottom of the water column, and thus, lower water temperatures can also be related to the presence of cyanobacteria in the lower depths of the water column because of additional processes not considered in the index model (mixing, cyanobacteria behavior with regards to light intensity, cloud cover, etc.). It has been reported that cyanobacteria move towards the bottom of the water column at night when water temperatures are lower [START_REF] Walsby | The role of gas vacuoles and carbohydrate content in the buoyancy and ertical distribution of Anabaena minutissima in Lake Rotongaio, New Zealand[END_REF]). Thus, when cyanobacteria are abundant, the dominant risk factors for a drinking water intake are related more to transport phenomena than to cyanobacterial growth. However, over the course of a season, the temperature becomes important when it falls outside of the optimal range for cyanobacterial growth. The model was constructed such that it would be possible to modify parameters according to local conditions. For example, various water bodies will have different dominant CB species. As such, optimal water temperatures could be adapted according to CB species present. However, this approach remains to be tested for systems other than Lake Champlain.

Definition of parameter functions for f(p)

Following the selection of the variables, f(p) was defined as follows:

( )

) g(p T g ) f(p i Wtemp i × = (2) direction) g(Wind Speed) g(Wind ) g(p i × = (3)
where Wtemp T represents the water temperature as measured by the probe installed at the drinking water intake, ( )

Wtemp T g
is the function describing the influence of water temperature on CB growth, i p is the meteorological parameter and ) g(p i is the function representing the influence of wind speed and wind direction on CB accumulation and growth. Considering meteorological parameters X1, X2, X3…XN for a given time series, the function ) g(p i at a given time can be written in the form of:

MAX i K(i)/K ) g(p = (4) with [ ] ∑ × × × × = ) F(X ... ) F(X ) F(X ) F(X K(i) N 3 2 1 (5) and max(K(i)) K MAX = (6)
From equation ( 4) we can see that ) g(p i will have a value between 0 and 1.

Hence, the expression for the growth and accumulation rate of CB integrating meteorological parameters and water temperature is:

( ) ) g(p T g µ µ i Wtemp 0 × × = (7)
and f(pi) is therefore equal to µ. The procedure for f(p) calculation is provided in Fig. SD-I2.

Function for water temperature

For conditions where water temperature exhibits a large range of variability, the description of the function ( ) Wtemp T g can adopt the following form [START_REF] Lehman | The assumptions and rationales of a computer model of phytoplankton population dynamics[END_REF]USEPA, 1985;[START_REF] Bouarab | Croissance en autotrophie et en mixotrophie de la microalgue Micractinium pusillum Fres. isolée d´un lagunage naturel : influence de la lumière et de la température[END_REF]: As Tmin < Topt <Tmax, these temperature levels will define four distinct intervals on which i) negligible growth occurs below Tmin or above Tmax; ii) increasing growth occurs from Tmin or Tmax to Topt ; and iii) decreasing growth occurs from Topt to Tmin or to Tmax (USEPA,1985).

( )                 - - × - =
By adjusting parameters, equations 8 and 9 can be applied for the environmental conditions at Missisquoi Bay or other source waters where there is a large variation of temperatures. Based on studies of [START_REF] Mooij | Predicting the effect of climate change on temperate shallow lakes with the ecosystem model PCLake[END_REF], [START_REF] Tang | Cyanobacterial dominance of polar freshwater ecosystems: are high-latitude mat-formers adapted to low temperature[END_REF] and Seaburg et al. (1982) minimum, optimal and maximum growth temperature are approximately 0°C, 25°C and 35°C, respectively. These studies have demonstrated that out of this temperature range, the growth is undetectable and that the majority of cyanobacterial species have an optimal growth temperature of 25°C. Additional information on the cyanobacterial species present in the Bay and their optimal growth temperatures are provided in supplemental information (Table SD-1).

Functions for wind speed and wind direction

The function describing the effect of wind direction (RDIR) is

( ) [ ] C a exp RDIR F × - = (10) ( ) ) Wind )(Wind Wind (Wind Wind Wind C dir_X dir_Max dir_Min dir_X 2 dir_opt dir_X - - - = (11)
where Winddir_X is the wind direction (°), Winddir_Min and Winddir_Max are its minimal and maximal values, Winddir_opt is its optimal value for CB accumulation, and a is a dimensionless calibration parameter.

The forms of equations 10 and 11 are similar those of equations 8 and 9. Similar to water temperature, wind direction is expected to have an optimal value in relation to the position of the drinking water treatment plant intake (the treatment plant is situated in the north-east region of the bay, at around 30 degrees). A similar procedure is used for describing the effect of wind speed. The function describing the effect of wind speed (RWND) is defined as

( ) [ ] B a exp RWND F × - = (12)
( )

) Wind )(Wind Wind (Wind Wind Wind B speed_X speed_Max speed_Min speed_X 2 speed_opt speed_X - - - = (13)
where Windspeed_X is the wind speed (m/s), Windspeed_Min and Windspeed_Max are its minimal and maximal values, Windspeed_opt is its optimal value for CB accumulation, and a is a dimensionless calibration parameter.

Weak or intense hydrodynamic fluctuations influence cyanobacteria [START_REF] Wu | Dynamics of cyanobacterial bloom formation during short-term hydrodynamic fluctuation in a large shallow, eutrophic, and wind-exposed Lake Taihu, China[END_REF]. Additional information on the influence of wind speed on cyanobacterial drift for the parameterization of the model is provided in supplemental information (Section S2). Based on the reported values in the literature, model calibration was performed to obtain the appropriate values for wind speed (Table 1). Equations 8, 11, and 13 rely on optimal ranges for the various variables. If measured values of a given factor fall outside of its optimal range, the index would predict no risk to the drinking water intake. For example, if wind speed were too low, minimal transport of CB towards an intake would occur. Thus, when variables fall outside of their optimal ranges, the risk that CB would be present at an intake is expected to be low and controlled by other factors not considered in the index such as light intensity, nutrients, and the vertical movement of CB in the water column.Model parameters were calibrated with data for 2007 and 2008 and were used for validation with data for 2009, 2010 and 2011. Logistic regression was performed on all data (2007 to 2011) in order to evaluate the probability that a CB bloom occurs at the drinking water intake in relation to the f(p) factor.

Results and Discussion

Meteorological parameters and cyanobacterial biovolume: model parameters

In this paper, our hypothesis is that exogenous factors dominate and drive the appearance of CB blooms at a drinking water intake where nutrients are not limiting.

A simple model was proposed to evaluate easily monitored exogenous variables on CB blooms occurrence in a water body via a meteorological index which assembles different weather variables into a unique index. This study proposes a simple, novel, but highly applicable index approach to evaluate drinking water intake exposure to cyanobacterial bloom occurrences due to their position in the water body.

As described in Section 2.4, the years 2007 and 2008 were used for calibration in order to determine the best values for the minimum, maximum and optimum values of water temperature, wind speed and direction and the factor 'a' to include in equations 8 to 12. A low value of f(p) was expected for most of 2007 when CB blooms were not observed at the drinking water intake (CBV < 1 mm 3 /L), in contrast to 2008 when many CB blooms were observed. Detailed data with regards to taxonomic counts for 2007 and 2008 are available in [START_REF] Mcquaid | Use of in vivo phycocyanin fluorescence to monitor potential microcystin producing cyanobacterial biovolume in a drinking water source[END_REF]. Approximately 98% of the CBV in 2007 and 2008 consisted of various species of Microcystis and Anabaena. The year 2007 was an exceptional year for Missisquoi Bay because no CB blooms had been observed until September [START_REF] Fortin | Detection of microcystinproducing cyanobacteria in Missisquoi Bay, Quebec, Canada, using quantitative PCR[END_REF] and no CB were present in large densities at the drinking water intake situated on the north-east side of the Bay. The probe was installed in the north-eastern region of the study site, thus even if a bloom occurred in Missisquoi Bay, the risk to the north-eastern region, where the probe was installed was low. The wind direction was such that CB would be transported to the outlet (south-western side of the bay) when CB blooms occurred and wind from the north and east were dominant. This flushing phenomenon observed in many lakes [START_REF] Dillon | The phosphorus budget of Cameron Lake, Ontario: The importance of flushing rate to the degree of eutrophy of lakes[END_REF] was dominant during the entire monitored period of 2007. In contrast, in summer 2008, the factor f(p) increased during (or one day prior) to high CBV values. The relationship of the meteorological factor f(p) with CBV at the drinking water intake is shown in Fig. 1a andb. It should be noted that CBV values below the threshold of 1mm 3 /L represent low concentrations of CB and thus can be considered as periods of low risk to the drinking water treatment plant despite greater uncertainty with regards to their exact concentrations [START_REF] Mcquaid | Use of in vivo phycocyanin fluorescence to monitor potential microcystin producing cyanobacterial biovolume in a drinking water source[END_REF]. The evolution of all the three functions calculated (for water temperature, wind speed and direction) was compared with CBV in order to observe the influence of each meteorological parameter (Fig. 1a andb). The relative importance of the factors was not considered in the analysis and all factors were given equal weighting. No single function could explain a CBV, but through f(p), all of these variables demonstrated an effect on the variation of CBV at the drinking water intake.

Few explanations have previously been available with regards to the non-existence of CB blooms observed during the summer 2007 at the drinking water intake at Missisquoi Bay. Based on this study, we can propose as an explanation that, for 2007 the value of f(p) which was very low except for a few days, clearly demonstrated that 2007 had meteorological conditions that were not conducive for CB growth or accumulation near the drinking water intake. This supports our theoretical model of the importance of wind (speed and direction) and water temperature as explained in supplemental information (Section S2 and S3).

Given that 2007 was not conducive for CB growth or accumulation, it is interesting to note that the flushing effect could also have prevented large blooms from forming, dying, settling to the bottom sediments and releasing additional nutrients that would have become available to create a new cycle of CB blooms. Redox conditions in the sediments [START_REF] Pant | Phosphorus sorbtion characteristics of estuarine sediments under different redox conditions[END_REF] for 2007 at Missisquoi Bay showed that anaerobic conditions did not occur as is typical for other years when large cyanobacterial blooms occur [START_REF] Smith | Relating sediment phosphorus mobility to seasonal and diel redox fluctuations at the sediment-water interface in a eutrophic freshwater lake[END_REF]. Through the use of a combination of different variables in one index (f(p)), it is possible to improve explanations of environmental conditions that contribute to CB bloom occurrences. This has been highlighted by other studies using statistical methods to explain CB bloom occurrence [START_REF] Hu | Weather Variability, Sunspots, and the Blooms of Cyanobacteria[END_REF][START_REF] Liu | Modelling the Effect of Weather Conditions on Cyanobacterial Bloom Outbreaks in Lake Dianchi: a Rough Decision-Adjusted Logistic Regression Model[END_REF].

The investigating results demonstrated that the effects of wind and water temperature cannot be ignored with regards to CB bloom occurrence at any considered position in the water body as the spatial distribution of CB is affected by wind direction and speed and water temperature [START_REF] Kanoshina | The influence of weather conditions (temperature and wind) on cyanobacterial bloom development in the Gulf of Finland (Baltic Sea)[END_REF]. The importance of wind has been emphasized by [START_REF] Izydorczyk | Measurement of phycocyanin fluorescence as an online early warning system for cyanobacteria in reservoir intake water[END_REF] although they did not quantitatively demonstrate the wind effect on cyanobacteria. The wind speed can induce mixing in the water column and influence nutrient availability (e.g. [START_REF] Macintyre | Boundary Mixing and Nutrient Fluxes in Mono Lake, California[END_REF][START_REF] Livingstone | Impact of secular climate change on the thermal structure of a large temperate central European lake[END_REF]O'Reilly et al., 2003;[START_REF] Schmittner | Decline of the marine ecosystem caused by a reduction in the Atlantic overturning circulation[END_REF].

Meteorological parameters and cyanobacterial biovolume: Model application

Model application for 2009, 2010 and 2011 show a similar trend between f(p) and CBV. As seen in Fig. 2a, b and c, the f(p) index is low but CBV is high in some cases (false negatives). However, the false negatives were frequently related to bloom events that occurred over multiple days when the meteorological index decreased faster than the CBV (returning to a state without a bloom). In other instances, false negatives were likely related to endogenous and exogenous factors optimal for cyanobacterial growth but not considered in the index (e.g. solar radiation, nutrients). The benefit of the current index approach is its potential applicability to cyanobacteria impacted drinking water intakes with data that can be easily and reliably obtained in real time. Thus, drinking water treatment plants can evaluate their vulnerability (the probability of bloom occurrence, a maximum probability occurs with a maximum f(p)) using meteorological conditions and water temperature without the need to wait for satellite imagery. CB occurrence is a complex phenomenon due to the non-linearity of the variables, as described previously. Thus, a logistic regression model was used to estimate the probability of CB bloom occurrence given by the calculated index (f(p)). In order to develop a logistic regression model, a CBV of 1 mm 3 /L was considered as a threshold value for when the drinking water treatment intake is at risk. The 1 mm 3 /L is a locally adapted threshold value and was proposed as the monitoring threshold (alert level) for Missisquoi Bay by [START_REF] Mcquaid | Use of in vivo phycocyanin fluorescence to monitor potential microcystin producing cyanobacterial biovolume in a drinking water source[END_REF][START_REF] Mcquaid | Use of in vivo phycocyanin fluorescence to monitor potential microcystin producing cyanobacterial biovolume in a drinking water source[END_REF] found that 1mm 3 /L corresponded to a maximum potential microcystin concentration of 13µg/L in Missisquoi Bay. The [START_REF] Mcquaid | Use of in vivo phycocyanin fluorescence to monitor potential microcystin producing cyanobacterial biovolume in a drinking water source[END_REF] threshold falls between alert level 1 (biovolume: 0.2mm 3 /L; maximum potential microcystin concentration (MPMC): 2.6µg/L) and alert level 2 (biovolume: 10 mm 3 /L; MPMC: 130µg/L) proposed by Bartram et al. (1999) and also represents the lowest biovolume value for which phycocyanin probes provided reliable estimates in Missisquoi Bay. For CB blooms to develop, it is likely that conditions must be favorable for more than one consecutive day. The logistic model developed was based on f(p) and CBV data from the five years (2007 to 2011). All the parameters values of the model are significant with p-value equal to 0.006 with an odds ratio of 6.03 (Table 2). As expected, meteorological conditions were observed to have a large influence on CBV increases at the drinking water intake. The odds ratio found for a logistic regression model demonstrated a significant relationship between f(p) and the probability of CB occurrence above the CBV threshold of 1 mm 3 /L. Thus, increasing values for f(p) coincide with CB blooms at the drinking water intake. As seen in Fig. 3, a maximum probability of 68% of CB bloom occurrence can be explained by the increase of f(p). An f(p) value of 1 means that the combined effects of wind speed and direction, and water temperature are optimal for CB accumulation at the drinking water intake. It is also interesting to note that at low f(p) values, there is a low probability (less than 26%) of a CB bloom at the intake. Thus, the probability of a false negative is less than 26%. As described in Section 3.2, false negatives were related to the index decreasing faster than the CBV, and both false positives and negatives also arose potentially from the exclusion of factors such as light intensity and nutrients from the index. Our index approach can be discussed within the context of other statistical approaches that have included regression models [START_REF] Smith | Predictive models for the biomass of blue-green algae in lakes[END_REF][START_REF] Smith | Predicting the Summer peak biomass of four species of blue-green algae (cyanophyta:cyanobacteria) in Swedish lakes[END_REF], correlations [START_REF] Varis | Associations between lake phytoplankton community and growth factors-a canonical correlation analysis[END_REF], classification and regression trees [START_REF] Hu | Weather Variability, Sunspots, and the Blooms of Cyanobacteria[END_REF], artificial neural networks [START_REF] Teles | Time series forecasting of cyanobacteria blooms in the Crestuma reservoir (Douro River, Portugal) using artificial neural networks[END_REF][START_REF] Wei | Use of artificial neural network in the prediction of algal blooms[END_REF][START_REF] Maier | Use of artificial neural networks for modelling cyanobacteria Anabaena spp. in the River Murray, South Australia[END_REF][START_REF] Recknagel | Artificial neural network approach for modelling and prediction of algal blooms[END_REF], and conditional probabilities [START_REF] Peretyatko | Assessment of the risk of cyanobacterial bloom occurrence in urban ponds: probabilistic approach[END_REF]. For many statistical approaches it is not possible to distinguish factors that are drivers of CB blooms from those that are merely correlated with CB blooms. For example, although pH is typically correlated with CBV and has often been included in regression models, it cannot be predictive because changes in pH are driven by changes in CBV (Zamyadi et al. 2012b). By studying the probability of CBV conditioned on environmental factors (not including history of CBV), [START_REF] Peretyatko | Assessment of the risk of cyanobacterial bloom occurrence in urban ponds: probabilistic approach[END_REF] reached a maximum probability of 0.5. In addition, some of the environmental factors, such as pH, were not drivers of CB blooms. In contrast, our index approach found a maximum probability of 0.68 with the maximum index value using only variables that can be viewed as potential drivers of CB accumulation at a drinking water intake. Thus, the index approach provides a robust estimate of the potential for CB accumulation at a drinking water intake without relying on variables that cannot be monitored continuously, or that are correlated with CB because the CB themselves were the drivers of the change. Although artificial neural networks and classification and regression trees are powerful techniques for highly non-linear problems, they are black box approaches that do not necessarily encode the available physical knowledge of the system. Using the index approach, we were able to include information and factors based on the physical environment (e.g., the specific location of the drinking water intake) to obtain reliable results for operational decision-making at drinking water treatment plants.

For our model based on the assumption that nutrients are not limiting for CB, a better estimate is not expected given that only limited meteorological variables were considered and major endogenous factors related to growth like nutrients were not included in the model. The unexplained part of CB occurrence is related to endogenous factors such as nutrient availability and limitation. Thus, we explored the use of a nutrient mass ratio that could be used to augment continuous phycocyanin monitoring at the drinking water intake. CB abundance was compared to the DIN/TP mass ratio. Total, rather than soluble reactive phosphorus was used in the ratio because the algae and cyanobacteria population can store sufficient phosphate to increase 10-fold, even if no dissolved phosphate can be detected [START_REF] Chorus | Cyanobacteria and Algae: Monitoring Bathing Waters -A Practical Guide to the Design and Implementation of Assessments and Monitoring Programmes: WHO[END_REF]. The mass ratio DIN/TP provides additional information on nutrient limitation that will impact cyanobacterial growth [START_REF] Bergstrom | The use of TN:TP and DIN:TP ratios as indicators for phytoplankton nutrient limitation in oligotrophic lakes affected by N deposition[END_REF][START_REF] Ptacnik | Performance of the Redfield Ratio and a Family of Nutrient Limitation Indicators as Thresholds for Phytoplankton N vs. P Limitation[END_REF]. We found a significant negative correlation between mass ratio DIN/TP and extracted CB phycocyanin pigment (r 2 =-0.48, p-value < 0.05). This correlation suggests that CB abundance is likely to be limited by either the DIN or TP. Two hypotheses are proposed 1) CB dominate under nitrogen limitation due to their efficiency to use trace levels of DIN, and 2) CB are more competitive for phosphorus than other algal species as CB are efficient at nutrient uptake and storage [START_REF] Sorokin | Dynamics of phosphorus in the Venice lagoon during a picocyanobacteria bloom[END_REF]. Low heterocytes in Missisquoi Bay during summer bloom periods suggests that nitrogen fixing cyanobacteria do not rely on atmospheric N2 [START_REF] Mccarthy | Implications of water column ammonium uptake and regeneration for the nitrogen budget in temperate, eutrophic Missisquoi Bay, Lake Champlain (Canada/USA)[END_REF]. Fig. 4 shows the nutrient mass ratio, DIN/TP, from weekly duplicate samples in 2011 from all sites described in Supplemental Information (Figure S1) and its correlation with CB phycocyanin pigment. It shows that when the mass ratio DIN/TP is low, the phycocyanin concentration, related to CBV, can reach high values. Furthermore, it shows that nutrient limitation does occur and could explain the false positives in the index model. To improve index predictions further, it would be useful to consider the addition of continuous nutrient and light intensity data that were not available in this study. In Fig. 4, the critical zone where CB abundance in Missisquoi Bay can be high, is for a DIN/TP mass ratio below 3 and it descends exponentially, similar to values reported between the range of 1-2 by [START_REF] Bergstrom | The use of TN:TP and DIN:TP ratios as indicators for phytoplankton nutrient limitation in oligotrophic lakes affected by N deposition[END_REF] and [START_REF] Ptacnik | Performance of the Redfield Ratio and a Family of Nutrient Limitation Indicators as Thresholds for Phytoplankton N vs. P Limitation[END_REF]. Thus, DIN/TP is potentially a useful indicator of nutrient limitation, although nutrient availability may be driven by phytoplankton dynamics [START_REF] Mccarthy | Implications of water column ammonium uptake and regeneration for the nitrogen budget in temperate, eutrophic Missisquoi Bay, Lake Champlain (Canada/USA)[END_REF]. The critical value of 3 obtained in our study could be used as an additional monitoring index for CB abundance, albeit not in real time and not necessarily for prediction. However, it could be useful when the meteorological index indicates a high vulnerability to CB occurrence at the intake.

The meteorological index proposed in this study is most important for determining initial conditions that could lead to a bloom rather than for explaining when the bloom will disappear. This is related to practical considerations at the drinking water plant intake with regards to operational decision-making. For example, drinking water treatment plant operators could take samples to determine optimal coagulants and doses [START_REF] Zamyadi | Speciesdependence of cyanobacteria removal efficiency by different drinking water treatment processes[END_REF], or they could opt to close their drinking water intake and use stored water should CBV increase rapidly. At the study site, current practice involves the use of phycocyanin monitoring exclusively for assessing onsite CB vulnerability. With the use of real-time phycocyanin monitoring and the meteorological index, the operators would have important information to determine whether or not optimal bloom conditions are occurring or are likely to occur that will enable them to plan and to closely monitor the CBV at the drinking water intake.

Although the demonstration of the approach was conducted at Missisquoi Bay, this index approach for surface waters other than Lake Champlain could be applied and validated using locally adapted threshold values and variables. The different variables used to calculate this index can easily be collected for any water body affected by CB and used to evaluate the risk of CB blooms affecting a zone of interest. Findings from the meteorological index model in our research strongly reinforce our further investigation on a coupled hydrodynamic-cyanobacterial growth model to fully consider meteorological variables and their influence on lake hydrodynamics and cyanobacterial transport.

Conclusions

An early warning system for drinking water treatment plants should consider meteorological factors in addition to in vivo phycocyanin monitoring in order to adequately respond to CB events. Wind (speed and direction) and water temperature are related to CB blooms, both directly with regards to accumulation near a drinking water intake and indirectly on factors related to CB growth. This study shows that all the variables have an effect on CBV at drinking water intake and must be analysed together to understand their effect on CBV at drinking water intakes. As nutrient limitation is related to CB occurrence, the DIN/TP mass ratio could be monitored in conjunction with phycocyanin and meteorological variables.

Favourable meteorological conditions (high f(p) value) was related to maximum probability of 68% that a bloom event with CBV exceeding 1 mm 3 /L would occur. As CB occurrence at drinking water treatment intakes are largely driven by exogenous factors, the use of hydrodynamic models could be used to more accurately predict CB events at drinking water intakes and to propose solutions to avoid breakthrough of cyanobacteria into drinking water treatment plants.

  a is a dimensionless calibration parameter. This approach takes into account three different levels of water temperature: (1) minimum growth temperature Tmin ( o C), (2) optimum growth temperature Topt ( o C) and (3) maximum growth temperature Tmax ( o C).

Fig. 1 :

 1 Fig.1: Functions defined from weather parameter versus CBV (model development): a) for 2007 and b) for 2008

Fig. 2 :

 2 Fig.2: Meteorological parameter f(p) versus CBV (model application): a) for 2009, for b) for 2010 and c) for 2011

Fig. 3 :

 3 Fig.3: CB blooms occurrence probability versus meteorological factor f(p)

Fig. 4 :

 4 Fig.4: Changes in the mass ratio of dissolved inorganic nitrogen (DIN) to total phosphorus (TP) concentration, and cyanobacteria phycocyanin (PC) concentration Missisquoi Bay in 2011 (mesured in lab)

Table 1 :

 1 Calibrated Parameters values for 2007 and 2008

	Model		Calibrated Parameters
	parameters	a	Minimum	Optimum	Maximum
	Wind direction ( o )	7	190	200	310
	Wind speed (m/s)	0,1	2	3	6
	Temperature ( o C)	7	10	20	35

Table 2 :

 2 logistic regression between cyanobacterial occurrence bloom and f(p) during the

	five years (2007 to 2011) (	t	=	log	  	1	-P	P	  	; P=probability of occurrence)
	Year	Model Equation			p-value		Odds ratio	Wald's Chi-square	Weather datas sources
	2007 to 2011	t = -0.018+1.845×f(p)	0.006 (n=393)	6.03	11.76	Freligsburg weather station

Acknowledgments

This study has been financially support by Engineering Research Council of Canada (NSERC), Canada Foundation for Innovation (CFI) , Fonds de Recherche du Québec -Nature et Technologies (FRQNT), Canada Research Chairs (CRC), Canadian International Development Agency (CIDA). The authors acknowledge the support of Yves Fontaine, Julie Philibert, Jacinthe Mailly and the operators of drinking water treatment plant at Missisquoi Bay for their technical help.