Yasir Hamid

Ludovic Journaux

John Aldo Lee

Lucile Sautot

Bushra Nabi

M Sugumaran

Large-scale nonlinear dimensionality reduction for network intrusion detection

come

Large-scale nonlinear dimensionality reduction for network intrusion detection

Yasir Hamid 1 , Ludovic Journaux 2 , John A. Lee 3 , Lucile Sautot 4 , Bushra Nabi 5 , M. Sugumaran 1

1-Department of C.S.E., Pondicherry Engineering College, Pondicherry-14, India 2-LE2I UMR6306, CNRS, Univ. Bourgogne Franche-Comté, AgroSup Dijon, France 3-Univ. catholique de Louvain, UCL/SSS/IREC/MIRO, Bruxelles, Belgium 4-AgroParistech, UMR TETIS , Maison de la télédétection, Montpellier, France 5-Department of Botany, University of Kashmir Abstract. Network intrusion detection (NID) is a complex classification problem. In this paper, we combine classification with recent and scalable nonlinear dimensionality reduction (NLDR) methods. Classification and DR are not necessarily adversarial, provided adequate cluster magnification occurring in NLDR methods like t-SNE: DR mitigates the curse of dimensionality, while cluster magnification can maintain class separability. We demonstrate experimentally the effectiveness of the approach by analyzing and comparing results on the big KDD99 dataset, using both NLDR quality assessment and classification rate for SVMs and random forests. Since data involves features of mixed types (numerical and categorical), the use of Gower's similarity coefficient as metric further improves the results over the classical similarity metric.

Introduction

With internet now being ubiquitous in many aspects of daily life, keeping the networks secure has become critical and requires continuous effort [START_REF] Stallings | Network security essentials: applications and standards[END_REF]. In this context, an Intrusion Detection System (IDS) aims to discriminate the legitimate network connections against potentially harmful ones [START_REF] Rowland | Intrusion detection system[END_REF].

Effective and popular processing techniques for IDS are feature selection and extraction [START_REF] Liao | Intrusion detection system: A comprehensive review[END_REF]. Regarding the former, works about IDS have explored a lot of feature selection techniques [START_REF] Xiang | Design of multiple-level hybrid classifier for intrusion detection system using bayesian clustering and decision trees[END_REF]. For the latter, though, the well known PCA remains largely dominant [START_REF] Heba | Principle components analysis and support vector machine based intrusion detection system[END_REF]. PCA reduces data dimensionality based on variance preservation but this simple criterion is not guaranteed to retain relevant information in data. In an IDS, PCA generally improves the detection rate of the dominant classes, i.e., Probe and Denial of Service (DOS), but it does not impact much on the detection rate of more marginal classes like User to Root (U2R) and Remote to Local (R2L) [START_REF] Lin | CANN: An intrusion detection system based on combining cluster centers and nearest neighbors[END_REF]. PCA is one of the oldest linear methods of dimensionality reduction (DR) and many nonlinear methods have been developed over the past decades [START_REF] Lee | Nonlinear dimensionality reduction[END_REF]. A usual assumption of DR methods is that the intrinsic dimensionality of data is much lower than the ambient space dimensionality, e.g., because data lives on a submanifold. Recent nonlinear DR methods have proved to be effective in many applications, like visualization and exploratory analysis [START_REF] Van Der Maaten | Visualizing data using t-sne[END_REF]. Scalability issues of these methods have been addressed lately [START_REF] Yang | Optimization equivalence of divergences improves neighbor embedding[END_REF].

In classification tasks like IDS, preprocessing with DR remains debatable, with possible gains (mitigation of the curse of dimensionality, noise attenuation) and losses (blindness with respect to the actual supervised task, degraded separability due to collapsing classes). This paper shows that cluster magnification in some DR methods [START_REF] Van Der Maaten | Visualizing data using t-sne[END_REF] lead to a positive outcome. Another challenge is to exploit both qualitative and quantitative features in the dataset, which is tackled by replacing the Euclidean distance with Gower's dissimilarity coefficient [START_REF] Gower | A general coefficient of similarity and some of its properties[END_REF].

The rest of the paper is organized as follows. Section 2 presents our material and methods (DR, quality assessment of DR, and Gower's coefficient). Section 3 presents the experimental results. Finally, Section 4 draws the conclusions.

Materials and methods

Real-world KDD99 dataset

The KDD99 dataset is the benchmark that is predominantly used for network intrusion detection1 . It comes with mixed-type features, some being numeric num_failed_logins, num_file_creations, num_access_files and others being nominal like Protocol, type_of_service. A total of M =4 1features are present for each record, while the 42 nd is the class label. The dataset includes a total of 23 attacks spanning across four groups, User to Root, Remote to Local, DoS, and Probe. In addition to four attack groups, a class is associated with normal network connections. So, let Ξ =[ξ i ] 1≤i≤N be the M -by-N dataset. For our experiments with KDD99, a subset of N = 61906 records was selected, a size that raises already significant scalability issues. Moreover this subset size also allow acceptable computing time.

Gower's similarity coefficient

Gower's similarity coefficient is often used in an ecological study or in a modeling work. It is designed to measure (dis)similarity between two items or individuals that are characterized by both numerical and categorical features [START_REF] Gower | A general coefficient of similarity and some of its properties[END_REF]. Gower's coefficient is defined as a weighted average of similarity scores, assigned in the pairwise comparison of items. The scores are computed differently depending on the type of the features: categorical (dichotomous, nominal, interval, or ratioscale, ...) or numeric (real or integer quantities). Let us consider two items ξ i and ξ j with M features. Gower's coefficient is calculated as

M k=1 w k S k (ξ i , ξ j ),
where w k is a binary weight of the kth feature. The weight is one except if the kth feature value is not defined for ξ i or ξ j ; S k (ξ i , ξ j ) then depends on type V k of the kth feature:

S k (ξ i , ξ j )= δ(ξ ik ,ξ jk ) if V k is qualitative |ξ ik -ξ jk | max V k -min V k if V k is quantitative , ( 1
)
where δ is Dirac's function, which is one iff ξ i and ξ j are from the same category.

Dimensionality reduction

Dimensionality reduction is widely used for exploratory data analysis, as it can help users to visualize data in low-dimensional spaces. Many DR methods have been proposed in the literature [START_REF] Lee | Nonlinear dimensionality reduction[END_REF], like principal component analysis (PCA) [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF], multidimensional scaling[1], Sammon's nonlinear mapping [START_REF] Sammon | A nonlinear mapping for data structure analysis[END_REF], stochastic neighbour embedding [START_REF] Hinton | Advances in neural information processing systems, chapter Stochastic neighbor embedding[END_REF], t-distributed stochastic neighbour embedding [START_REF] Van Der Maaten | Visualizing data using t-sne[END_REF], neighbour retrieval visualizer [START_REF] Venna | Information retrieval perspective to nonlinear dimensionality reduction for data visualization[END_REF], Jensen-Shannon embedding(JSE) [START_REF] Lee | Type 1 and 2 mixtures of kullback-leibler divergences as cost functions in dimensionality reduction based on similarity preservation[END_REF], and recently Multi-Scale Jensen-Shannon embedding (Ms.JSE) [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF]. Most of these approaches, though, scale rather poorly in a 'big data' context, namely, if N gets larger than a few thousand items, like in the KDD99 dataset. Except for PCA, these DR methods entail the computation of pairwise (dis)similarities in large N -by-N matrices, with quadratic time complexity. To overcome this limitations, pairwise relationships in SNE-like methods are aggregated and approximated with space partitioning trees [START_REF] Yang | Optimization equivalence of divergences improves neighbor embedding[END_REF]. In our experiments, we use in particular fast approximations of Elastic Embedding, SNE, symmetric SNE, NeRV, t-SNE, and weighted t-SNE2 . Afterwards, we focus on the best-performing method and replace the Euclidean distance with Gower's coefficient in order to compare results.

Quality criterion to compare DR methods objectively

To compare the DR results, we use the quality assessment tools (QA) that quantify neighborhood preservation [START_REF] Lee | Quality assessment of dimensionality reduction: rank-based criteria[END_REF][START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF]. If ν K i if the K-ary neighborhood of ξ i and n K i is the K-ary neighborhood of its low-dimensional counterpart x i found by DR, then

Q NX (K)= N i=1 ν K i ∩n K i KN is the average neighborhood preservation. Rescaling into R NX (K)= (N -1)Q NX (K)-K N -1-K
locates the considered embedding between a random one (Q NX (K) ≈ K, R NX (K)=0 ) and a perfect one (Q NX (K)=R NX (K)=1 ). A synthetic scalar score over all sizes K is the AUC of R NX (K) in plot with log abscissae, namely,

AUC(R NX (K)) = N -2 K=1 R NX (K)/K N -2
K=1 1/K . Like in DR methods, space-partitioning trees and incomplete sort (quick select) can improve scalability of DR QA.

Results and discussion

All DR methods listed in the previous section are applied to KDD99 to get 3D embeddings. For easier visualization, 2D embeddings corresponding to the best methods are also shown and discussed before reporting the classification results.

Quality assessment of DR. Figure 1 shows the quality curves (R NX (K)) and their AUC (in the legend box) for the 3D embeddings of all considered DR methods. Looking at all curves except the one for Gower t-SNE, we conclude that t-SNE outperforms the other methods, with the tallest and widest bump between K =10and 200, as well as highest AUC. Among those methods, t-SNE is the only method using discrepant Gauss/Student neighbourhoods in the highand low-dimensional spaces, respectively. Thereby, it induces an exponential stretch of distances and cluster magnification [START_REF] Lee | On the role and impact of the metaparameters in t-distributed stochastic neighbor embedding[END_REF]. The curve of Gower t-SNE is not straightforward to interpret, since QA tools do not account for the change of metric occurring with Gower's coefficient. As discussed below, its use is expected to improve classification results, since it exploits qualitative features that can increase class and cluster separability, as illustrated below. Best emdeddings of KDD99 with t-SNE. For visualization purposes, and since t-SNE outperforms other methods, Figure 2 shows 2D embeddings after running t-SNE with the Euclidean distance and Gower's coefficient. Some clusters seem to be monodimensional dashes, while disc-shaped ones would reflect a larger intrinsic dimensioanlity. Visual inspection shows that Gower t-SNE separates clusters even better than regular t-SNE. This is indirectly confirmed in the 3D embeddings by the improved classification results below.

Classification results. Table 1 reports the classification performance metrics of SVMs and RFs, applied to all 3D embeddings of KDD99. Due to the clear superiority of t-SNE in terms of cluster separation and longer computation times, we used Gower's similarity coefficient only with this DR method. SVMs with Gaussian kernel was tried with five different bandwidths Γ and 10-fold crossvalidation. The best model was kept. A random forest was constructed as an ensemble of random trees without any restriction on the depth of the tree and a batch size set to 100. The total number of iterations was explicitly set to 100 and the results were checked for 2 decimal points. To assess classification quality, we relied on: %Classification accuracy, Precision, Recal l, and ROC AUC.

Conclusion

In the context of network intrusion detection, this paper proposes a methodology in two steps. First, data are embedded in a lower-dimensional space. In the case of KDD99, a scalable method of dimensionality reduction like Barnes-Hut t-SNE mitigates the curse of dimensionality and magnifies cluster separation (as well as class separation, indirectly), provided the target dimensionality is not too low, in order to prevent classes to collapse on each other. To extract maximal information from the KDD99 mixed-type features, we also combined t-SNE with Gower's similarity coefficient. Experimental results demonstrate the effectiveness of the proposed methodology. Challenges for the future are: getting even better scalability, in order to process the entire KDD99 set; making the methodology fully parametric (the DR method and Gower's coefficient have no straightforward out-of-sample ex-tension); and applying the proposed methodology to other datasets in order to generalize and develop a new approach of mixed data analysis.

Figure 1 :

 1 Figure 1: Quality curves (R NX (K)) and their AUC (legend) for all methods.

  (a) t-SNE with the Euclidean distance (b) t-SNE with Gower's similarity coefficient

Figure 2 :

 2 Figure 2: Scatter plots of the t-SNE embeddings, without and with Gower.

Table 1 :

 1 Classification performance metrics

		% Classification	Precision	Recall	ROC Area	M.A.Error	RMSE
	RD Methods	SVM	RF	SVM RF SVM RF SVM	RF	SVM	RF	SVM	RF
	Normal Data	96.97	97.81	0.97	0.97	0.97	0.97	0.98	0.98 0.0012 0.0019 0.0319	0.0089
	EE	98.92	98.55	0.98	0.98	0.98	0.98	0.97	0.99 0.0009 0.0055 0.0306	0.0379
	SNE	99.17	98.81	0.99	0.98	0.99	0.98	0.99	0.99 0.0013 0.0058 0.3621	0.0366
	NerV	97.81	97.95	0.97	0.97	0.97	0.97	0.98	0.998 0.0019 0.0089 0.0436	0.0496
	t-SNE	99.60	99.54	0.99	0.99	0.99	0.98	0.99	0.99 0.0004 0.0005 0.0198	0.0169
	WTSNE	98.87	99.54	0.98	0.99	0.98	0.99	0.99	0.99 0.0015 0.0009 0.3116	0.0184
	t-SNE with Gower 99.96	99.97	1	1	1	1	0.99	1	0	0.0002 0.0051 0.0056

http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html

http://research.cs.aalto.fi/pml/software/ne/


Jie Chen
email: jchen@ee.ucr.edu

Peilin Fu

Silviu-Iulian Niculescu
email: silviu.niculescu@lss.supelec.fr

will become unstable if a critical zero enters the open righthalf plane, and otherwise will remain stable if all the critical zeros remain in the left-half plane. The essential problem in this approach, consequently, is to characterize the analytical properties of the imaginary zeros.

When Will Zeros of Time-Delay Systems Cross Imaginary Axis?

Jie Chen, Peilin Fu and Silviu-Iulian Niculescu Abstract-A time-delay system may or may not be stable for different periods of delay. When will then a delay system be stable or unstable, and for what ranges of delay? This paper attempts to answer these questions. We show that by finding a set of critical delay values, for which the system's characteristic quasipolynomial has zeros on the imaginary axis, it is possible to determine its stability in the full range of the delay parameter by characterizing the analytical behaviors of the zeros. This characterization is facilitated by an operator perturbation approach, which is both conceptually attractive and computationally efficient. The entire procedure, which first identifies the critical zeros on the imaginary axis and next determines whether the zeros cross the imaginary axis, requires only solving a generalized eigenvalue problem.

Index Terms-Time-delay, asymptotic stability, critical zeros, asymptotic behavior, matrix pencil.

I. INTRODUCTION

In this paper we study stability properties of linear timedelay systems. We consider specifically retarded systems with commensurate delays. For systems of this type, it is known that the stability problem can be tackled efficiently by computing the critical delay values that result in imaginary zeros of the system's characteristic quasipolynomial. This computation can be executed in a number of different ways, notably, by, e.g., solving a generalized eigenvalue problem [START_REF] Gu | Stability of Time-Delay Systems[END_REF], [START_REF] Niculescu | Stability and hyperbolicity of linear systems with delayed state: A matrix pencil approach[END_REF], [START_REF] Niculescu | Delay effects on stability. A robust control approach[END_REF]. Since the critical delay values form consecutive intervals and the zeros vary continuously with respect to the delay parameter, whether the system is stable over a specific delay interval can be determined by checking the stability corresponding to any fixed delay value within that interval, thus in principle, determining the system's stability over the entire range of delay values. One should note, nevertheless, that even for a fixed delay, the testing of stability for a time-delay system is not a simple task.

A more potent approach in accomplishing this goal is to analyze the asymptotic behavior of the critical zeros on the imaginary axis. The idea is rather simple. Specifically, at each critical delay value corresponding to a critical zero of the characteristic quasipolynomial, one may seek to determine whether the zero may cross the imaginary axis, from one half plane into another; for example, the system The problem alluded to above, known in the literature as the stability switch problem, has been studied for quite some time (see, e.g., [START_REF] Boese | Stability with respect to the delay: On the paper of K[END_REF], [START_REF] Cooke | Discrete delay, distributed delay and stability switches[END_REF], [START_REF] Cooke | On zeroes of some transcendental equations[END_REF], [START_REF] Hertz | Simplified analytical stability test for systems with commensurate time delays[END_REF], [START_REF] Lee | On the τ -decomposition method of stability analysis for retarded dynamical systems[END_REF], [START_REF] Olgac | An exact method for the stability analysis of time-delayed LTI systems[END_REF], [START_REF] Walton | Direct method for TDS stability analysis[END_REF]). Recently, the authors developed an operator perturbation approach [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF] to facilitate this analysis, which appears both analytically attractive and computationally appealing. The approach seeks to reformulate the zero asymptotic analysis problem as one of eigenvalue perturbation [START_REF] Kato | Perturbation Theory for Linear Operations[END_REF], and it amounts to performing a first-order asymptotic expansion of the zeros around their critical values. Whether a critical zero will cross the imaginary axis can then be checked by computing the eigenvalues of a constant matrix. It is worth noting that the asymptotic expansion can be carried out in a rather efficient manner; in fact, it can be considered a direct consequence resulted from the computation of the critical delay values and the imaginary zeros. This paper continues the work of [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF] and gives a full development of the asymptotic analysis up to the second order. The results, along with the first-order analysis given in [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF], thus provide a more complete tool which can aid in stability analysis when the first-order tests cease to be useful; this scenario arises when the first-order coefficient in the asymptotic series is imaginary, which renders the firstorder test inconclusive. Similar to their first-order counterparts, the second-order tests developed herein also require computing the eigenvalues of a constant matrix only.

The remainder of this paper is organized as follows. In Section 2, we first develop general eigenvalue perturbation results, which serve as our primary technical machinery. We then apply these results to time-delay systems; we first consider in Section 3 systems in the state-space form, and next in Section 4 systems described by differential-difference equations. Conditions are derived for characterizing whether a critical zero will cross the imaginary axis. Section 5 presents an illustrative example and Section 6 concludes the paper.

We use the following notation throughout the paper. Let IR be the set of real numbers, C the set of complex numbers, and IR + the set of nonnegative real numbers. Denote the open right half plane by C + := {s : (s) > 0}, the closed right half plane by C + , and the imaginary axis by ∂C + . Similarly, denote the open unit disc by ID, the unit circle by ∂ID, and the closed exterior of the unit disc by ID c . For a matrix A, denote its spectrum by λ(A), and the ith eigenvalue by λ i (A). For a matrix pair (A, B), denote the set of all its generalized eigenvalues by σ(A, B), i.e., σ(A, B) := {λ ∈ C : det(A -λB) = 0} .

The operation A⊕B denotes the Kronecker sum, and A⊗B the Kronecker product, of the matrices A and B.

II. PRELIMINARY RESULTS

This section introduces the operator perturbation theory for matrix eigenvalue problems. We shall first give an account of the first-order results developed in [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF], and subsequently perform a second-order analysis. The development is based on the classical work by Kato [START_REF] Kato | Perturbation Theory for Linear Operations[END_REF] and the authors' own extension in [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF].

A. First-Order Analysis

Consider a matrix operator T (x) of a real variable x. Suppose that in the neighborhood of x = 0, the perturbed operator T (x) is holomorphic, or equivalently, can be expanded into the power series,

T (x) = T (0) + xT (0) + x 2 T (0) + • • • . ( 1
)
It is known that near the origin, any semisimple eigenvalue of T (x), namely the semisimple root of the characteristic equation

det (T (x) -ξI) = 0, (2)
is an analytic functions of x and can be expressed as a power series in x; here by a semisimple eigenvalue, we mean a repeated but diagonalizable eigenvalue. With no loss of generality, let λ (0) be ordered as the first eigenvalue of T (0) with multiplicity m. Then T (0) can be decomposed as

T (0) = QΣR = Q 1 Q 2 Σ 1 0 0 Σ 2 R 1 R 2 , (3)
where Σ 1 is diagonal with diagonal entries as

λ (0) , R = Q -1 = r T 1 • • • r T n T , and Q = q 1 • • • q n consist
of the eigenvectors of T (0). The following lemma gives the asymptotic series for the eigenvalues of T (x) originating from λ (0) . Lemma 1 [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF] Let λ (0) be a semisimple eigenvalue of T (0) with multiplicity m. Then the corresponding eigenvalues of T (x) are analytic in x and have the form

µ i (x) = λ (0) + λ (1) i x + o(x 2 ), i = 1, . . . , m, (4) where λ (1) i , i = 1, . . . , m, are the eigenvalues of R 1 T (0)Q 1 .
Consider next the case that λ (0) is not a semisimple but repeated eigenvalue of T (0) with multiplicity m. In this case, T (0) admits a Jordan decomposition in which Σ is block diagonal with diagonal Jordan blocks, and Q 1 consists of the generalized eigenvectors associated with λ (0) . In particular,

Σ 1 =       λ (0) 1 • • • 0 0 . . . . . . 0 . . . . . . 1 0 • • • • • • λ (0)       .
The eigenvalue of T (x) can no longer be expanded in the form of ( 4), but instead as a Puiseux series.

Lemma 2 [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF] Let λ (0) be a non-semisimple eigenvalue of T (0) with multiplicity m. Then the corresponding eigenvalues of T (x) have the form

µ i (x) = λ (0) + γ (1) i 1 m x 1 m + • • • , i = 1, . . . , m, ( 5
)
where γ

(1) i

= r m T (0)q 1 . As a consequence, Lemma 1 and Lemma 2 provide a complete first-order characterization of the asymptotic expressions for the eigenvalues of T (x) near the origin.

B. Second-Order Asymptotic Expansion

The first-order asymptotic series introduced in the preceding subsection can be further developed to include terms of higher orders. In this section we derive formulas for computing the coefficients of the asymptotic series up to the second order.

Define the operator-valued function

Υ(ξ) = (T (0) -ξI) -1 , ( 6
)
which is known as the resolvent of T (0). It is obvious that the singularities of Υ(ξ) are the eigenvalues of T (0). Let λ (0) be an eigenvalue of T (0). Then Υ(ξ) can be expanded as a Laurent series at ξ = λ (0) , that is,

Υ(ξ) = -(ξ -λ (0) ) -1 P - ∞ n=1 (ξ -λ (0) ) -n-1 D n + ∞ n=0 (ξ -λ (0) ) n S n+1 , (7)
where P , D n and S n+1 are corresponding coefficient matrices. Evidently, the matrix P , known as the eigenprojection for λ (0) , can be found as

P = - 1 2πj Γ Υ(ξ)dξ = 1 2πj Γ (ξI -T (0)) -1 dξ,
where Γ is a positively-oriented closed contour enclosing λ (0) but no other eigenvalues of T (0). It was found in [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF] that

P = Q 1 R 1 .
The holomorphic part in the Laurent expansion is called the reduced resolvent of T (0) with respect to the eigenvalue λ (0) , denoted as

S(ξ) = ∞ n=0 ξ -λ (0) n S n+1 .
Let S = S(λ (0) ), namely the value of the reduced resolvent of T (0) at ξ = λ (0) . Then it is obvious that

S = S 1 = 1 2πj Γ Υ(ξ) ξ -λ (0) dξ.
Lemma 3 For any matrix T (0) decomposed in the form of [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF], where Σ 1 is in Jordan form with diagonal entries as λ (0) , the reduced resolvent at ξ = λ (0) is equal to

S = Q 0 0 0 (Σ 2 -λ (0) I) -1 R = Q 2 (Σ 2 -λ (0) I) -1 R 2 .
The following Lemma given in [START_REF] Kato | Perturbation Theory for Linear Operations[END_REF] provides the result on the second order perturbation of T (x) when all the eigenvalues of T (0) are semisimple. Lemma 4 [START_REF] Kato | Perturbation Theory for Linear Operations[END_REF] Let λ (0) be a semisimple eigenvalue of T (0), λ [START_REF] Boese | Stability with respect to the delay: On the paper of K[END_REF] i be a semisimple eigenvalue of P T (0)P with the eigenprojection P [START_REF] Boese | Stability with respect to the delay: On the paper of K[END_REF] i , that is

P (1) i = Γ i (ξI -P T (0)P ) -1 dξ, ( 8
)
where Γ i is a positively-oriented closed contour enclosing λ

(1) i but no other eigenvalues of P T (0)P . Then T (x) has d = dim P (1) i repeated eigenvalues of the form

µ ip (x) = λ (0) +xλ (1) i +x 2 µ (2) ip +o(x 2 ), p = 1, • • • , d, (9)
where µ [START_REF] Chen | Linear System Theorem and Design[END_REF] ip are the repeated eigenvalues of

P (1) i T (2) P (1) i = P (1) i (T (0) -T (0)ST (0)) P (1)
i .

The eigenvalues of P

(1) i T (2) P (1) i can be computed in a manner similar to that in the first-order analysis.

Lemma 5 Let λ

(1) i be an eigenvalue of P T (0)P . Let also T (0) be decomposed as in [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF], and R 1 T (0)Q 1 be decomposed as

R 1 T (0)Q 1 = Q (2) Σ (2) R (2) = Q (2) 1 Q (2) 2 Σ (2) 1 0 0 Σ (2) 2 R (2) 1 R (2) 2 ( 10
)
where Σ

(2)

1 is the Jordan block corresponding to the eigenvalue λ [START_REF] Boese | Stability with respect to the delay: On the paper of K[END_REF] i . Then,

P (1) i = Q 1 Q (2) 1 R (2) 1 R 1 ,
and the eigenvalues of

P (1) i T (2) P (1) i are those of the matrix R (2) R 1 T (2) Q 1 Q (2) 1 .
We note that when λ (1) i is a simple eigenvalue, then

P (1) i = q 1 r 1 .
Finally, in the case that λ

(i)
1 is a repeated non-semisimple eigenvalue of P T (0)P , we can state a result similar to Lemma 2. We state this result without proof. Lemma 6 Let λ (0) be a semisimple eigenvalue of T (0), and λ

(1) i be a non-semisimple eigenvalue of P T (0)P with multiplicity d. Then the d corresponding eigenvalues of T (x) can be expanded in the form [START_REF] Gosiewski | The effects of feedback delays on the performance of multivariable linear control systems[END_REF] where rd and q1 are the dth right and the first left eigenvectors of P T (0)P , respectively.

µ i (x) = λ (0) + xλ (1) i + rd T (2) q1 1/d x 1+ 1 d + • • • , i = 1, . . . , d,

III. STATE-SPACE MODELS

In this section we consider linear time-delay systems described by the state-space equation

ẋ(t) = A 0 x(t) + q k=1 A k x(t -kτ ), τ ≥ 0, ( 12
)
where τ is the delay parameter and A k ∈ IR n×n are the given system matrices. The characteristic quasipolynomial associated with this system is given by p(s, e -τ s ) = det sI -

q k=0 A k e -kτ s . ( 13
)
For a fixed τ ≥ 0, the system is asymptotically stable if and only if all the zeros of the quasipolynomial p(s, e -τ s ) lie in C -(see, e.g., [START_REF] Gu | Stability of Time-Delay Systems[END_REF]).

The critical delay values and the critical zeros of the above system can be computed in the following way. Define the matrices

U =        I 0 • • • 0 0 0 I • • • 0 0 . . . 0 0 • • • I 0 0 0 • • • 0 B 2q        , V =        0 I 0 • • • 0 0 0 I • • • 0 . . . 0 0 0 • • • I -B 0 -B 1 -B 2 . . . -B 2q-1        , where B m ∈ IR n 2 ×n 2 , m = 0, 1, • • • , 2q are given by B q-m = I ⊗ A T m , B q = A 0 ⊕ A T 0 , B q+m = A m ⊗ I.
The following lemma, obtained in [START_REF] Chen | A new method for computing delay margins for stability of linear delay systems[END_REF] (see also [START_REF] Fu | Stability of linear neutral timedelay systems: exact conditions via matrix pencil solutions[END_REF]), shows that the critical delay values and zeros can be determined by solving a generalized eigenvalue problem.

Lemma 7

The characteristic quasipolynomial p(s, e -τ s ) has a critical zero on the imaginary axis if and only if the following conditions are satisfied:

(i) σ(V, U ) ∩ ∂ID = ∅; (ii) For some z i ∈ σ(V, U ) ∩ ∂ID, σ q k=0 A k z k i ∩ jIR + = ∅.
The imaginary number jω i ∈ σ(

q k=0 A k z k i ),
where ω i > 0, is a critical zero. The corresponding critical delay forms the set

T (ω i ) = Log(z i ) jω i + 2π ω i > 0, = 1, 2, • • • ,
where Log(•) represents the Cauchy principal value and (.) the complex conjugate.

Note that while T (ω i ) is a countable set consisting of infinitely many elements, the number of critical zeros is finite.

In addition, since the zeros of p(s, e -τ s ) are conjugate symmetric, it suffices to consider only the critical zeros with ω i > 0. The entire range of delay values can be partitioned into intervals (τ i , τ i+1 ), and the stability is invariant in each of such intervals. Let τ * and jω * , ω * ∈ IR + , ω * = 0 be a pair of critical delay and critical zero of p(s, e -τ s ), i.e., p(jω * , e -jω * τ * ) = 0. The asymptotic behavior of the critical imaginary zero jω * can be examined by casting the problem into one of eigenvalue perturbation. Indeed, introduce a new real variable x = τ -τ * , and define

T (x) := q k=0
A k e -jω * kτ * e -jω * kx . (

Clearly, T (x) is holomorphic. Furthermore, since

det(T (0) -jω * I) = p(jω * , e -jω * τ * ) = 0,
where

T (0) = q k=0 A k e -jω * kτ * , ( 15
)
jω * is an eigenvalue of T (0). Without loss of generality, let jω * be ordered as the first eigenvalue of T (0), with multiplicity m. Let also T (0) be partitioned as in [START_REF] Chen | Asymptotic behavior of imaginary zeros of linear systems with commensurate delays[END_REF]. Then, the following asymptotic expansion can be obtained.

Theorem 1 Let jω * be a semisimple eigenvalue of T (0) with multiplicity m, and let λ

(1) i be a semisimple eigenvalue of P T (0)P with multiplicity d. Then for any τ sufficiently close to τ * , the characteristic zeros corresponding to jω * can be expanded into the power series

jω * + λ (1) i (τ -τ * ) + µ (2) ip (τ -τ * ) 2 + o (τ -τ * ) 3 , (16) with λ (1) i = λ i (R 1 T (0)Q 1 ) , i = 1, 2, • • • , m, µ (2) ip = λ p R (2) 1 R 1 (T (0) -T (0)ST (0)) Q 1 Q (2) 1 , p = 1, 2, • • • , d,
where

T (0) = - q k=1 jkω * A k e -jω * kτ * , T (0) = - q k=1
(ω * k) 2 A k e -jω * kτ * , and S is given in Lemma 3. (i) For τ sufficiently close to τ * but τ > τ * , there are m repeated characteristic zeros entering the right-half plane (or the left-half plane) if the following eigenvalues satisfy the condition

Re {λ i (R 1 T (0)Q 1 )} > 0 (< 0). (ii) If Re {λ i (R 1 T (0)Q 1 )} = 0, , i = 1, 2, • • • , m,
then for τ sufficiently close to τ * but τ > τ * , there are d repeated characteristic zeros entering the right-half plane (or the left half plane) if the eigenvalues µ

(2)

ip satisfy the condition

Re λ p R (2) 1 R 1 (T (0) -T (0)ST (0)) Q 1 Q (2) 1 > 0 (< 0).
Thus, upon finding the critical delays and critical zeros, the additional computation for the asymptotic expansion given in Theorem 1 requires only computing the eigenvalues in the expressions of λ (1) i and µ

(2) ip . In particular, when jω * is a simple critical zero, this computation can be further simplified. The following corollary is an immediate consequence of Theorem 1.

Corollary 1 Let jω * be a simple eigenvalue of T (0). Then for any τ sufficiently close to τ * but τ > τ * , jω * enters the right-half plane (or the left-half plane) if

Re {r 1 T (0)q 1 } > 0 (< 0), where q 1 and r 1 are the right and left eigenvectors associated with jω * . Additionally, if

Re {r 1 (T (0)) q 1 } = 0, then for any τ sufficiently close to τ * but τ > τ * , jω * enters the right-half plane (or the left half plane) if

Re {r 1 (T (0) -T (0)ST (0)) q 1 } > 0 (< 0), (17)
The next two results concern the cases that jω * is not a semisimple but repeated eigenvalue, or that λ (1) i is such an eigenvalue.

Theorem 2 Let jω * be a repeated eigenvalue of T (0) with multiplicity m. Suppose that jω * is not semisimple. Then for any τ sufficiently close to τ * but τ > τ * , the m characteristic zeros corresponding to jω * can be expanded into the Puiseux series

jω * + {-r m T (0)q 1 } 1 m (τ -τ * ) 1 m + • • • , Hence, for k = 0, 1, • • • , m -1

, the kth branch of the eigenvalue enters the right half plane (or the left half plane) if

cos 2kπ + π + θ m > 0 (< 0), where θ ∈ [0, 2π] is the phase angle of r m T (0)q 1 .
Theorem 3 Let jω * be a semisimple eigenvalue of T (0) with multiplicity m and λ

(1) i = λ i (R 1 T (0)Q 1 )
be a repeated eigenvalue of P T (0)P with multiplicity d. Suppose that λ

(1) i is not semisimple. Then for any τ sufficiently close to τ * but τ > τ * , the characteristic zeros corresponding to jω * can be expanded by the Puiseux series

jω * + λ (1) i (τ -τ * ) + -r d T (2) q1 1 d (τ -τ * ) 1+ 1 d + • • • , ( 18
)
where

T (2) = R (2) 1 R 1 (T (0) -T (0)ST (0)) Q 1 Q (2) 1 ,
and rd , q1 are the dth right and the first left eigenvectors of P T (0)P . Hence, if

Re{λ i (R 1 T (0)Q 1 )} = 0, then for k = 0, 1, • • • , d -1
, the kth branch of the eigenvalue enters the right half plane (or the left half plane) if

cos 2kπ + π + θ d > 0 (< 0),
where θ ∈ [0, 2π] is the phase angle of rm T (2) q1 .

Theorem 2 and Theorem 3 reveal a fundamental difference between the asymptotic behaviors of a semisimple eigenvalue and one that is not. Take Theorem 2 for example, for a repeated eigenvalue that is not semisimple, the result shows that, precluding the case that Re{r m T (0)q 1 } = 0, whether the zero will enter the right-half plane is solely determined by the multiplicity of the zero. In this case, the branches of the zero loci will almost generically enter the right-half plane when the multiplicity is greater than two. Note that in the degenerate case m = 1, Theorem 2 also reduces to Corollary 1.

IV. DIFFERENTIAL-DIFFERENCE EQUATION MODEL

Alternatively, we also consider the differential-difference equation

y (n) (t) + n-1 i=0 q k=0 a ki y (i) (t -kτ ) = 0, τ ≥ 0. ( 19
)
Systems described by ( 19) can be represented equivalently by the quasipolynomial

a s, e -τ s = q k=0 a k (s)e -kτ s , τ ≥ 0, (20)
where

a 0 (s) = s n + n-1 i=0 a 0i s i , a k (s) n-1 i=0 a ki s i , k = 1, • • • , q.
We note that a similar generalized eigenvalue-based solution (see, e.g., [START_REF] Chen | A new method for computing delay margins for stability of linear delay systems[END_REF]) exists for the critical zeros jω * of a (s, e -τ s ) at the critical delay value τ * , at which a jω * , e -jω * τ * = 0. This result is available from [START_REF] Chen | A new method for computing delay margins for stability of linear delay systems[END_REF], [START_REF] Fu | Stability of linear neutral timedelay systems: exact conditions via matrix pencil solutions[END_REF].

Lemma 8 Define H n := 0, T n := I, and for i = 0, , 1,

• • • , n -1, T i : =      a 0i 0 • • • 0 a 1i a 0i • • • 0 . . . . . . . . . . . . a q-1,i a q-2,i • • • a 0i      , H i : =      a qi a q-1,i • • • a 1i 0 a qi • • • a 2i . . . . . . . . . . . . 0 0 • • • a qi      .
Define further

P i := (j) i T i (j) i H i (-j) i H T i (-j) i T T i , i = 0, 1, • • • , n, P :=      0 I • • • 0 . . . . . . . . . . . . 0 0 • • • I -P -1 n P 0 -P -1 n P 1 • • • -P -1 n P n-1      . F (s) :=      0 1 • • • 0 . . . . . . . . . . . . 0 0 • • • 1 -a 0 (s) -a 1 (s) • • • -a q-1 (s)      , G(s) := diag (1 • • • 1 a q (s)) .
The quasipolynomial a(s, e -τ s ) has a critical zero on the imaginary axis if and only if the following conditions are satisfied:

(i) σ(P ) ∩ ∂IR + = ∅ and σ(P ) ∩ ∂IR + = {0}; (ii) For some ω i ∈ σ(P ) ∩ ∂IR + , σ(F (jω i ), G(jω i )) ∩ ∂ID = ∅.
The imaginary number jω i ∈ σ(P ), where ω i > 0, is a critical zero. The corresponding critical delay forms the set

T (ω i ) = Log(z i ) jω i + 2π ω i > 0, = 1, 2, • • • , where z i ∈ σ(F (jω i ), G(jω i )).
It is clear that Theorem 1-Theorem 3 can be directly applied to determine the asymptotic behavior of the critical zero jω * , by making use of the realization

A 0 =      0 1 • • • 0 . . . . . . . . . . . . 0 0 • • • 1 -a 00 -a 01 • • • -a 0,n-1      , ( 21
)
A k =      0 0 • • • 0 . . . . . . . . . . . . 0 0 • • • 0 -a k0 -a k1 • • • -a k,n-1      , k > 0. (22)
Additionally, with the specific structures of this realization, it is not difficult to derive that the eigenvalues of T (0) are also the roots of quasipolynomial a s, e -jω * τ * . This makes it possible to state the results directly in terms of the quasipolynomial a (s, e -τ s ). We first show that with a realization given by the above companion form, T (0) cannot have semisimple repeated eigenvalues.

Lemma 9 Any companion matrix

A =      0 1 • • • 0 . . . . . . . . . . . . 0 0 • • • 1 -α 0 -α 1 • • • -α n-1      cannot have semisimple repeated eigenvalues.
The implication of Lemma 9 is rather clear. With the above realization of A k , the matrix T (0) is in the companion form with

α i = q k=0 a ki e -jkω * τ * , i = 0, 1 • • • , n -1.
As such, when applying the results in Section III, one can preclude the case that T (0) has semisimple repeated eigenvalues; instead, the critical zero jω * is either a simple eigenvalue, or a repeated, non-semisimple eigenvalue of T (0). Theorem 4 and Corollary 2 given below describe the asymptotic behavior of imaginary poles in these two cases, in terms of a (s, e -τ s ).

Theorem 4 Let jω * be a repeated zero of a s, e -jω * τ * with multiplicity m. Then for any τ sufficiently close to τ * but τ > τ * , the m zeros corresponding to jω * can be expanded into the Puiseux series

jω * +   -m! d a(jω * , e -jω * τ ) dτ | τ =τ * d m a(s, e -jω * τ * ) ds m | s=jω *   1 m (τ -τ * ) 1 m + • • • ,
If jω * is a simple zero of a s, e -jω * τ * , the following corollary is an immediate consequence of Theorem 4. It is useful to note that in the case of a simple critical zero, results similar to Corollary 2 have been reported in [START_REF] Boese | Stability with respect to the delay: On the paper of K[END_REF], [START_REF] Cooke | Discrete delay, distributed delay and stability switches[END_REF], [START_REF] Cooke | On zeroes of some transcendental equations[END_REF], [START_REF] Hertz | Simplified analytical stability test for systems with commensurate time delays[END_REF], [START_REF] Walton | Direct method for TDS stability analysis[END_REF], while without much technical deliberation. Our derivation herein thus provides an independent, rigorous justification to some of the previously available results.

Theorem 4 and Corollary 2 correspond to the firstorder analysis of the quasipolynomial. When these results become inapplicable, the following second-order result can be employed.

V. AN ILLUSTRATIVE EXAMPLE

In this section we use an example to demonstrate the second-order zero-crossing conditions.

Example Consider the following second-order delay system1 [START_REF] Hertz | Simplified analytical stability test for systems with commensurate time delays[END_REF] (see also [START_REF] Gosiewski | The effects of feedback delays on the performance of multivariable linear control systems[END_REF]):

ẋ1 (t) ẋ2 (t) = 0 1 -1 -1 x 1 (t) x 2 (t) + 0 0 0 -1 x 1 (t -τ ) x 2 (t -τ ) . (23)
The corresponding characteristic quasipolynomial is given by

p(s, e -τ s ) = s 2 + s + se -sτ + 1. (24)
By a direct application of Lemma 7 or Lemma 8, a pair of critical zero and delay are found to be jω * = j and τ * = π. Following the state-space computations in Theorem 1 and Corollary 1, we first find

T (0) = 0 1 -1 -1 + 0 0 0 -1 e -jπ = 0 1 -1 0 ,
which can be decomposed as

T (0) = 1 1 j -j j 0 0 -j 1/2 -(1/2)j 1/2 (1/2)j . Furthermore, T (0) = 0 0 0 -j , T (0) = 0 0 0 -1 .
It follows instantly that r 1 T (0)q 1 = -(1/2)j. As such, it is necessary to invoke the second-order condition. For this purpose, we find S = 1 -j (-j -j) -1 1/2 (1/2)j = -1 4 j -1 1 j .

A trivial calculation shows that r 1 (T (0) -T (0)ST (0)) q 1 = -

1 2 1 + 1 4 j .
This means that the zeros ±j enter the left-half plane. Since the system is stable for τ ∈ [0, π), we conclude that the critical zeros ±j do not cross the imaginary axis. We may also use Corollary which, according to Theorem 5, indicates that the imaginary zeros ±j enter the left half plane, resulting in the same conclusion based on the state-space computation. Figure 1 confirms indeed that the zeros only "touch" the imaginary axis, but do not cross it.

VI. CONCLUDING REMARKS

In this paper we have studied the asymptotic behavior of the critical characteristic zeros of linear time-delay systems with commensurate delays, for systems posed both in statespace form and as a quasipolynomial. We have shown that in both cases the asymptotic behavior of the critical zeros on the imaginary axis can be characterized in a simple manner, either by computing the eigenvalues of a constant matrix, or by computing the derivatives of the quasipolynomial. We demonstrated that when the first-order analysis ceases to be conclusive, the second-order condition can be used to effect.

Corollary 2

 2 Let jω * be a simple zero of a s, e -jω * τ * . Then for any τ sufficiently close to τ * but τ > τ * , jω * enters the right-half plane (or the left-half plane) if Re    d a(jω * , e -jω * τ ) dτ | τ =τ * d a(s, e -jω * τ * ) ds

Theorem 5 ∂ 2

 52 Let jω * be a simple zero of a s, e -jω * τ * . If Re    d a(jω * , e -jω * τ ) dτ | τ =τ * d a(s, e -jω * τ * ) ds | s=jω *    = 0, then for any τ sufficiently close to τ * but τ > τ * , jω * enters the right-half plane (or the left-half plane) if Re    d 2 a(jω * , e -jω * τ ) dτ 2 d a(s, e -jω * τ * ) ds -  d a(jω * , e -jω * τ ) dτ d a(s, e -jω * τ * ) ds a(s, e -jω * τ ) ∂s∂τ d a(s, e -jω * τ * ) ds at s = jω * , τ = τ * .

∂ 2

 2 2 and Theorem 5 to check the behavior of the imaginary zero jω * = j. In this vein, we first compute dp(jω * , e -jω * τ ) dττ =τ * = -(jω * ) 2 e -jω * τ | τ =π = -1, dp(s, e -jω * τ * ) ds s=jω * = 2s + 1 + e -jω * τ * | s=j = 2j.This gives rise tod p(jω * , e -jω * τ ) dτ | τ =τ * d p(s, e -jω * τ * )ds 2 is rendered inconclusive. To proceed to the second-order condition, we findd 2 p(jω * , e -jω * τ ) dτ 2 τ =τ * = (jω * ) 3 e -jω * τ | τ =π = j, ∂ 2 p(s, e -jω * τ ) ∂s∂τ s=jω * ,τ =τ * = -jω * e -jω * τ | τ =π = j.It then follows that   d 2 p(jω * , e -jω * τ ) dτ 2 d p(s, e -jω * τ * ) ds jω * , e -jω * τ ) dτ d p(s, e -jω * τ * ) ds p(s, e -jω * τ ) ∂s∂τ d p(s, e -jω * τ * ) ds

Fig. 1 .

 1 Fig. 1. Eigenvalues as functions of the delay τ

Note that the state-space form of this system was incorrectly given in[START_REF] Hertz | Simplified analytical stability test for systems with commensurate time delays[END_REF].

This research is supported in part by CNSF/China, CNRS/France, and NSF/USA.


M Bourrel

P Aumont

P L Azemar

Thuy Vu

B Le Toan

Lubac

Influence of recent

Keywords: Surface water storage, Rainfall, ENSO, PDO, Monsoon, E V I, Mekong Basin, Multi-satellite

  L'archive ouverte pluridisciplinaire

Introduction

Climate variability and human population growth have led to inland waters being identified as some of the most severely endangered ecosystems [START_REF] Rosenzweig | Assessment of observed changes and responses in natural and managed systems[END_REF]. Lakes and reservoirs are widely distributed around the world and act as integrators of changes occurring in their watershed, the surrounding landscape and the atmosphere [START_REF] Schindler | Lakes as sentinels and integrators for the effects of climate change on watersheds, airsheds, and landscapes[END_REF]Williamson et al., 2009a). As their physical, chemical, and biological properties respond rapidly to climate-induced changes, they offer a valuable source of information for assessing and monitoring the effects of climate change (e.g., [START_REF] Magnuson | Historical Trends in Lake and River Ice Cover in the Northern Hemisphere[END_REF][START_REF] Lemoalle | Recent changes in Lake Chad: Observations, simulations and management options (1973-2011)[END_REF]. Lakes and reservoirs are commonly seen as sentinels of climate change as their physical, chemical, and biological responses to climate variations integrate the changes occurring in their catchment [START_REF] Adrian | Lakes as sentinels of climate change[END_REF][START_REF] Schindler | Lakes as sentinels and integrators for the effects of climate change on watersheds, airsheds, and landscapes[END_REF]C. E. Williamson et al., 2009).

Tonle Sap Lake in Cambodia is the largest lake and wetland ecosystem in the Mekong River Basin (MRC, 2005). It is considered as the "heart of the Lower Mekong" due to the exceptional biodiversity of its ecosystem [START_REF] Campbell | Species diversity and ecology of Tonle Sap Great Lake, Cambodia[END_REF] and is also one of the most productive ecosystems in the world [START_REF] Lamberts | The Tonle Sap Lake as a Productive Ecosystem[END_REF] with its high productivity being linked to the flood pulse of the Mekong River [START_REF] Kummu | Ecosystem Management of the Tonle Sap Lake: An Integrated Modelling Approach[END_REF]. The hydrological cycle in the Mekong River Basin is governed by the monsoon regime, itself being strongly influenced by El Niño Southern Oscillation (ENSO). El Niño events are responsible for a decrease in the rainfall, discharge and annual flood period [START_REF] Räsänen | Spatiotemporal influences of ENSO on precipitation and flood pulse in the Mekong River Basin[END_REF]. Conversely, La Niña increases the rainfall, discharge and annual flood period. The combination of climaterelated and anthropogenic (hydropower, irrigation) changes in the Mekong River Basin are expected to cause large modifications to the flood pulse of the Tonle Sap [START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF][START_REF] Västilä | Modelling climate change impacts on the flood pulse in the Lower Mekong floodplains[END_REF]. The frequency and intensity of drought and flood events in the region have increased over the past few decades, impacting irrigated fields [START_REF] Cruz | Asia. Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF][START_REF] Kamoshita | Field level damage of deepwater rice by the 2011 Southeast Asian Flood in a flood plain of Tonle Sap Lake, Northwest Cambodia[END_REF]. Dam construction and operations in the Chinese Upper

Mekong Basin are also altering the hydrological regime in the Lower Mekong Basin [START_REF] Stone | Severe Drought Puts Spotlight on Chinese Dams[END_REF][START_REF] Lu | Observed changes in the water flow at Chiang Saen in the lower Mekong: Impacts of Chinese dams?[END_REF]. Development of new hydropower dams in the Chinese Upper Mekong

Basin combined with possible development of thirty-eight new dams in the Vietnamese

Lower Mekong Basin will have a major impact on the hydrological regime of this region (see [START_REF] Lu | Observed changes in the water flow at Chiang Saen in the lower Mekong: Impacts of Chinese dams?[END_REF] for the impacts on the Chinese Upper Mekong Basin; [START_REF] Arias | Impacts of hydropower and climate change on drivers of ecological productivity of Southeast Asia's most important wetland[END_REF] for the impacts on the Vietnamese Lower Mekong Basin).

Remotely sensed observations offer a unique opportunity to continuously monitor floods in large watersheds. Previous studies have determined spatio-temporal flooding extent by applying a threshold approach to land surface reflectances, or spectral indexes from the MODerate resolution Imaging Sensor (MODIS) radiometer [START_REF] Sakamoto | Detecting temporal changes in the extent of annual flooding within the Cambodia and the Vietnamese Mekong Delta from MODIS time-series imagery[END_REF][START_REF] Arias | Quantifying changes in flooding and habitats in the Tonle Sap Lake (Cambodia) caused by water infrastructure development and climate change in the Mekong Basin[END_REF][START_REF] Siev | Seasonal Changes in the Inundation Area and Water Volume of the Tonle Sap River and Its Floodplain[END_REF][START_REF] Fayne | Flood mapping in the lower Mekong River Basin using daily MODIS observations[END_REF] in the Tonle Sap watershed. Gravity Recovery and Climate Experiment (GRACE) products have recently been used along with MODIS images to estimate total water storage and identify flood events in the Tonle Sap watershed [START_REF] Tangdamrongsub | Assessing total water storage and identifying flood events over Tonlé Sap basin in Cambodia using GRACE and MODIS satellite observations combined with hydrological models[END_REF]. Estimates of surface water volume changes in the Tonle Sap watershed are sparse and have only been for short time-periods (e.g., 1997[START_REF]Overview of the Hydrology of the Mekong Basin[END_REF][START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF]and 2003[START_REF]Overview of the Hydrology of the Mekong Basin[END_REF][START_REF] Siev | Seasonal Changes in the Inundation Area and Water Volume of the Tonle Sap River and Its Floodplain[END_REF].

The purpose of this study is to estimate Tonle Sap lake storage variations over the two last decades and analyze the link between its annual variations with atmospheric forcing, especially climatic oscillations. Surface water storage variations of the lake are computed over the 1993-2017 time span, through a synergistic use of flood extent derived from MODIS land surface reflectances (see sections 3.1 and 4.1) and altimetry-based water levels (see sections 3.2 and 4.2). During these two decades, three extreme droughts (1998, 2015, and 2016) and three exceptional floods (2000, 2002, and 2011) occurred in the Tonle Sap watershed. In the following, relationships with climate variability are discussed through comparisons with climate indices (sections 3.3 and 5.2). A first impact assessment of the extreme drought of 2015 on the Tonle Sap ecosystem is also presented.

Study area

Ninety-five percent of the Tonle Sap Lake drainage area is located in Cambodia, with the remaining five percent in Thailand. It is the largest lake and wetland in the Mekong River

Basin and the largest lake in Southeast Asia (MRC, 2005;[START_REF] Campbell | Chapter 10 -Tonle Sap Lake, the Heart of the Lower Mekong[END_REF]. Its whole drainage basin occupies an area of 85,800 km² (i.e., 11% of the whole surface of the Mekong basin, MRC 2003) (Fig. 1). The lake-floodplain system extends from 2,400 km² during the low water period to 16,000 km² at peak flood during wet years [START_REF] Lim | Diversity and spatial distribution of freshwater fish in Great Lake and Tonle Sap river (Cambodia, Southeast Asia)[END_REF][START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF].

The catchment incorporates five different land cover classes, as described by Arias et al.,

 Open water (under water from 9 to 12 months a year and covering an area of 3,027 km²), Tonle Sap Lake is an important ecosystem of Southeast Asia due to its exceptional biodiversity [START_REF] Campbell | Species diversity and ecology of Tonle Sap Great Lake, Cambodia[END_REF]. It has been internationally recognized through three Biosphere Reserve core areas under the United Nations Educational, Scientific and Cultural Organization (UNESCO) Programme on Man and the Biosphere (UNESCO, 2006), as well as one Ramsar site under the Ramsar Convention on wetlands (see www.ramsar.org). The Tonle Sap lake ecosystem also supports important fisheries and aquaculture activities. These industries sustain over one million people [START_REF] Rainboth | Fishes of the Cambodian Mekong[END_REF][START_REF] Sverdrup-Jensen | Fisheries in the Lower Mekong Basin: Status and Perspectives[END_REF][START_REF] Hortle | Consumption and the Yield of Fish and Other Aquatic Animals from the Lower Mekong Basin[END_REF].


The annual monomodal flood pulse in Tonle Sap Lake is governed by the flow in the Mekong River [START_REF] Junk | The comparative biodiversity of seven globally important wetlands: a synthesis[END_REF][START_REF] Lamberts | The Tonle Sap Lake as a Productive Ecosystem[END_REF]. The Tonle Sap River connects the lake to the Mekong River. In the May to September wet season its flow direction is inversed due to water level increase in the Mekong River [START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF]. Annual average inflow and outflow are 79.0 and 78.6 km 3 , respectively. 57% of the inflow comes from the Mekong mainstem, either through Tonle Sap River discharge (52%) or by overland flooding (5%)

(MRCS/WUP-FIN, 2006).
The contributions of Tonle Sap Lake tributaries and direct rainfall into the lake are estimated at around 30% and 13%, respectively (Kummu et al., 2014). Water receding from the lake mostly flows into the Mekong River (88%) or evaporates (12%) (Kummu et al., 2014).

Datasets

Surface reflectances from MODIS

The MODIS spectroradiometer is part of the payload of the Aqua and Terra satellites (launched in 2002 and 1999, respectively). The MODIS sensor acquires radiances in 36 spectral bands. In this study, the MOD09A1 product (8-day binned level 3, version 6) derived from Terra raw radiance measurements were downloaded from the United States Geological Survey (USGS) Earthexplorer website (http://earthexplorer.usgs.gov). It consists in gridded surface reflectances acquired in 7 bands from 620 nm (Red) to 2155 nm (SWIR2) at a resolution of 500 m and corrected for atmospheric effect. This product is obtained by combining the best surface reflectance data of every pixel acquired during an 8-day period for each wavelength. Each MODIS tile covers an area of 1200 km by 1200 km.

In this study, 815 composite images from February 2000 to December 2017 were used.

Radar altimetry data

Our radar altimetry data comes from the acquisitions on their nominal orbit of the following missions:

Topex-Poseidon (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002), Jason-1 (2002Jason-1 ( -2008)), Jason-2 (2008Jason-2 ( -2016)), Jason-3 (since 2016), ERS-2 (1995ERS-2 ( -2003)), ENVISAT (2003[START_REF] Stone | Severe Drought Puts Spotlight on Chinese Dams[END_REF]), and SARAL (2013-2016). Topex/Poseidon and Jason-1 to -3 missions are on the same 10-day repeat period orbit, whereas ERS- Topography Mission data ("The Shuttle Radar Topography Mission," 2007) for the surrounding areas to complete the DBM (see ("Water balance analysis for the Tonle Sap Lake-floodplain system," 2014)for more details).

SRTM Digital Elevation Model (DEM)

The Shuttle Radar Topography Mission (SRTM) DEM was derived from the C and X Synthetic Aperture Radar (SAR) images acquired for interferometry purposes on board the space shuttle Endeavour from 11-22 February 2000 [START_REF] Farr | The Shuttle Radar Topography Mission[END_REF]. In this study, the SRTM

MODIS-based Global Land Cover Climatology

This climatology is based on the Collection 5.1 MODIS Land Cover Type (MCD12Q1) product at 500 m spatial resolution and obtained from analysis of MODIS images acquired between 2001 and 2010.

A re-analysis attributing a confidence score for each land cover type at each pixel was performed, as substantial interannual variability in land cover types was found to affect 40% of the pixels. The final land cover classification was obtained by choosing the type with the highest confidence score between 2001 and 2010 for each pixel [START_REF] Broxton | A Global Land Cover Climatology Using MODIS Data[END_REF]. This dataset is available from USGS at: https://landcover.usgs.gov/global_climatology.php.

Monsoon Indices

The Indian Monsoon Index (IMI) and the Western North Pacific Monsoon Index (WNPMI) are defined as the difference of 850 hPa zonal wind over the areas (5°N-15°N, 40°E-80°E) and (20°N-

Methods

The lake and floodplain surface extent (section 4.1) and elevation (section 4.2) are first estimated using MODIS and radar altimetry measurements, respectively. They are then combined to compute the lake storage change (section 4.3). Finally, a rainfall standardized anomaly index is computed (section 4.4) to investigate the link between this storage change and rainfall.

Surface water extent

The approach developed by [START_REF] Sakamoto | Detecting temporal changes in the extent of annual flooding within the Cambodia and the Vietnamese Mekong Delta from MODIS time-series imagery[END_REF] and simplified by [START_REF] Normandin | Quantification of surface water volume changes in the Mackenzie Delta using satellite multi-mission data[END_REF] to monitor flood extent in the Mekong Basin was applied to the Tonle Sap drainage area (processing steps are presented in Figure S1 in the supplementary material section). It is based on the thresholding of the Enhanced Vegetation Index (EVI), the Land Surface Water Index (LSWI), and the Difference

Value between EVI and LSWI (DVEL) to determine the status (non-flooded, mixture, flooded, permanent water body) of any pixel in an 8-day composite MODIS image of surface reflectance.

These two spectral indexes were selected as the Short-wave infrared (SWIR) band is very sensitive to soil and vegetation moisture content and because the Near infrared (NIR) band presents low surface reflectance values over water bodies and high values over densely vegetated areas, whereas surface reflectance values in the red band are higher for vegetation than over water bodies. Thus, these two indexes are well suited to discriminate between water and vegetation.

The two indexes used in this approach are defined as follows [START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF]:

(1) (2)
where ρ NIR is the surface reflectance value in the near infrared (841-875 nm, band 2), ρ R is the surface reflectance value in the red (621-670 nm, band 1), ρ B the surface reflectance value in the blue (459-479 nm, band 3), and ρ SWIR the surface reflectance in the short-wave infrared (1628-1652 nm, band 6).

Following [START_REF] Sakamoto | Detecting temporal changes in the extent of annual flooding within the Cambodia and the Vietnamese Mekong Delta from MODIS time-series imagery[END_REF], the approach summarized in Figure S1 for [START_REF] Normandin | Quantification of surface water volume changes in the Mackenzie Delta using satellite multi-mission data[END_REF] was applied to the series of MODIS images. First, the pixels identified as cloudcovered were filtered out (ρ B ≥ 0.2). Then, a linear interpolation was performed to limit information loss over the study area. Radar altimetry was initially developed to provide accurate measurements of the sea surface topography and is now commonly used for the monitoring of inland water levels (see [START_REF] Crétaux | Hydrological applications of satellite altimetry: rivers, lakes, man-made reservoirs, inundated areas[END_REF]2017 for a recent review). Variations of altimeter height from one cycle to another can be associated with water level changes.

In this study, we used the Multi-mission Altimetry Processing Software (MAPS) for processing altimetry data over land and ocean (e.g., (Frappart et al., 2015c); (Frappart et al., 2015a); [START_REF] Biancamaria | Satellite radar altimetry water elevations performance over a 200 m wide river: Evaluation over the Garonne River[END_REF]; [START_REF] Vu | Multi-Satellite Altimeter Validation along the French Atlantic Coast in the Southern Bay of Biscay from ERS-2 to SARAL[END_REF], that allows to build time-series of water levels at a virtual station (intersection between satellite track and water body), after a refined selection of the valid altimeter data. Data processing is composed of four main steps:

i)

a rough delineation of the cross-section between the altimeter tracks and Tonle Sap using Google Earth (Fig. 1b);

ii)

loading of the altimetry over the study area and computation of altimeter heights from raw data contained in the GDR;

iii) a refined selection of the valid altimetry data through visual inspection;

iv)

computation of the water level time-series as the median of the selected water levels for each cycle.

A detailed description of altimetry data processing using MAPS can be found in (Frappart et al., 2015b) Measurement bias between missions and difference in satellite track locations were computed and removed from each time series to obtain a consistent time series over the lake between all missions.

The final water elevation time series (in meters wrt EGM08 geoid) is computed with reference to ENVISAT track 107 (the most eastward ENVISAT track on the lake, see Fig. 1b).

Specific processing of Topex-Poseidon and Jason time series

As shown in Figure 1b, Topex-Poseidon and Jason track over Tonle Sap Lake is very close to its outlet. Consequently, and as visually observed on Landsat images, during low flow period the observed area is representative of the Tonle Sap River, which connects the lake to the Mekong River, whereas during high flow period it fully includes the lake extent. Correlation between Jason-2 time series and ENVISAT over the period in common (more than 2 years) still remains high (0.9936), with an RMSE of 0.23 m, after removing the common bias between the two time series. For low water elevations (below 5.5 m), Jason-2 water elevations rise before those of ENVISAT and decrease later.

This could be due to hydraulic effects: water elevations vary differently in the lake and outlet river because of bathymetry changes. To account for this observation in a simple way, we arbitrarily considered that difference between the two time series can be corrected applying a time lag to Jason-2 time series, depending on the water elevations, to get a better match with ENVISAT time series. After different tests on time lags and water elevation thresholds, the best correlation (0.9986) and lowest RMSE between ENVISAT and Jason-2 (0.11 m) are obtained when Jason-2 water elevations below 5.5 m are delayed by 9 days, whereas water elevations above this threshold are delayed by 2 days.

More information on this issue is shown in Figures S2 to S4 in the supplementary information section. The Topex-Poseidon, Jason-1 and Jason-3 time series used in this study are from the same track and therefore the same processing has been applied (after removing the bias between these time series and Jason-2 time series).

Surface water volume variations

MODIS-based flood extents were combined with altimetry-derived water levels to estimate lake volume variations. Due to the extent of flooding in the watershed, our analysis was limited to the pixels which are i) under water filling SRTM DEM with the water stage of the lake (see the supplementary information, section 2, for method) and ii) connected to the lake using a four-connected neighborhood. As flood extent is derived from a shorter time period than water stage (1993-2017 instead of 2000-2017), a polynomial relationship between water stage (h) and lake surface (S) was fitted to estimate surface water extent during the whole altimetry period using Eq. 3:

( ) ∑ (3)
where a i is the coefficient of degree i of the polynomial and n the maximum degree of the polynomial. A confidence interval of 95%, a RSME of 74 km² and R² of 0.99 were found using a second degree polynomial with the following values: a 0 =6013.0 km², a 1 =-1518.0 km and a 2 =-176.9. As in [START_REF] Baup | Combining high-resolution satellite images and altimetry to estimate the volume of small lakes[END_REF] or [START_REF] Crétaux | Global surveys of reservoirs and lakes from satellites and regional application to the Syrdarya river basin[END_REF]2016), the volume variations of the lake (ΔV), between two consecutive time steps (t-1 and t) were estimated using the assumption of a regular morphology of the lake and a pyramidal shape:

( ( ) ( )( ( ) ( ) √ ( ) ( ))
(4)

Standardized Precipitation Anomaly Index

The Standardized Precipitation Anomaly Index (SPAI) was defined by [START_REF] Lamb | Persistence of Subsaharan drought[END_REF] to study interannual rainfall variations:

( ) ∑ ( ) 〈 ( )〉 ( ( )) (5)
where R(i,t) is the annual rainfall at TRMM 3B43v7 gridpoint i for hydrological year t, 〈 ( )〉 and ( ( )) are the annual average and standard deviation over the reference period T at gridpoint i, respectively, and n is the total number of gridpoints in the study area. The hydrological year is defined between November of year t-1 and October of year t.

Results and discussion

Validation of surface water extent and storage

Variations of lake surface extent were estimated between February 2000 and December 2017

by applying the threshold approach to MODIS images. The average number of days per year with missing data due to cloud presence were computed for each pixel, as well as the associated standard deviation over the whole study period (see Fig. months, but can also occur in the eastern part of the lake and along its north to south flowing tributary in the eastern part of the basin. Longer inundation, lasting between 4 and 5 months but of smaller extent, can be observed in the lake's outlet and also in the northeastern part (Fig. 2a). These inundations are characterized by a large interannual variability with standard deviations between 40 and 50 days (Fig. 2b). In 2011, important inundations (>100 days)

were observed in all eastern and northwestern parts of the lake (Fig. 2c), whereas almost no flooding was detected in 2015 (Fig. 2d). In addition, considering the mixed pixels (i.e., the pixels partly inundated), the patterns are similar to those presented in Fig. 2, but with a larger extent and longer presence (Fig. 3). Floods lasting up to 6 months were detected in the western parts of the lake (Figs. 3a and3c), and are in good agreement with results from previous studies [START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF][START_REF] Arias | Quantifying changes in flooding and habitats in the Tonle Sap Lake (Cambodia) caused by water infrastructure development and climate change in the Mekong Basin[END_REF]Kummu et al., 2014). Large differences were obtained between the flood extent using only the flooded pixels and using both the flooded and mixed pixels in terms of maximum extent and flood duration inside the maximum inundation extent defined filling SRTM DEM with altimetry-based water stage (Fig. S6). Even if SRTM DEM exhibits a high frequency noise (variations of 1 or 2 m from one pixel to another, see Fig. S7), large-scale elevation variations are well retrieved considering the minimum elevation in the MODIS pixels (Fig. S8) and do not prevent its use for determining temporal variations of maximum flood extent. In the following, surface water extents were obtained by considering both the flooded and mixed pixels which are i) under water filling SRTM DEM with the water stage of the lake and ii) connected to the lake using a four-connected neighborhood. good agreement was also found between MODIS-based and [START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF] surface water extent (R=0.89, RMSE=982 km² and RRMSE=9.9 %). If the very good agreement between in situ and altimetry-based water levels can confirm observations from Topex/Poseidon in a previous study with a RMSE of 0.23 cm (Frappart et al., 2006b), the small discrepancies observed in surface water extent and volume can be attributed to differences of approach in [START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF] and our own study. The composite DBM including old measurements of the floodplain's topography is liable to introduce some errors to [START_REF] Kummu | Impact of the Mekong River Flow Alteration on the Tonle Sap Flood Pulse[END_REF] approach, whereas the delineation of inundation extent using images at the moderate resolution of 500 m and the SRTM DEM is responsible for overestimation of the flooded area due to the presence of mixed pixels. Nevertheless, both approaches present very consistent results and can be used to estimate water balance in the drainage area.

Over the study period, interannual variations of surface water volume were computed from monthly lake water volume with the removal of monthly climatological averages. The range is between -35.99 km 3 in 2015 to 35.53 km 3 in 2000 (Fig. 5a). The variations were compared with SPAI (see section 4.4) from TRMM/TMPA 3B43 averaged over the Tonle Sap However, high positive anomalies (~2) were observed in the dry season and contributed to the excess in surface water storage for the same year (see Fig. 5d). This was also observed in the year 2000, which was the largest flood between 1995 and 2016. Inverse conclusions can be made for the large droughts of 1998, 2010 and 2015 with SPAI between -1.5 and -2 (Fig. 5c)

during the wet season. Lower rainfall during the dry season seems to have a low impact as small anomalies were observed during the lowest drought in 2015 (Fig. 5d). The 2002 flood was higher than that of 2001 in spite of a much lower positive SPAI. This was due to a higher than usual lake level at the beginning of the 2002 hydrological year. Conversely, low rainfall deficit was responsible for low water surface storage in 2016 due to the very low lake level at the end of the 2015 hydrological year. For a comparison, the anomaly of lake water volume was 6.8 km 3 in December 2001 and -21.9 km 3 in December 2015. In addition, daily in-situ data of water levels recorded at the Prek Kdam station in the Tonle Sap River (which joins the Mekong River to Tonle Sap Lake) are displayed on the Mekong River Commission website during both dry and flood seasons since 2013 (ffw.mrcmekong.org/historical_rec.htm). A visual inspection showed very similar water stage amplitude between in situ and altimetry time series. Extreme events were defined by considering deviations larger (lower) than two thirds of the maximum (minimum) estimated over the observation period (i.e., around ±24 km 3 ). Four positive and three negative extrema can be observed. For positive extrema:

between July and November 2000, peaking at 36.36 km 3 in September 2000; 24.22 km 3 in September 2001; in September and October 2002, peaking at 35.37 km 3 in September 2002;

and in September and October 2011, peaking at 30.77 km 3 in October 2011. For negative extrema: between September and November 1998, peaking at -31.42 km 3 in October 1998;

between September and November 2015, peaking at -36.07 km 3 in October 2015; and -24.22 km 3 in September 2011 (Fig. 6a). In the following, we compare these exceptional droughts and floods to climatic indexes to see if they are related to climate variability.

Surface water storage and climate variability

Positive anomalies of surface water volume in the Tonle Sap watershed in September-October correspond to negative values of ENSO3.4 index defined by Rayner (2003) (Fig. 6b). On the contrary, negative anomalies of surface water volume in the watershed correspond to positive values of ENSO 3.4. Larger anomalies occur when ENSO index and the Pacific Decadal Oscillation (PDO) index (as defined by [START_REF] Mantua | A Pacific Interdecadal Climate Oscillation with Impacts on Salmon Production[END_REF]) exhibit large deviations and are in phase. Large positive (respectively negative) anomalies of the lake volume are obtained for moderate (between 1 and 1.5 °C) to very strong (above 2°C) negative (respectively positive) values of both ENSO (El Niño, respectively La Niña, Fig. 6b) and PDO (Fig. 6c) during the previous winter. This phenomenon is caused by a significant intensification of the negative relationship between ENSO and the monsoon when ENSO and PDO are in phase. It was previously described regarding the East Asian monsoon [START_REF] Kim | Combined effect of El Ni?o-Southern Oscillation and Pacific Decadal Oscillation on the East Asian winter monsoon[END_REF]. Positive values of ENSO and PDO are associated with anomalous warm temperature patterns over East Asia that cause rainfall deficit the following summer. Negative values of ENSO and PDO are linked with anomalous cold temperature patterns, leading to excess rainfall the following summer.

On the contrary, when ENSO and PDO are out of phase, interannual variations are lower. The super strong ENSO event in 2015/2016 was an exception, as it did not result in huge rainfall in either the Tonle Sap watershed or the Mekong Basin. This was most likely due to its rapid decay in boreal spring 2016 causing a weakening of the western Pacific subtropical anticyclone [START_REF] Liu | Why was the western Pacific subtropical anticyclone weaker in late summer after the 2015/2016 super El Niño?: WEAKENED WPSA FOLLOWING THE 2015/2016[END_REF]. Correlation between annual amplitude of volume variations and climatic indices were also estimated, and several time-periods were tested. Stronger correlations were found when the ENSO index was averaged over December to February (DJF) as in [START_REF] Räsänen | Spatiotemporal influences of ENSO on precipitation and flood pulse in the Mekong River Basin[END_REF] and when the PDO index was averaged over November to March (NDJFM) as in [START_REF] Newman | ENSO-Forced Variability of the Pacific Decadal Oscillation[END_REF]. For ENSO DJF , R=-0.70 with a pvalue equal to 10 -4 , whereas for PDO NDJFM , R=-0.68 with a p-value equal to 2.10 -4 . This confirms earlier findings about relationships in the Mekong Basin between rainfall, discharge, ENSO [START_REF] Räsänen | Spatiotemporal influences of ENSO on precipitation and flood pulse in the Mekong River Basin[END_REF] and PDO [START_REF] Delgado | A climate-flood link for the lower Mekong River[END_REF]). An even higher correlation was found between annual anomalies of surface water storage and ENSO 3.4 index than between annual flow and ENSO 3.4 [START_REF] Räsänen | Spatiotemporal influences of ENSO on precipitation and flood pulse in the Mekong River Basin[END_REF] and this study provides a quantification of the relationship between flood pulse and PDO in the Mekong

Basin. The optimal least-squares linear combination of both indices provides a correlation of -0.75. The coefficients of this linear combination are equal to 1.0 and 0.84 for ENSO 3.4 and PDO, respectively. Time series of annual amplitude of lake volume (Fig. 7a), of the ENSO 3.4 index averaged over DJF (Fig. 7b) and of the PDO index averaged over NDJFM (Fig. 7c) We propose to quantify the monsoonal impact on lake storage, as the flood pulse in the Mekong Basin is driven by the monsoon. The Mekong Basin is under the influence of both the Indian and the Western North Pacific Monsoons [START_REF] Delgado | A climate-flood link for the lower Mekong River[END_REF]. Correlation between anomalies of annual surface water storage and Western North Pacific Monsoons (WNPMI, [START_REF] Wang | Interannual Variability of the Asian Summer Monsoon: Contrasts between the Indian and the Western North Pacific-East Asian Monsoons*[END_REF] and Indian Monsoon (IMI, Wang and Fan 1999) provides better results (R=0.55, p-value=0.006, and R=0.49, p-value=0.018) for July and August, respectively. This is similar to the relationship between water discharge and WNPMI in the downstream Mekong Basin discovered by [START_REF] Delgado | A climate-flood link for the lower Mekong River[END_REF].

illustrate

5.3) Impact of the 2015 extreme drought on vegetation activity

Vegetation activity is highly related to flood duration [START_REF] Arias | Quantifying changes in flooding and habitats in the Tonle Sap Lake (Cambodia) caused by water infrastructure development and climate change in the Mekong Basin[END_REF]. In the following, impact of the 2015 drought on the Tonle Sap ecosystem is analyzed using MODIS-derived EVI time series and a land cover map. The 1998 and 2015 extreme droughts were the worst observed during the study period, but MODIS images are onlyavailable since February 2000.

Due to this, our analysis is limited to the 2015 drought. Figure 8a shows land cover types from the 500 m MODIS-based Global Land Cover Climatology [START_REF] Broxton | A Global Land Cover Climatology Using MODIS Data[END_REF] anomalies of -0.11, -0.13 and -0.06 (Fig. 8c), corresponding to an annual variation of amplitude around 0.10, 0.25 and 0.15, respectively. The same timing is observed for Savannas, Grasslands and Permanent Wetland with EVI anomalies of -0.10, -0.05 and -0.12 (Fig. S9b), corresponding to an annual variation of amplitude around 0.25, 0.10 and 0.30, respectively. They reached their lowest values in May 2016 for Mixed Forests and Croplands with EVI anomalies of -0.15 and -0.04 (Fig. S9b), corresponding to annual variation of amplitude around 0.10 and 0.15, respectively. The spatial distribution of anomalies of EVI for April and May are presented in Figures 9a and9b, respectively. They are quite similar to lowest anomaly values in April 2016. The largest negative anomalies are located on the Woody Savannas in the northeast and southwest of the watershed, as well as on the Evergreen Broadleaf Forest located in the east and southwest, and the grasslands of the downstream part.

Croplands were generally less affected thanks to irrigation. Some rare positive anomalies were also found around the lake. For instance, the extreme drought of 2015 modified vegetation activity during the first semester of 2016. The temporal evolution of EVI of the Evergreen Broadleaf Forest, Woody Savannas and Cropland/Natural Vegetation Mosaic (Fig. 10) as well as Savannas, Grasslands and Croplands (Fig. S10) exhibit a longer and lower than usual period of low activity (until June instead of April) where it presents minimum activity between January and June instead of a maximum for Mixed Forests and Wetlands (Fig. S10).

Normal and even slightly above normal EVI values are found for the rest of 2016 and also for 2017 (Figs. 8b, 8c, S9a and S9b).

Conclusion

Surface

8.

Figures:

Figure 1

  Gallery forest (annually flooded for 9 months and covering an area of 197 km²)  Seasonally flooded habitats, composed of shrublands and grasslands (annually flooded from 5 to 8 months on average and covering an area of 5,409 km²)  Transitional habitats dominated by abandoned agricultural fields, receding rice/floating rice, and lowland grasslands (annually flooded from 1 to 5 months on average and covering an area of 3,658 km²)  Rainfed habitats mostly composed of wet season rice fields and village crops (flooded less than 1 month per year and covering an area of 8,641 km²)

  3.3.5. Climate indices 3.3.5.1 El Niño 3.4 Sea Surface Temperature index El Niño Southern Oscillation (ENSO) oceanic indices are defined using Sea Surface Temperature (SST) anomalies in different regions of the Pacific Ocean. El Niño 3.4 SST index is computed over a region located in the eastern Pacific, between longitude 120° W and 170° W and latitude 5° N and 5° S ("Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century," 2003). It is available on a monthly time step starting from January 1870 at: https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/Nino34/. 3.3.5.2 Pacific Decadal Oscillation index The Pacific Decadal Oscillation (PDO) index is defined as the principal component of monthly SST variability in the North Pacific, poleward of 20° N ("A Pacific Interdecadal Climate Oscillation with Impacts on Salmon Production," 1997). It is available on a monthly time step starting from January 1900 at: https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/PDO/.

  S5 in supplementary information). The annual lack of data due to cloud presence over Tonle Sap Lake and its floodplains is generally less than 8 days (1 image) and can reach up to 24 days (3 images), with standard deviations less than 16 days (2 images). This lack of data (not occurring on consecutive dates) does not prevent the use of MODIS images for monitoring Tonle Sap seasonal flooding. Maps of annual average duration of inundation extent and associated standard deviation are presented in Figures2a and 2b, respectively, along with maximum and minimum inundation extents during the flood peak in 2011(Figs. 2c) and 2015 (Figs. 2d), respectively, considering only pixels classified as inundated. Fig.3is similar to Fig.2, but uses both pixels classified as inundated and the ones classified as mixed pixel. Extensive inundation is generally in the western part of the lake during time spans shorter than 3

Figure

  Figure 4a presents temporal variations of surface water extent in the Tonle Sap drainage

  watershed and over the entire Mekong Basin from 1998-2017 (Fig.5b). Correlation coefficients between the interannual variations of surface water volume and SPAI (defined during the hydrological year between November of year n-1 and October of year n) are equal to 0.61 (p-value=5.1.10 -3 ) and 0.77 (p-value=10 -4 ) for the Tonle Sap watershed and the Mekong Basin, respectively. These results are related to the hydrological regime of the Tonle Sap watershed, whose flood pulse is driven by the flow of the Mekong River. Annual precipitation anomalies can account for the large droughts of 1998, 2010 and 2015 and large floods of 2000, 2001 and 2011 with SPAI lower than -1 and higher than 1, respectively, for both the Tonle Sap watershed and the Mekong Basin, but cannot explain the large flood of 2002 and large drought of 2016 according to TRMM/TMPA 3B43 dataset. Similar comparisons were performed during the dry and rainy seasons. Correlation coefficients were found equal to 0.31 (p-value=0.19) and 0.84 (p-value close to 0) during the dry and rainy seasons, respectively, for the Mekong Basin and to 0.45 (p-value=4.5.10 -2 ) and 0.58 (p-value=7.5.10 -3 ) during the dry and rainy seasons, respectively, for the Tonle Sap watershed. These results confirm that surface water volume of the Tonle Sap watershed is strongly influenced by rainfall during the rainy seasons in the Mekong Basin. During the large flood events of 2000, 2001 and 2011, SPAI during the rainy season of between 1 and 1.5 were obtained in the whole Mekong Basin (Fig. 5c). Larger floods were recorded in 2000 and 2011 as they were enhanced by SPAI between 1 and 1.5 in the Tonle Sap watershed. A very low positive SPAI was observed in the Tonle Sap watershed during the wet season of 2001.

  differences in the impact of extreme events on lake volume. The drought of 1998 and the floods of 2000 and 2011 are directly related to large positive and negative anomalies of both ENSO 3.4 and PDO. The 2002 flood and 2015/2016 droughts can be attributed to both a large value in module of either ENSO 3.4 or PDO (and a lower value of the other index), and the existence of an excess (since 2000) or a deficit (since 2014) of water in the lake. This shows that the possibility of occurrence of a hydrological extreme is not only conditioned by climate factors but also by base water storage from previous years.

  water storage variations of Tonle Sap were determined combining remotely sensed inundation extent data derived from MODIS reflectance images and altimetry-based water levels from 1993-2017. Their interannual anomalies were found to be more related to the rainfall occurring in the Mekong Basin (R=0.77 for the interannual variations during the whole hydrological year) than in the Tonle Sap watershed (R=0.61) as flooding in the Tonle Sap watershed is mostly driven by the flow of Mekong River. Nevertheless, extreme events (exceptional floods and droughts) occurred when large excess/deficit rainfall occurred during the rainy season in both the Mekong Basin (R=0.84) and Tonle Sap watershed (R=0.58). Excess rainfall occurring in both dry and rainy seasons in the Mekong Basin combined with a rainfall deficit in the Tonle Sap watershed during the rainy season was only seen to cause large flooding in 2001. Important memory effects were also detected as in the extreme flood of 2002. This was mostly due to excess in surface water storage of the exceptional flood of 2000 and the large flood in 2001, whereas the extreme drought of 2016 was caused by the deficit in surface water by the drought of 2014 and the extreme drought of 2015. A quasiequal influence of the West North Pacific and Indian Monsoons was also found with correlations between the interannual variations of surface water storage and the WNPMI in July of 0.55 and the IMI in August of 0.49. Interannual variations of surface water storage in the Tonle Sap watershed are strongly influenced by climate variability. They exhibit large positive and negative anomalies corresponding to extreme drought (1998, 2015, and 2016) and flood (2000, 2001, 2002, and 2011) events. They were related to the very strong El Niño events of 1997-1998 and 2015-2016 and the moderate La Niña events of 1999-2000 and 2010-2011 enhanced by the PDO. Correlation coefficients between the interannual variations of surface water storage of the Tonle Sap watershed and the average ENSO 3.4 index over DJF and PDO over NDJFM were found equal to -0.70 and -0.68, respectively. R of -0.75 was found using a linear combination of the two former indices. Below normal EVI values for the types of vegetation were observed in the first semester of 2016 as a consequence of the 2015 drought.The NASA/CNES Surface Water and Ocean Topography (SWOT) mission is to be launched in 2021 and will improve flood monitoring in Tonle Sap by providing water level maps through SAR interferometry techniques at better spatial and temporal resolutions with low impact of cloud presence and vegetation. The data will provide a better understanding of the strong relationship between climatic indices and Tonle Sap lake volume extreme changes, suggested by[START_REF] Räsänen | Spatiotemporal influences of ENSO on precipitation and flood pulse in the Mekong River Basin[END_REF] for the whole Mekong basin and therefore help in forecasting the risk of extreme droughts or floods.7.(VolTransMesKong)" project and to the "Centre de Topographie de l'Océan et de l'Hydrosphère (CTOH)". It is based on observations with nadir radar altimeters on board Topex/Poseidon, Jason-1, Jason-2, Jason-3, ERS-2, ENVISAT, and SARAL satellites and with MODIS instruments on board Terra and Aqua satellites. C. Normandin is supported by a PhD grant from Ministère de l'Enseignement Supérieur et de la Recherche. The Mekong River Commission (MRC) is also warmly thanked for providing on its website (www.mrcmekong.org) many resources concerning the Mekong basin. We thank Dr. S. Calmant, Dr. J.-F. Crétaux and two anonymous reviewers for their very useful comments that helped improve this manuscript.

  Figure 1: a) The Tonle Sap watershed (cyan) is located in the Mekong River Basin (black), the Mekong River and its major tributaries appear in blue; b) 10-day (Jason-1, 2 and 3 in red) and 35-day (ERS-2, ENVISAT and SARAL in green) repeat period tracks over the Tonle Sap Lake watershed. Mekong River Basin boundary and drainage networkcome from the Major River Basins of the World from the Global Runoff Data Centre -Federal Institute of Hydrology (BfG) available at http://www.bafg.de/GRDC/EN/02_srvcs/22_gslrs/221_MRB/riverbasins_node.html. Tonle Sap boundaries were provided by Open Development Cambodia (http://www.opendevelopmentcambodia.net/maps/downloads/). Altimetry tracks were provided by AVISO+ (https://www.aviso.altimetry.fr/en/data/tools/pass-locator.html). Background images come from Google Earth.

Figure 2 :

 2 Figure 2: a) Maps of annual flood duration averaged over 2000/2017 time period and b) associated standard deviations expressed in days, using only pixels classified as "inundated". c) Maps of flood duration in 2011 (year with the maximum flood extent over the whole time period) and d) in 2015 (year with the minimum flood extent), considering only pixels classified as inundated. Background images come from Google Earth.

Figure 3 :

 3 Figure 3: a) Maps of annual flood duration averaged over 2000/2017 time period and b) associated standard deviations expressed in days, using pixels classified as "inundated" and pixels classified as "mixed pixels". c) Maps of flood duration in 2011 (year with the maximum flood extent) and d) in 2015 (year with the minimum flood extent), considering inundated and mixed pixels. Background images come from Google Earth.

Figure 4 :

 4 Figure 4: Time series of a) surface water extent (in km 2 ) in the Tonle Sap basin based on the processing of MODIS images, over 2000-2017, b) altimetry-based water level of the Tonle Sap Lake (in m wrt EGM08 geoid), combining information from Topex/Poseidon(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002),ERS-2 (1995ERS-2 ( - 2003)),Jason-1 (2002Jason-1 ( -2008)),ENVISAT (2003[START_REF] Stone | Severe Drought Puts Spotlight on Chinese Dams[END_REF]), Jason-2 (2008-2016), SARAL (2013-2016) and Jason-3 (since 2016), c) water volume variations of the Tonle Sap lake obtained applying the surfaceheight relationship over 1993-2017.

Figure 5 :

 5 Figure 5: a) Interannual variations of annual surface water storage. b) Annual, c) dry and d) wet seasons SPAI in the Mekong Basin and in the Tonle Sap watershed.

Figure 6 :

 6 Figure 6: Interannual variations of monthly a) surface water volume of the Tonle Sap watershed, b) anomalies of ENSO 3.4 and c) PDO indices over 1993-2017.

Figure 7 :

 7 Figure 7: Interannual variations of annual a) surface water storage between 1993 and 2017, b)

Figure 10 :

 10 Figure 10: Monthly variations of EVI averaged over 2000-2017 (continuous line) and EVI for 2016 (dotted line) for the Evergreen Broadleaf Forest (green), Woody Savannas (orange) and Cropland/Natural Vegetation Mosaic (black).

  The pixels are classified into two major classes: non-

	4.2.	Altimetry-based water levels
		4.2.1 Post-processing of altimeter data

flooded (EVI > 0.3 or EVI ≤ 0.3 but EVI -LSWI > 0.05) and water-related pixels (EVI ≤ 0.3 and EVI -LSWI ≤ 0.05 or EVI ≤ 0.05 and LSWI ≤ 0). This latter class is divided into 3 subclasses: mixture pixels when 0.1< EVI ≤ 0.3, flooded pixels EVI ≤ 0.1, and long-term water bodies (e.g., rivers and lakes) when a pixel is determined flooded more than 250 days per year. Spatio-temporal variations of the flood were determined in the Tonle Sap drainage basin from February 2000 to December 2017.

  for the Tonle Sap Basin. Time series of Enhanced Vegetation Index (EVI) and their interannual variations (estimated by removing the climatological monthly mean to the monthly values) were computed for each land cover type. In the following, only land cover types corresponding to vegetation are considered. They are presented for the three main land cover types which cover 80.3% o of the Tonle Sap basin surface (Evergreen Broadleaf Forest, Woody Savannas and Cropland/Natural Vegetation Mosaic) as seen in Figures8b and cfor the time series and associated anomalies, respectively, and in FiguresS9a and bfor the five other classes of vegetation which cover 14.6% of the surface of the Tonle Sap Basin surface

	(Mixed Forests, Savannas, Grasslands, Permanent Wetlands and Croplands). Lower than
	usual minimum values of EVI are observed in spring 2016 for all land cover types as shown
	in Figures 8b and S9a. The anomalies reached their lowest values in April 2016 for Evergreen
	Broadleaf Forest, Woody Savannas and Cropland/Natural Vegetation Mosaic with EVI

November to March (NDJF).

Acknowledgements

This work was supported by the CNES, through TOSCA grants attributed to the "Variabilité hydro-sédimentaire du bassin hydrologique du MEKONG par télédétection