Sabine Szunerits

Rabah Boukherroub

GRAPHENE BASED BIOELECTROCHEMISTRY AND BIOELECTRONICS: A CONCEPT FOR THE FUTURE?

Introduction

Clinical analysis benefits world-wide from a variety of diagnostic tests. These analytical tests should be fast and highly accurate to help in establishing a treatment protocol that is appropriate for the patient. The interest in the development of new clinical tests is not only driven by the demand to sense new analytes, but also to reduce costs, complexity and lengthy analysis times of current techniques. Among the myriad of possibilities available today, electrochemical and field-effect based biosensors are prominent players [START_REF] Bunyakul | Combining Electrochemical Sensors with Miniaturized Sample Preparation for Rapid Detection in Clinical Samples[END_REF][START_REF] Kaisti | Detection principles of biological and chemical FET sensors[END_REF]. These analytical platforms require little and cheap instrumentation, can provide low detection limits once optimized, and be easily miniaturized, offering many of the desirable attributes for pointof-care tests. The best known example of an electrochemical based diagnostic device is the glucose biosensor, used widely in commercial glucometers and pH electrodes. Due to interest of patients as well as clinicians alike in these efficient, fast and simple to use sensing platforms, the number of biosensors developed in research laboratories integrating the market are gradually increasing.

A typical bioelectronics sensor contains two basic functional units: (i) an electrical transducer and (ii) a recognition element. The sensor selectivity is governed by surface-linked recognition elements (oligonucleotides, antibodies, enzymes, etc) as well as anti-fouling molecules (polyethylene glycol derivatives, serum molecules, etc) to limit non-specific signal, while the sensitivity is conferred by the transducer, converting the recognition event into a measureable electrical signal. One of the challenges in this field, concerns the sensitive and selective recording of biological recognition events in complex media such as human serum and salvia where picomolar (pM) and even lower detection limits for biological analytes are often requested. For advanced biosensors, the choice of the electrical transducer, the To whom correspondence should be send to: sabine.szunerits@univ-lille1.fr recognition element and the surface linking strategy, have to be considered simultaneously, making the construction of viable biosensing platforms of at great challenge. The attractive properties of graphene based nanomaterials have paved the way for the fabrication of a wide range of electrical and electrochemical based biosensors with improved analytical performance [START_REF] Ibau | Current advances and future visions on bioelectronic immunosensing for prostate-specific antigen[END_REF][START_REF] Pandikumar | Graphene and its nanocomposite material based electrochemical sensor platform for dopamine[END_REF][START_REF] Ambrosi | Gold-graphene nanocomposites for sensing and biomedical applications[END_REF][6]. Different architectures are proposed to integrate graphene-based materials with a range of detection techniques implemented to record the biological recognition events.

The materials science aspect has become one of the main research focus in the last years for graphene based bioelectronics platforms and the examples of the type of graphene materials used are almost countless (Figure 1).

This short review provides insights into the various benefits of graphene based bioelectrical and bioelectrochemical sensors. Situated at the interface of materials science, chemistry and the life sciences, graphene based electronics and electrochemistry offer a broad palette of opportunities for researchers and clinicians for targeted theranostics in biology and medicine.

Ligand modified graphene

Electrochemical based sensing: Aspects about electron transfer on different graphene modified interfaces

The underlying principle of an electrochemical sensor is to convert a biological binding event into a measurable electrical signal. This approach is represented by electron transfer between the transducer, being graphene coated electrical interfaces in our case, and an electroactive species, which can be either the molecule to be analyzed or a species, which electrochemical signal correlates with the presence of the target analyte. Aspects governing electron transfer on graphene coated electrodes become of importance to obtain large output signals, which will allow reaching the low detection limits needed when analyzing clinical samples. Given the enormous variety of graphene-based materials (Figure 1), the question of where electron transfer occurs and what are the parameters influencing electron transfer rate are often difficult to answer, and mostly neglected in the construction of bio-electrochemical sensing interfaces. Impressive electrocatalytic performance was obtained on multi-layered graphene flakes coated Si wafers in response to biomolecules such as dopamine, ascorbic acid, attributed to the large area of edge planes that are available on such an interface to allow rapid heterogeneous electron transfer [START_REF] Shang | Catalyst-Free Efficient Growth, Orientation and Biosensing Properties of Multilayer Graphene Nanoflake Films with Sharp Edge Planes †[END_REF]. Simple solution processable reduced graphene oxide (rGO), which enables large scale production, results in a graphene material with structural defects due to the presence of different oxygen functional groups, which are never completely removed. While these properties might impact the electrical readout, these defects are actually beneficial to achieve fast heterogeneous electron transfer and promote the possibility for covalent immobilization of surface ligands to target specific sensing applications [START_REF] Ambrosi | Gold-graphene nanocomposites for sensing and biomedical applications[END_REF][START_REF] Zhou | Electrochemical Sensing and Biosensing Platform Based on Chemically Reduced Graphene Oxide[END_REF]. Kamat and co-workers have lately looked into the aspect of controlled electron capture, transport and discharge of rGO modified glassy carbon electrodes (GCE), showing that cooperative effects of two redox couples (in this case ferrocene and methyl viologen) on a rGO platform can increase the sensitivity of electrochemical detection, something still not exploited for sensing, but opening up the possibility of detection at low levels [START_REF] Bridewell | Electrocatalytic Sensing with Reduced Graphene Oxide: Electron Shuttling between Redox Couples Anchored on a 2-D Surface[END_REF].

The question is ultimately whether single-layer two-dimensional graphene, reduced graphene oxide (rGO), porous rGO, doped and functionalised rGO derivatives exhibit more rapid electron transfer over other electrodes and highly ordered pyrolytic graphite (HOPG).

Comparison is a very difficult task as the way the graphene materials are coated onto the electrical interface will influence strongly the electrochemical behaviour. We opted on highlighting the complex issue by recording cyclic voltammetric profiles in ferrocyanide/KCl solutions of a gold electrode modified with CVD graphene [10], and by drop casting different suspensions of rGO (1 mg mL -1 in water) (Figure 2A) [START_REF] Chekin | Reduced Graphene Oxide Modified Electrodes for Sensitive Sensing of Gliadin in Food Samples[END_REF][START_REF] Vasilescu | Porous reduced graphene oxide modified electrodes for the analysis of protein aggregation. Part 1: Lysozyme aggregation at pH 2 and 7.4[END_REF][START_REF] He | Copper oxide supported on three-dimensional ammonia-doped porous reduced graphene oxide prepared through electrophoretic deposition for non-enzymatic glucose sensing[END_REF][14]. Deposition of CVD graphene results in a slight increase in the voltammetry peak height; the largest increase was achieved on N-doped porous reduced graphene oxide, with the peak-to-peak potential separation remaining unchanged. The increase in the electron transfer rate on chemically transferred CVD graphene nanosheets is mostly due to the enhancement of the surface area as in fact a single layer is not completely flat and has a tendency to fold, buckle and corrugate. The presence of oxygen-related edge and plane functions as is the case in rGO based materials together with the increased active surface area results in superior electron transfer rates, surpassing those of gold and gold/CVD graphene [START_REF] Lim | Direct Voltammetric Detection of DNA and pH Sensing on Epitaxial Graphene: An Insight into the Role of Oxygenated Defects[END_REF]. The intrinsic chemical doping with heteroatoms is an effective way for modulating and tailoring the physico-chemical properties of graphene. Substitution with electron-rich atoms, such as nitrogen, which is able to form strong valence bonds with carbon atoms, can lead to an increase in free charge carriers in the graphene's aromatic ring network and therefore an enhanced conductivity. As seen in Figure 2A, porous nitrogen-doped graphene nanostructures lead to a further increase in the electrochemical electroactivity; this underlines the general trend that the sample with the highest defect density displays the fastest electron transfer kinetics [START_REF] Zhou | Electrochemical Sensing and Biosensing Platform Based on Chemically Reduced Graphene Oxide[END_REF][START_REF] Lim | Direct Voltammetric Detection of DNA and pH Sensing on Epitaxial Graphene: An Insight into the Role of Oxygenated Defects[END_REF]. This highlights that electron transfer rates depend strongly on the type, morphology, as well as on the film thickness of the graphene coatings, where each parameter acts differently on the electrochemical activity of graphene. The materials science part of graphene and its derivatives is thus of fundamental importance to overcome the ambiguity about the electrochemical properties of graphene. To further demonstrate how the electron transfer kinetics can affect the sensitivity and linear range of an electrochemical sensor, a toponin I (TnI) specific sensor was constructed lately upon covalent integration of aptamer Tro4 (Figure 2B). Addition of TnI results in a decrease of the differential pulse voltammetry (DVP) signal of ferrocyanide, an inner-sphere redox couple sensitive to surface functions (Figure 2C), which scales linearly with the TnI concentration down to 1 pg mL -1 over a wide linear range of 5 orders of magnitude between 1 pg mL -1 and 100 ng mL -1 with a slope of 43 µA cm -2 /decade (r 2 =0.9992) (Figure 2D).

Graphene based biological field-effect transistors (GFET)

Utilizing a 2D channel material such as graphene in back or liquid gated bioFETs has several advantages over classical FET technology. For most semiconductor-based transistor sensors, local electric field changes at the channel surface have little effect deeper in the device channel, limiting the response sensitivity. With a GFET (Figure 3A), the graphene channel is only one atom thick, meaning the entire channel is effectively on the surface and directly exposed to the environment. Any molecule attached to the surface of the channel impacts electronic transfer through the entire depth of the device and can be sensed. When target molecules bind to the receptors on the graphene surface, the redistribution of electronic charge generates a change in the electric field across the GFET channel region, which changes the electronic conductivity in the channel and the overall device response as seen in the transfer curve of a liquid-gated GFET for the analysis of DNA hybridization events (Figure 3B). Most importantly, for sensing applications using GFETs, the Dirac point is sensitive to the presence of immobilized DNA probes as well as to the hybridization of complementary target DNA [START_REF] Dong | Electrical detection of DNA hybridization with single-base specificity using transistors based on CVD-grown graphene sheets[END_REF].

The main concern with biological GFETs remains the effect of ionic, also called double-layer screening (Figure 3C). This screening effect is dependent on the distance between the FET sensing surface and the point of observation and is characterized by the Debye screening length ( D ), defined by the distance where the electrical signal decays to 1/e of its original value. Unless very diluted solutions or salt free solutions are considered, typical screening lengths are in the order of 1 nm with stronger screening effects for larger molecules. The use of small molecular sized ligands (antigen binding fragments, aptamers) are thus currently reviling classical FET based immunosensors (Figure 3D) [START_REF] Cheng | Field Effect Transistor Biosensor Using Antigen Binding Fragment for Detecting Tumor Marker in Human Serum[END_REF]. (Reproduced with permission from Ref. [START_REF] Dong | Electrical detection of DNA hybridization with single-base specificity using transistors based on CVD-grown graphene sheets[END_REF]; (C) Electrical double-layer length in the presence of different targets (dimensions are not scaled) (Reprinted with permission from Ref. [START_REF] Huang | Electrochemical processes and mechanistic aspects of field-effect sensors for biomolecules[END_REF]), (D) Influence of ligand size on FET based sensing.

OPINION ON THE FUTURE OF GRAPHENE BASED ELECTRICAL AND ELECTRIOCHEMICAL SENSING

The examples of graphene based electrochemical and electrical based biosensors are countless and vary not only in the way the graphene based interface is used, but also in the manner electrical or electrochemical read out is achieved (Table 1). While some efforts have been made here by the integration of polyethylene modified pyrene ligands as well as blocking with human serum itself, more efforts have to be put into this direction to optimize the selectivity of the sensor when used in complex biological media. The other limiting factor for commercialisation is the current limitation of large-scale reproducibility of the fabrication of graphene biosensor interfaces.

In the case of bioelectronics, the overbearing sensitivity-limiting factor from the Debye screening length must be mitigated to allow for sensing in high ionic strength solutions such as serum, saliva, etc. A detailed understanding of aptamer/analyte affinity interaction and strength becomes a fundamental issue here to ensure promising bioelectronics detection of various proteins. This is also important to convince industry to change from the well-known and highly specific antibody strategies to less selective aptamer technology.

CONCLUSION AND OUTLOOK

Graphene and its derivatives prove to be a rich source for the construction and use of bioelectronics and bioelectrochemical sensors. Routinely, such sensors achieve a picomolar detection limit, with some even reaching the low femtomolar concentration range. The possibility that a large range of different detection methods can be employed with graphene based sensors is of high advantage, as depending on the looked after final application, sensor size and read out can be customised at will. The collaboration between materials scientists, chemists, physicists as well as engineers and medical personal is of fundamental importance to drive this field further and to propose graphene based biosensors as point-of-care alternatives for patients. The success of any new biosensor material lies in addition on its reproducibility and possible industrial scale production. The emergence of a number of companies providing mono-and bilayered graphene nanosheets on several interfaces, graphene oxide (GO), reduced graphene oxide (rGO) and even modified matrixes has been an additional motivation for using graphene for biosensor applications. When it comes to in vivo application of some of the sensing concepts, graphene based biosensors are still in their infancy. Toxicity and biocompatibility issues still need to be addressed carefully to avoid any undesired secondary health effects. Current in vivo and in vitro assessments of the biostability of the sensors are encouraging and promising for further technological transfer. Finally, beside various integration techniques insuring highly sensitive sensing, the prevalence of the current mobile area has introduced the emergence of portable and smartphone based analytical quantification read outs using wireless technology. It will be exciting to witness the futures of these wireless-based sensing devices which will finally make real personalized diagnostic possible.

Figure 1 :

 1 Figure 1: Different components for the construction and use of graphene based bioelectronics and bioelectrochemical sensors.

Figure 2 :

 2 Figure 2: (A) Current vs. potential plots recorded by cyclic voltammetry on Au (black), Au/graphene (chemically transferred of CVD grown graphene); Au/rGO (red), Au/prGO (blue) and Au/N-doped prGO using [Fe(CN) 6 ] 4-(10 mM)/PBS (0.1 M), scan rate = 100 mV s - 1 ; (B) Schematic presentation of the construction of an immunosensor for TnI; (C) Change in differential pulse voltammogram signal upon addition of TnI on N-doped pRGO modified electrode; (D) Calibration curve for TnI on different interfaces.

Figure 3 :

 3 Figure 3: (A) Image of a liquid-gated GFET together with change in drain-source current upon a binding event; (B) Transfer characteristics of a GFET before adding DNA, after immobilization of probe DNA, and after reaction with (b) complementary or (c) one-base mismatched DNA molecules with the concentration ranging from 0.01 to 500 nM.(Reproduced with permission from Ref.[START_REF] Dong | Electrical detection of DNA hybridization with single-base specificity using transistors based on CVD-grown graphene sheets[END_REF]; (C) Electrical double-layer length in the presence of different targets (dimensions are not scaled) (Reprinted with permission from Ref.[START_REF] Huang | Electrochemical processes and mechanistic aspects of field-effect sensors for biomolecules[END_REF]), (D) Influence of ligand size on FET based sensing.

Table 1 :

 1 Selected recent examples of bioelectronic and bioelectrochemical graphene based sensors. CEA: carcinoembryonic antigen; DPV: differential pulse voltammetry; ESI: Electrochemical impedance spectroscopy; EPD: electrophoretic deposition; FAB: folic acid protein; G: pristine graphene, PA: 1-pyrenebutyric acid, PASE: 1-pyrenebutanoic acid succinimidyl ester, prGO: porous reduced graphene oxide, PSA: prostate specific-antigen; rGO: reduced graphene oxide; SNP: single nucleotide polymorphism; SWP: single sweep potential.While advances in graphene materials have shown to be highly crucial for improving electrical and electrochemical biosensor characteristics, different challenges are still to overcome. There is still an urgent need for moving beyond research by developing new concepts for achieving even better sensitivity and selectivity, in order to bring some of the current sensors into real biomedical applications. Even though a large amount of graphene based sensors reported in the literature exhibit good sensing characteristics, the performance in real biological samples, where next to a variety of other proteins high salt concentrations are present, is often not considered and information is largely lacking. As non-specific interactions of these proteins are of primary concern on graphene based interfaces, this lack of information is one of the two most crucial limiting factors for current commercialization.

	analyte	electrode	method	Linear range	LoD	comments		Ref
	CEA	G	FET		0.5 pM	PASE+anti-CEA	[19]
	micro-RNA	rGO	FET		10 fM	Au-NPs modified	[20]
						PNA
	SNP	G	FET	100 nM-100 µM	100 nM	+1-pyCCOH		[21]
	nucleic acid	Thiofluoro-	EIS			nucleic	acid	[22]
		graphene				adsorption
	Interleukin-8 rGO	DVP	500 fg mL -1 -4 ng	72.7 pg mL -1 Au NPs+anti-IL8	[23]
				mL -1
	PSA	rGO	DPV	90 ng mL -1 -0.1 pg	10 fg mL -1	+aptamer		[24]
				mL -1
	lysozyme	prGO	DVP	0.05-7.5 µM	50 nM	Direct oxidation	[12]
	FAP	rGO (EPD)	DVP	1-100 pM	1 pM	folic acid modified [13]
	gliadin	prGO	DVP	1.2-34 ng ml -1	1.2 ng ml -1	PA+anti-gliadin	[11]
	PSA	NH 2 -graphene EIS	2 pg ml -1 -2 µg ml -1	0.46 pg ml -1	+ anti-PSA and Au	[25]
						NP amplification
	thrombin	Laser-scribed	DVP	1 pM-0.1 nM	1 pM	PA+thrombin		[26]
		G-electrode				aptamer
	VEGF165	rGO	SWP	10 fg ml -1 -1 ng ml -1 8 fg mL -1	Cucurbituril and	[27]
						N 3 modified rGO
						+alkynyl-
						DNA+aptamer
	Dopamine	rGO/PEI	DVP		50 nM	microsystem		[28]
		(EPD)
	Folic acid	rGO/MoS 2	DVP	0.01-100 µM	10 nM	Direct oxidation	[29]
	glucose	rGO/N-doped	CA	0.25µM-6 mM	0.25 µM	EPD		[14]
		prGO/CuO
	H 2 O 2	rGO/Pt NPs	CA	0.2 nM-3 µM	20 nM	chitosan-ferrocene	[30]
						carboxylic acid
	H 2 O 2	rGO-TiO 2	CA	0.1-360 µM	10 nM	+hemoglobin		[31]

Acknowledgements

This work has emanated in part from research conducted with the financial support from the Centre National de la Recherche Scientifique (CNRS), the University of Lille, the Hauts-de-France region, the CPER "Photonics for Society", the Agence Nationale de la Recherche (ANR) through FLAG-ERA JTC 2015-Graphtivity and the EU through the Marie Sklodowska-Curie action (H2020-MSCA-RISE-2015, PANG-690836).


G L Ramillien

F Frappart

S Gratton

X Vasseur

Sequential estimation of surface water mass changes from daily satellite gravimetry data

Keywords: GRACE satellite gravimetry, Continental hydrology, Kalman filtering, Regional solutions

We propose a recursive Kalman filtering approach to map regional spatio-temporal variations of terrestrial water mass over large continental areas, such as South America. Instead of correcting hydrology model outputs by the GRACE observations using a Kalman filter estimation strategy, regional 2-by-2 degree water mass solutions are constructed by integration of daily potential differences deduced from GRACE K-band range rate (KBRR) measurements. Recovery of regional water mass anomaly averages obtained by accumulation of information of daily noise-free simulated GRACE data shows that convergence is relatively fast and yields accurate solutions. In the case of cumulating real GRACE KBRR data contaminated by observational noise, the sequential method of step-by-step integration provides estimates of water mass variation for the period 2004-2011 by considering a set of suitable a priori error uncertainty parameters to stabilize the inversion. Spatial and temporal averages of the Kalman filter solutions over river basin surfaces are consistent with the ones computed using global

Introduction

Launched in March 2002, the Gravity Recovery and Climate Experiment (GRACE) mission has globally mapped the temporal variations of the Earth's gravity field at an unprecedented millimetre precision in terms of geoid height thanks to the accurate inter-satellite measurements made by the onboard K-band range (KBR) system [START_REF] Tapley | The gravity recovery and climate experiment: mission overview and early results[END_REF], and by now for more than 12 years. Pre-processing of Level-1 GRACE data consists of removing the effects of known a priori gravitational accelerations such as static gravity field, atmosphere and ocean mass changes, pole and oceanic tides from the measurements to produce Level-2 solutions that consist of "reduced" Stokes coefficients (i.e., dimensionless spherical harmonic coefficients of the geopotential-see theoretical aspects in Hofmann-Wellenhof and Moritz 2006) up to degree 90, or equivalently, of spatial resolutions of 200-300 km [START_REF] Bettadpur | CSR Level-2 processing standards document for level-2 product release 04, GRACE. The GRACE project[END_REF][START_REF] Fletchner | GFZ Level-2 processing standards document for level-2 product release 04[END_REF][START_REF] Chambers | Evaluation of Release 05 time-variable gravity coefficients over the ocean[END_REF][START_REF] Dahle | GFZ RL05: an improved time series of monthly GRACE gravity field solutions[END_REF]. These global Level-2 solutions correspond to the gravity signatures of not modelled surface mass variations such as sudden earthquakes and continuous glacier melting, but mainly water mass transport for detecting more localized hydrological events lasting a few days.

Previously, [START_REF] Kurtenbach | Deriving daily snapshots of the Earth's gravity field from GRACE L1B data using Kalman filtering[END_REF][START_REF] Kurtenbach | Improved daily GRACE gravity field solutions using a Kalman smoother[END_REF] and [START_REF] Sabaka | Improving global mass flux solutions from Gravity Recovery and Climate Experiment (GRACE) through forward modeling and continuous time correlation[END_REF] have proposed a recursive Kalman filter scheme for estimating daily Level-1B GRACE data taking into account statistical information on process dynamics and noise from geophysical models to gain in temporal resolution. [START_REF] Kurtenbach | Deriving daily snapshots of the Earth's gravity field from GRACE L1B data using Kalman filtering[END_REF] have applied a Kalman filtering after having removed annual and semi-annual parts to obtain gravity variation time series of Stokes coefficients at daily intervals. These corrections have enabled them to consider a stationary isotropic process noise in time and space to derive an exponential covariance function of a first-order Markov process. In describing this process noise, [START_REF] Kurtenbach | Deriving daily snapshots of the Earth's gravity field from GRACE L1B data using Kalman filtering[END_REF] have used a priori time and space covariances on Stokes coefficients. Unfortunately, daily spherical harmonics solutions remain smooth in space, due to the poor geographical coverage of the daily satellite tracks.

In this article, no a priori model is corrected by the GRACE observations, but the estimate is built iteratively by cumulating information of daily GRACE satellite tracks. As there are no daily GRACE satellite data to cover the surface of the Earth efficiently, we demonstrate that estimating regional time-cumulated solutions by the successive injection of daily GRACE tracks is achievable. The progressive integration of satellite observations has the advantage of avoiding inversion of large systems of equations in the determination of time-constant maps of water mass variation. The sequential method we propose is inspired from previous methods. However, our Kalman filter estimation is applied to regional solutions instead of global Stokes coefficients. First, the methodology of the two stages of the Kalman filter estimation (i.e., projection and correction) is presented, and secondly, its application for inverting simulated GRACE satellite data is shown. Thirdly, both the analysis of the convergence and the detection of short-term events are made by noise-free simulations from hydrology model outputs. This method is applied to invert along-track potential anomalies derived from real GRACE KBRR measurements to recover multi-year series of 2-by-2 degree solutions of water mass variations over South America. In the following, we use the term "regional averaged (or equivalently, cumulated) solutions computed at successive daily time steps", but this does not imply a daily time resolution. Then, these estimates are confronted to the official global monthly GRACE solutions provided by CSR, GFZ and JPL (Release 5) filtered using the independent component analysis approach [START_REF] Frappart | An independent Component Analysis filtering approach for estimating continental hydrology in the GRACE gravity data[END_REF], and 10-day regional solutions from [START_REF] Ramillien | Constrained regional recovery of continental water mass time-variations from GRACE-based geopotential anomalies over South America[END_REF], 10day global solutions from GRGS [START_REF] Bruinsma | CNES/GRGS 10-day gravity field models (release 2) and their evaluation[END_REF], and GLDAS-NOAH [START_REF] Rodell | The global land data assimilation system[END_REF]) and WGHM [START_REF] Hunger | Value of river discharge data for globalscale hydrological modeling[END_REF] outputs, as well as rapid events described by river discharge observations for validation. The possibility of detecting short-term water mass events is finally over continental areas. North-south striping is particularly visible in the tropics where coverage of satellite tracks is poor due to (1) sparse GRACE track sampling in the latitudinal direction, (2) propagation of errors of the a priori correcting model accelerations [START_REF] Han | Time-variable aliasing effects of ocean tides, atmosphere, and continental water mass on monthly mean GRACE gravity field[END_REF][START_REF] Thompson | Impact of short period, non-tidal, temporal mass variability on GRACE gravity estimates[END_REF][START_REF] Ray | Tide model errors and GRACE gravimetry: towards a more realistic assessment[END_REF]) and (3) numerical correlations generated while solving the undetermined systems of equations for high-degree Stokes coefficients [START_REF] Swenson | Post-processing removal of correlated errors in GRACE data[END_REF]).

Short-term mass variability with periods from hours to days of ocean tides and atmosphere is removed using dealiasing techniques of correcting model outputs. While the atmosphere pressure fields from ECMWF allow a reasonable de-aliasing of high frequency caused by non-tidal atmospheric mass changes, errors due to tide model appear in the GRACE solutions, especially from diurnal (S1) and semidiurnal (S2) tides [START_REF] Han | Time-variable aliasing effects of ocean tides, atmosphere, and continental water mass on monthly mean GRACE gravity field[END_REF][START_REF] Ray | Tide model errors and GRACE gravimetry: towards a more realistic assessment[END_REF][START_REF] Forootan | Comparisons of atmosphere mass variations derived from ECMWF reanalysis and operational fields, over 2003 to 2011[END_REF]. [START_REF] Thompson | Impact of short period, non-tidal, temporal mass variability on GRACE gravity estimates[END_REF] showed that the degree error increased by a factor ∼20 due to atmospheric aliasing, ∼10 due to ocean model and ∼3 due to continental hydrology, e.g., aliasing of the S2 tide has a strong impact on the determination of the C 20 spherical harmonic coefficient of the geopotential [START_REF] Seo | GRACE's spatial errors[END_REF]. See, e.g., [START_REF] Guo | Non-isotropic Gaussian smoothing and leakage reduction for determining mass changes over land and ocean using GRACE data[END_REF] to know about the reduction of atmosphere aliasing by Gaussian smoothing.

A regional approach alternative to the classical global one has been more recently proposed to improve geographical localization of patterns of water storage change. This energy integral method consists of recovering equivalentwater thicknesses of juxtaposed 2-degree tiles from GRACE inter-satellite velocity residuals, by considering matrix regularization techniques for solving ill-posed problems, but without using spherical harmonics for representing water mass variations (see [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF][START_REF] Ramillien | Constrained regional recovery of continental water mass time-variations from GRACE-based geopotential anomalies over South America[END_REF]. These multi-year series of 10-day regional solutions are comparable to independent datasets, such as local water level in the Amazon Basin, and thus they provide realistic amplitudes of water mass change from seasonal to inter-annual timescales over South America.

So far, global and regional solutions have been produced as weekly, 10-day and monthly averages from GRACE data, leading to a loss of resolution in time. In other words, global or regional averaging suppresses events with a gravity signature of a few days and thus does not enable us to capture them. Moreover, determining the regional equivalent-water heights for each separate time interval would neglect the temporal correlation between successive independent intervals imparted by the temporal dynamics of the water mass processes. While regional solutions computed as averages on constant 10-day intervals proved to be realistic snapshots of the surface water mass variability (Frappart et al. 2013a, b;[START_REF] Seoane | Regional GRACE-based estimates of water mass variations over Australia: validation and interpretation[END_REF], the next challenge is to attempt to improve the resolution of these regional solutions in time discussed and tested by simulating localized water mass anomalies.

Datasets

Independent hydrological datasets

WGHM land water storage

The WaterGAP Global Hydrology model (WGHM) [START_REF] Döll | A global hydrological model for deriving water availability indicators: model tuning and validation[END_REF][START_REF] Hunger | Value of river discharge data for globalscale hydrological modeling[END_REF]) is a conceptual model that simulates the water balance on continental areas at a spatial resolution of 0.5 • . It describes the continental water cycle using several water storage compartments that include interception, soil water, snow, groundwater and surface water (rivers, lakes and wetlands). We consider the sum of all these contributions and called it Total Water Storage (TWS) to be comparable to GRACE observations, as these latter data correspond to the integrated continental hydrology change without distinguishing each compartment. WGHM has been widely used to analyse spatio-temporal variations of water storage globally and for large river basins [START_REF] Günter | A global analysis of temporal and spatial variations in continental water storage[END_REF]. In this study, we use daily TWS grids from the model version WGHM 2.1f described by [START_REF] Hunger | Value of river discharge data for globalscale hydrological modeling[END_REF] to simulate the GRACE hydrology-related geopotential anomalies, and recover the starting TWS variation grids, as accurately as possible from these modelled geopotential anomaly tracks over a region, to demonstrate the feasibility of the Kalman filter approach for estimating TWS change from GRACE observations.

GLDAS NOAH land water storage

The NOAH (NCEP, OSU, Air Force and Office of Hydrology) land surface model (LSM) simulates surface energy and water fluxes/budgets (including soil moisture) in response to near-surface atmospheric forcing and depending on surface conditions (e.g., vegetation state, soil texture and slope) [START_REF] Ek | Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational mesoscale Eta model[END_REF]. The data used in this study are soil moisture (SM) storage values from the NOAH LSM, with the NOAH simulations being driven (parameterization and forcing) by the Global Land Data Assimilation System (GLDAS) [START_REF] Rodell | The global land data assimilation system[END_REF]). These SM estimates from GLDAS-NOAH version 2 have a spatial resolution of 1 • and a temporal resolution of 3 h. They are accessible via the Hydrology Data Holdings page at the Goddard Earth Sciences Data and Information Services Center, http://disc.sci.gsfc.nasa. gov/hydrology/data-holdings. They were cumulated on the four soil layers representative of the top 2 m of the soil from the GLDAS-NOAH model and resampled at a daily timescale over the time period 2003-2010.

Measurements of river discharge

Time series of daily water discharges from in situ gauges located in Obidos (Amazon), Ciudad Bolivar (Orinoco), Tucurui (Tocantins), and Chapeton (La Plata) are used for comparisons to daily, 10-day and monthly anomalies of GRACE-based TWS over 2003-2010. These

Global GRACE solutions

Three processing centres including the Center for Space Research (CSR), Austin, Texas, USA, the GeoForschungs-Zentrum (GFZ), Potsdam, Germany and the Jet Propulsion Laboratory (JPL), Pasadena, California, USA, and forming the Science Data Center (SDC) are in charge of the processing of the GRACE data and the production of Level-1 and Level-2 products. These products are distributed by the GFZ's Integrated System Data Center (ISDC-http://isdc.gfz-potsdam. de) and the JPL's Physical Oceanography Distributive Active Data Center (PODAAC-http://podaac-www.jpl.nasa.gov). Pre-processing of Level-1 GRACE data (i.e., positions and velocities measured by GPS, accelerometer data and KBR inter-satellite measurements) is routinely made by the SDC, as well as monthly global GRACE gravity solutions (Level-2). These latter solutions consist of time series of monthly averages of Stokes coefficients (i.e., dimensionless spherical harmonics coefficients of geopotential) developed up to a degree between 50 and 120 that are adjusted from along-track GRACE measurements. A dynamical approach, based on the Newtonian formulation of the satellite's equation of motion in an inertial reference frame, centred at the Earth's centre of mass combined with a dedicated modeling of the gravitational and non-conservative forces acting on the spacecraft, is used to compute the monthly GRACE solutions. During the estimation process, atmospheric and ocean barometric redistribution of mass variations are removed from the GRACE coefficients using ECMWF and NCEP reanalysis for atmospheric mass variations and ocean tides, as well as global ocean circulation models [START_REF] Bettadpur | CSR Level-2 processing standards document for level-2 product release 04, GRACE. The GRACE project[END_REF][START_REF] Fletchner | GFZ Level-2 processing standards document for level-2 product release 04[END_REF]. The GRACE coefficients are hence residuals that should include continental water storage change, and also signals from other geophysical phenomena and errors from correcting models and noise. The monthly GRACE solutions differ from one official provider to another due to the differences in the data processing, the choice of the correcting models and the data selection for computing the monthly averages.

A post-processing method based on independent component analysis (ICA) was applied to the Level-2 GRACE solutions from different official providers (i.e., CSR, GFZ, JPL), after a 400-km Gaussian pre-filtering. The separation is based on the assumption of statistical independence between the sources that compose the measured signals.

The estimated contributors to the observed gravity field are forced not to be correlated numerically by imposing diagonal cross-correlations. The efficiency of ICA for separating land hydrology-related signals from noise by combining Level-2 GRACE solutions has previously been demonstrated in [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF][START_REF] Frappart | An independent Component Analysis filtering approach for estimating continental hydrology in the GRACE gravity data[END_REF].

The GRGS-EIGEN GL04 models are derived from Level-1 KBRR measurements and LAGEOS 1 and 2 data for enhancement of lower harmonic degrees and using an empirical stabilization; thus, these solutions do not require any lowpass filtering to get rid of striping [START_REF] Lemoine | Temporal gravity field models inferred from GRACE data[END_REF][START_REF] Bruinsma | CNES/GRGS 10-day gravity field models (release 2) and their evaluation[END_REF] [START_REF] Carrère | Modeling the barotropic response of the global ocean to atmospheric wind and pressure forcingcomparisons with observations[END_REF]. These KBR Rate (KBRR) residuals represent the gravitational effects of non-modelled phenomena, and mainly the contribution of continental hydrology. They are easily converted into variations of along-track potential differences between the two GRACE satellites, following the energy integral method as proposed earlier by [START_REF] Jekeli | The determination of gravitational potential differences from satellite-to-satellite tracking[END_REF], [START_REF] Han | Precise estimation of in situ geopotential differences from GRACE low-low satellite-to-satellite tracking and accelerometer data[END_REF] and lately [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF]. Once corrected from known gravitational accelerations, along-track Residual Differences of Potential (RDP) mainly caused by hydrology variations in a selected continental region can reach ±0.1 m 2 /s 2 within a precision of ∼10 -3 m 2 /s 2 (see [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF]. These form the initial data set for our Kalman filter approach.

Methodology

The forward problem

We suppose that the continental region consists of M juxtaposed surface elements ( j = 1, . . . , M) of area S j located at longitude λ j and latitude θ j , and characterized by an equivalent-water height h j . The grid steps in longitude and latitude are λ and θ respectively, in the case of a geographical grid. Let Ŵ k be the N -by-M Newtonian matrix that relates the equivalent uniform water height h j and the GRACE-based potential differences Y i (i = 1, . . . , N ) for each daily period of observation k. The coefficients of this matrix are simply deduced from the inverses of the Cartesian distances ξ between the surface element heights X k and the successive positions of the GRACE satellites flying over the considered region (see [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF]. Then, the observation equation can be written as:

Ŵ k X k = Y k + v k (1)
where v k is a zero-mean process noise usually drawn from a zero-mean multivariate normal distribution with covariance matrix R k (i.e., v k ∼ N (0; R k )). In the construction of the Newtonian matrix Ŵ k , developing the inverse of the distance in sums of Legendre polynomials of n = 300-500 terms enables us to introduce the elastic Love numbers k n , to take compensation of surface water masses by elastic response of the Earth's surface into account [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF].

Additional information has to be included to find a stable mass variation estimate. In this study, we use a simple first-The K-band range (KBR) is the key science instrument of GRACE which measures the dual one-way range change of the baseline between the two coplanar, low-altitude satellites, with a precision of ∼0.1 µm/s on velocity difference, or equivalently, 10 µm in terms of line-of-sight (LOS) distance after integration versus time [START_REF] Bruinsma | CNES/GRGS 10-day gravity field models (release 2) and their evaluation[END_REF]. The average distance between the two GRACE vehicles is ∼220 km. The Level-1B KBR data represent the more precise measure of gravity variations sensed by the GRACE satellite tandem with a 10 -7 m/s accuracy, that gives access to surface water mass transfers. Coupled with the accurate 3-axis accelerometers measuring the effects of non-conservative forces acting on the satellites (i.e., atmospheric drag and solar pressure) and a priori models for correcting atmosphere and ocean mass, oceanic and solid tides and polar tides, KBR rate residuals are computed by least-squares dynamical orbit determination of 1-day-long and 5-second sampled tracks. A priori gravitational force models for numerical orbit integration of the GRACE satellites A and B prepared at the GRGS centre in Toulouse [START_REF] Bruinsma | CNES/GRGS 10-day gravity field models (release 2) and their evaluation[END_REF] are: (1) a static gravity field model EIGEN-GRGS.RL02.MEAN-FIELD to order Gauss-Markov process to model the evolution of the current estimate X k with time. More precisely, given k > 1, we use the simple prediction equation:

X k = X k-1 + w k (2)
where w k is a zero-mean process noise usually drawn from a zero-mean multivariate normal distribution with covariance matrix Q k (i.e., w k ∼ N (0; Q k )) describing the errors of the process. Although this dynamic is known to remain very crude, promising results will be shown in Sect. 4. The authors are well aware of the fact that this model can be improved by introducing an evolution model that might consist in exploiting outputs of hydrology models such as WGHM and GLDAS. Following [START_REF] Kurtenbach | Deriving daily snapshots of the Earth's gravity field from GRACE L1B data using Kalman filtering[END_REF], a more sophisticated dynamical equation than Eq. 2 has been considered through a linearized dynamic (i.e., a finite difference matrix obtained from WGHM) that is described by the prediction equation (Eq. 2). Unfortunately, the results were not significantly different from the ones shown in Sects. 5 and 6. These disappointing results may be due to the use of a linearization technique or to the presence in the observations of effects not taken into account by the model. As the Newtonian operator Ŵ k constructed from the relative Cartesian distances and applied to the unknown parameters (i.e., the equivalent-water heights) is completely linear (see [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF][START_REF] Ramillien | Constrained regional recovery of continental water mass time-variations from GRACE-based geopotential anomalies over South America[END_REF], for a given set of a priori uncertainty parameters, the Kalman filter estimate considering the entire GRACE-based RDP dataset is the same as the one obtained by summing all the separate Kalman filter estimates. This is the case if the complete RDP series are decomposed into pure dominant annual and semi-annual components plus any kind of RDP sub-dataset. In particular, the problem of recovering short-term hydrological events by tuning a priori uncertainty parameters is tackled in Sect. 6 by several simulations.

The inverse problem: sequential estimation-the Kalman filter equations

The observation equation (Eq. 1) and the prediction equation (Eq. 2) directly fit into the Kalman filter equation settings.

The Kalman filter is a recursive estimator where one needs the knowledge of the previous state and new measurements to determine the actual state [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF][START_REF] Kalman | New results in linear filtering and prediction theory[END_REF][START_REF] Evensen | Data assimilation. The ensemble Kalman filter[END_REF]). In our setting at iteration number k, the state of the estimator is represented by two variables: the current estimate X k (i.e., the vector containing the equivalentwater heights) and its error covariance matrix P k (i.e., the matrix of the uncertainties on the estimate).

The observations of the current state are used to correct the predicted variables to obtain a more precise estimate. First, the Kalman gain matrix is computed as:

K k = P k Ŵ T k Ŵ k P k Ŵ T k + R k -1 (3)
where R k is the covariance of the measurements errors considered as independent (i.e., R k = σ 2 d I , where σ 2 d is the error variance on the RDP Y k , and I represents the identity matrix). Secondly, the updated a posteriori state estimate is obtained as:

X * k = X k + K k [Y k -Ŵ k X k ] (4)
and the updated a posteriori error covariance matrix of the state estimate is computed as:

P * k = [I -K k Ŵ k ] P k (5)
At the end of this step, X k+1 is deduced according to Eq. 2: X k+1 = X * k and the resulting covariance matrix to be used at the next step is simply defined as:

P k+1 = P * k + Q k (6)
Note that for independent potential difference observations at time interval k, Q k is often chosen as the diagonal matrix:

Q k = σ 2 p I (7)
where σ 2 p is the a priori variance of process errors. The Kalman gain appears as a measure of the relative uncertainty of the measurements and the current state estimate, and it can be "tuned" to achieve particular performance. When the gain is high, it places more weight on the observations and thus follows them more closely, whereas when it is low, the Kalman filtering follows the model prediction, smoothing noise out. In the extreme case of gain of zero, the measurements are completely ignored.

Application

4.1 Inversion of potential differences simulated from hydrology model outputs

Recovery of a piece-wise time-constant water mass anomaly map

First, we need to validate the proposed sequential method of estimation by recovering "static" 30-day averaged water mass anomaly maps from GRACE-type potential differences simulated using daily outputs from WGHM, over large continental areas. For this purpose, regular 2-by-2 degree grids of equivalent-water heights over South America [90 • W-30 • W, 60 • S-20

• N] are averaged onto monthly periods from daily WGHM outputs. GRACE-type ascending and descending tracks of 5-second sampling are easily generated by the GINS software from the satellite ephemeris data and for the region of interest. Along-track potential differences are simply computed at each satellite position using Eq. 1 without noise (i.e., v k = 0). The test consists of recovering, as precisely as possible, the static 30-day averaged water mass anomaly map by applying the Kalman filter integration (Eqs. 2-7) of successive synthetic daily potential differences. The first guess is assumed to contain no hydrological signals (i.e., X k=0 = 0 for all elements of this vector for "cold start"), and in the case of no starting cross-error covariances between equivalentwater heights to recover, so we start with: P k=0 = σ 2 m I . Then, the solution is progressively built by accumulation of information from the daily GRACE RDP tracks.

As illustrated on the top row of Fig. 1, the process of estimation converges rapidly to a stable solution, which is identical to the starting water mass anomaly after 30 days of integration. Root mean square error is typically less than 10 mm after the first iteration, 1 mm after the fifth iteration, 0.1 mm after the 10th iteration, and finally 0.01 mm of equivalent-water height after a month of integration, compared to amplitude of ±300 mm of the hydrological patterns. Even short-wavelength details are revealed in the final Kalman filter solution, confirming that this noise-free recovery from simulated GRACE data is successful. A posteriori uncertainties on the fitted equivalent-water heights (i.e., square root of the diagonal elements of P k ) decrease, following the north-south direction of the tracks during the first iterations. Then they tend to be homogeneous on the whole region at the end of the fit. They are finally less than 0.1 mm, when the starting value is σ m = 1 mm.

Different intervals (1, 2 or 5 days) of integration yield to the same final solution in the case of recovering a noise-free 30-day constant map of water mass anomaly. Several combinations of σ d and σ m have been tested on the simulated case of recovery from 10 -9 to 10 -2 m 2 /s 2 , and 10 and 800 mm, respectively. A priori uncertainty on observations σ d acts as a regularization parameter and makes the inversion of the linear system for getting the Kalman gain (Eq. 3) possible. Large values of this uncertainty parameter enable substantial improvements of the solution at each integration step and thus accelerate the convergence to the final estimate. On the contrary, the convergence is slow and the final estimate is smooth when σ m is low (Fig. 2). In this latter case, low values of this parameter correspond to less weight of the GRACE observations in the Kalman gain K during the refreshing process of the solution. The final errors of a constant-time (monthly) map of water mass (i.e., difference between the estimate and water heights, starting from 1 mm down to <0.05 mm (bottom row). Final absolute error is around 0.01 mm after 30 steps (i.e., 30 days) of integration. Note the residual edge effect on the Southern boundary due to geographical truncation the reference WGHM water mass used for simulating RDP data) remain particularly small (Fig. 3), especially if σ m is large (e.g., 500 mm), suggesting that the recovery is successful.

Recovery of water mass change maps by daily updates

While the recovery of a water mass solution from RDP simulated from a 30-day constant WGHM map is successful (see previous part), the next test is to estimate cumulated solutions from RDP computed using daily-varying hydrology. GRACE RDP tracks passing over South America are simulated each day at 5-second sampled orbit positions and from the WGHM (or GLDAS) total water storage (TWS) outputs for the period 2005-2007, using Eq. 1 without noise (i.e, v k = 0). Complete series of regional solutions of water mass variation can be estimated by a Kalman filter integration strategy on daily sampling intervals and tuning a priori error parameters. Figure 4 presents regional solution when considering σ d = 0.01 m 2 /s 2 and σ m = 200 mm, and obtained by cumulating WGHM-simulated daily alongtrack differences of potential. These solutions reveal realistic seasonal amplitudes in the drainage basins of South America. Absolute errors are defined as the differences between input and recovered 2-by-2 degree water mass grids for the same day. In this case of considering very accurate RDP data (i.e., σ d very small) and in the absence of additional noise, the main hydrological structures of ±300 mm of EWH are retrieved and quite well located on the main river basins of South America. 0.001 m 2 /s 2 creates unrealistic meridian striping in these maps that increases with the number of days integrated into the current solution.

Tests of recovery of a "static" water mass grid, made in the previous Sect. 4.1.1, show that σ m controls the amplitude of information by the RDP tracks each day. In the case of recovering time-varying water mass maps with different starting values of σ m , the Kalman filter process progressively converges to the same series of daily sampling step estimates (Fig. 6), and the time taken by the integration to reach this common solution (or "spin up") from a cold start (i.e., no starting information: X k=0 = 0) is ∼3 months. This "spinup" is of 1 month for recovering a constant-time water mass In particular, Fig. 5 shows the water mass time series for the surface element corresponding to Manaus (60 • 01 ′ 32 ′′ W, 3 • 08 ′ 06 ′′ S), located in the centre of the Amazon basin. In our tests for Manaus tile and using model-simulated data, a priori error uncertainty ranges from σ d = 10 -6 m 2 /s 2 to σ d = 0.1 m 2 /s 2 while the parameter σ m = 200 mm is constant during the Kalman filter integration. As in the case of recovery of a time-constant map, the error of recovery of submonthly time-varying signals appears small when the RDP data are considered accurate (i.e., σ d very small). However, this strong assumption permits the development of numerical instabilities. Representing the series of estimated water mass maps shows that using a priori error uncertainty σ d less than Fig. 5 Multi-year time series of TWS for the surface tile number 357 centred over Manaus, that are obtained by integration of WGHM-simulated daily GRACE RDP data for several a priori error uncertainties, and a constant a priori error uncertainty on the parameters to be retrieved (i.e., the equivalent-water heights) first guess is X k=0 = 0 (i.e., "cold start") Fig. 6 Multi-year time series of TWS for the surface tile number 357 (Manaus) obtained for several a priori error uncertainties on the process (i.e., σ p ) and the equivalent-water heights (i.e., σ m ) map, as the simulated RDP data are more consistent to each other on the period of integration and reinforce the same "static" surface distribution of water mass to be retrieved (see Sect. 4.1.1). Considering daily WGHM-based water mass variations, RDP data partly contain unexpected time-varying signals that perturb the convergence of the integration process and make the "spin up" longer. As for the a priori error σ d on the RDP observations, important values of a priori process error σ p (i.e., >50 mm) generate unrealistic meridian stripes in the estimated water mass maps.

For approaching more realistic conditions of data acquisition, high-frequency random noise can be added to the modelsimulated GRACE RDP, but its effect is highly amplified if considered real and accurate signals (e.g., when a random noise of 0.001 m 2 /s 2 amplitude is added in the simulated GRACE data, the daily Kalman filter solutions are only slightly degraded but smoothed if σ d = 0.01 m 2 /s 2 or ∼10 % of hydrology-related RDP according to Ramillien et al. 2011). Unfortunately, the recovery errors reach tens of mm when the level of noise is greater than 0.01 m 2 /s 2 . To cancel the effect of noise amplification in the inversion, we found that the best compromise is to consider that σ d and σ m are about 10 -3 -10 -2 m 2 /s 2 and 100-200 mm, respectively, with value of a priori process error σ p as small as possible.

Recovery of maps from real GRACE RDP

Inversion of real potential differences for continental hydrology appears complicated because of contaminating instrumental noise, even if potential anomalies should be smooth at satellite altitude as the result of upward continuation. In return, downward continuation associated to the recovery of surface water mass anomaly amplifies the high frequencies of the signal, in particular noise of any kind. Moreover, important errors from pre-treatment such as correcting by imperfect models still remain in the orbit observations and aliased with space and time (see Sect. 2.2). Another problem arises as KBR rate residuals contain unrealistic long-term variations at fractions of the satellite revolution period. The strategy lately proposed by [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF][START_REF] Ramillien | Constrained regional recovery of continental water mass time-variations from GRACE-based geopotential anomalies over South America[END_REF] is to remove a linear trend to each potential difference track crossing the considered regions for each day before inversion, so that a recovery of medium-and high-frequency regional water mass variations is made, but at least without adding erroneous spatial long wavelengths. These missing wavelengths are added back from GRGS global solutions to complete the inverted signals afterwards.

2-by-2 degree water mass solutions after Kalman filter integration of daily GRACE RDP over South America for 2004-2010 are displayed on Fig. 7a,b. Starting parameters of Kalman filter estimation are σ d = 0.005 m 2 /s 2 and σ d = 0.01 m 2 /s 2 X k=0 = 0 (i.e., "cold start"), P k=0 = σ 2 m I with σ m = 200 mm of equivalent-water height, and no a priori process errors. When σ d = 0.005 m 2 /s 2 , meridian striping and edge effects rapidly dominate the Kalman filter solutions after ∼1 year of integration of daily RDP tracks and till the end of the total period (Fig. 7a). Increasing slightly this a priori parameter up to 0.01 m 2 /s 2 for an efficient regularization enables us to avoid the development of such unrealistic numerical instabilities in the time series of smoother solutions (Fig. 7b) where seasonal amplitudes are comparable to the ones of the global CSR solutions (Fig. 7c) (see also Tables 1,2). Besides, additional runs for testing different a priori errors have been made to compute multi-year series of regional Kalman filter solutions that show the seasonal alternating of the large water mass amplitudes in the Amazon and Orinoco river basins. Figure 8 illustrates the case of the time series for the surface tile centred on Manaus. It reveals realis-Fig. 7 Snapshots of 2-by-2 degree regional solutions of TWS over South America estimated from real GRACE RDP data and assuming σ m = 200 mm and σ p = 10 mm, and plotted at 3-month intervals revealing the dominant seasonal amplitudes of water mass over South America: Kalman filter solutions for a priori error uncertainty a σ d = 0.005 m 2 /s 2 and b σ d = 0.01 m 2 /s 2 , and c corresponding 400-km low-pass filtered CSR solutions (monthly averages) for the same periods for comparison Fig. 7 continued tic seasonal oscillations of 300-600 mm of EWH, which are clearly modulated by inter-annual variations. As expected, smooth estimates versus time are obtained when considering not precise GRACE RDP data (i.e., σ d = 0.1 m 2 /s 2 ), but we found that values of σ d lesser than 0.007 m 2 /s 2 significantly amplify the noise contained the real GRACE RDP. The time series of these estimated Kalman filter maps will be validated in Sect. 5.

Errors due to spatio-temporal aliasing

Figure 9 is a visualization of the extra information brought by the GRACE satellite tracks computed as two successive Kalman filter solutions. It shows that the Kalman filter estimate is daily updated in the very close neighbourhood under the satellite tracks-in a surface radius of about 600-800 km (e.g., see the numerical tests made in [START_REF] Ramillien | Constrained regional recovery of continental water mass time-variations from GRACE-based geopotential anomalies over South America[END_REF])just under the satellite tracks where the new information is brought. While the covariance function of the RDP data is smooth at satellite altitude (∼400 km), the one of the corresponding water mass (i.e., the source of anomaly) is also much localized on the Earth's surface due to downward continuation. This new along-track information represents a few tens of mm of EWH. Consequently, nothing is refreshed elsewhere, in areas which are not surveyed by the GRACE satellite during the considered day. This partial sampling explains that the reconstruction of a "static" map takes at least 10 days of data before the Kalman integration has a sufficient spatial coverage of the satellite tracks at the end (see previous Sect. 4.1.1). It also explains residual errors of recovery in the series of regional Kalman filter solutions, even if noisefree model-simulated RDP data are inverted (e.g., differences with model TWS values in Fig. 5), since off-track small and rapid hydrological features cannot be recovered, or partly retrieved, alternatively their signatures remain in the following Kalman filter solutions. Detection of sudden and localized hydrological events by tuning a priori error uncertainties and considering different cases of data coverage is explored with synthetic RDP data in the discussion of Sect. 6.

Validation of the Kalman filter solutions

Validation consists of confronting regional estimates of TWS obtained by integration of Kalman filter to existing GRACEbased products and independent datasets. For this purpose, we consider different sources of information presented in Sect. 2: daily WGHM and GLDAS land waters outputs, monthly global GRACE solutions computed by the official providers CSR, GFZ and JPL, as well as local records of river discharge. Comparisons with the 10-day regional solution of water mass variation over South America lately proposed by [START_REF] Ramillien | Constrained regional recovery of continental water mass time-variations from GRACE-based geopotential anomalies over South America[END_REF] have been also performed. Global Release 5 monthly solutions from CSR, JPL and GFZ have been low-pass filtered using a classical 400-km Gaussian filter [START_REF] Wahr | Time variability of the Earth's gravity field: hydrological and oceanic effects and their possible detection using GRACE[END_REF]) and an ICA approach [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF][START_REF] Frappart | An independent Component Analysis filtering approach for estimating continental hydrology in the GRACE gravity data[END_REF] to reduce striping. To make the sampling of Kalman filter solutions comparable with the other GRACE solutions in space and time, they have been averaged over both the 10-day and monthly intervals, and over the largest drainage basins of South America: Amazon (∼6 millions of km 2 ), Paranà (∼2.6 millions of km 2 ), Orinoco (∼1 million of km 2 ) and Tocantins (∼0.8 million of km 2 ).

Statistical results of the comparisons of averages are summarized in Table 1. It shows that the Kalman filter solutions are statistically closer to the 10-day regional solutions from [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF] than the monthly global solutions.

Fig. 9 Maps of the differences between successive Kalman filter solutions and the coverage of the corresponding GRACE satellite tracks used for the daily refreshing GLDAS-NOAH outputs at daily timescale. Significant time shifts and large differences of amplitude can be observed between GRACE-based and simulated TWS, especially on the Amazon, Orinoco and Tocantins basins. In these basins, a large part of the hydrological signal comes from slow reservoirs as floodplains (see [START_REF] Frappart | Surface freshwater storage in Amazon basin during the 2005 exceptional drought[END_REF] in the case of the Amazon) and groundwater (see [START_REF] Chen | Recent La Plata basin drought conditions observed by satellite gravimetry[END_REF] in the case of the La Plata basin, [START_REF] Gleeson | Water balance of global aquifers revealed by groundwater footprint[END_REF]Frappart et al. 2013b for their signatures in 10-day regional solutions) that are either not or not well modelled. As a consequence, significant differences and time-lags are likely to occur between model outputs and GRACE estimates [START_REF] Alkama | Global evaluation of the ISBA-TRIP continental hydrologic system, Part I: a twofold constraint using GRACE terrestrial water storage estimates and in-situ river discharges[END_REF] (not shown here). Larger differences in amplitude and phase Because of lower seasonal hydrological dynamics, there is less difference for the Paranà basin (i.e., <27 mm of EWH RMS) and much difference for the Tocantins river (i.e., up to 98 mm of EWH RMS). Figure 10 shows a superposition of the monthly averages of the Kalman filter solutions and the global GRACE solutions. This proves that the Kalman filter method succeeds in recovering amplitudes and phases of TWS. Indeed, they are consistent to the ones of 10-day and monthly global GRACE solutions (for instance, see Frappart et al. 2013a, b), when they are averaged over large areas. In the comparison, RMS differences are more important for the two hydrological models. Figure 11 presents the TWS time series of the Kalman filter solutions and of WGHM and are observed with GLDAS outputs, as this model (NOAH), contrary to WGHM, does not consider neither surface storage (water stored in the rivers and the floodplains) nor the routing of the surface runoff into the river network which contributes to surface storage in the downstream grid points. The contribution of surface storage to TWS represents 40-50 % in the Amazon basin (see [START_REF] Han | Dynamics of surface water storage in the Amazon inferred from measurements of inter-satellite distance change[END_REF][START_REF] Han | Movement of Amazon surface water from time-variable satellite gravity measurements and implications for water cycle parameters in land surface models[END_REF]Frappart et al. 2012;[START_REF] Paiva | Large-scale hydrologic and hydrodynamic modelling of the Amazon River basin[END_REF]. So incorporation of total runoff outputs is likely to reconcile the GLDAS TWS variations to the GRACE-based amplitudes.

Comparisons were also achieved between 10-day and monthly GRACE and daily river discharges in the four largest drainage basins of South America (i.e., Amazon, La Plata, Orinoco and Tocantins). They are presented in Fig. 12 andTable 2. An overall good agreement is found between GRACE-derived TWS and river discharge variations for all the basins (linear correlation greater than 70 % for the Kalman filter solutions for zero time-lag). Except in the case of Amazon basin where the values of correlation with discharge records reach 95 % for regional Kalman filter solutions, important linear correlations of 80-90 % are found at monthly and 10-day timescales. As TWS is the sum of the contributions all the hydrological reservoirs in a soil column (i.e., surface, soil moisture and groundwater storages), the rapid fluctuations of the discharge at daily timescale have a small impact on the TWS that contains slow Recovery from real GRACE-derived RDP confirms the low predicted seasonal amplitudes and time shifts with river discharge variations if σ d is high (Fig. 12): the a priori error matrix R dominates the other terms of Eq. 3, so the Kalman gain K is small and minimizes the weights of the input RDP data. Consequently, the current solution is constructed very slowly versus time, is not refreshed efficiently by the daily RDP data, and thus existing water mass structures in this solution are more persistent, creating delay, in other words time shifting. As the satellite tracks bring local information (Fig. 9), this impossibility of catching any information versus time has an important impact on spatial averaging over relatively small Orinoco and Tocantins River basins (see Fig. 10). changes in groundwater and other residual signals (e.g., from atmosphere and oceans mass corrections). This is why the correlations are slightly lower at daily timescale.

As they are sensitive to the input a priori error uncertainty parameters, the amplitudes of the regional Kalman filter solutions can be overestimated (e.g., if the RDP data are optimistically considered too accurate: σ d < 0.001 m 2 /s 2 ) or underestimated (see Fig. 13 for comparison of the energy spectra of solutions from daily WGHM-simulated RDP data, in particular for the smoothing parameter σ d = 0.01 m 2 /s 2 ). This sensitivity also explains important time shifts, as it is shown by previous results obtained with simulated RDP data (Figs. 5,6), when the a priori error uncertainty of the observations is high and gives smooth estimates (i.e., σ d > 0.01 m 2 /s 2 ). Tests on Kalman-type filtering of simulated error-free GRACE data show a quite rapid convergence from zero first guesses of surface mass density (i.e., an exact solution is found after a few days of data integration, as presented in Fig. 1 for South America). However, inversion of real alongtrack potential difference data remains difficult because of the presence of high-frequency noise and correcting model errors, effects of which have to be canceled in the inversion, or at least minimized. Thus, even if the convergence is fast to recover a 30-day constant map of equivalent-water heights, suitable a priori input values have to be chosen to build the Kalman gain (i.e., Eq. 3). The choice of parameters for constructing the Kalman gain is necessary to obtain the lower errors and to stabilize the final estimate to cope with noise and outliers in the real GRACE-derived RDP. Inversion of noise-free model-simulated RDP confirms that errors of recovery over South America are less than 1 mm of EWH RMS when using suitable a priori error parameters. In the presence of noise, these aliasing errors represent a few mm of EWH, as already shown by [START_REF] Encarnação | Influence of hydrology-related temporal aliasing on the quality of monthly models derived from GRACE satellite gravimetric data[END_REF]. Considering large a priori uncertainties on potential observations (e.g., σ d > 10 -2 m 2 /s 2 ) leads to recovery errors greater than 2 mm RMS in terms of equivalent-water height, as smoothing is important. Besides, considering artificially very accurate observations (e.g., σ d < 10 -4 m 2 /s 2 ) makes the Kalman gain, and thus the cumulated solution at step number k very unstable. The problem of instability is also due to the sparse sampling of the GRACE satellite tracks whose spatial coverage is not sufficient to access cumulated solutions with a daily

Discussion

While the long periods of water mass variations (>1 month) can be well recovered (see the previous section), we propose in the following discussion to perform numerical tests for exploring the ability of our sequential integration of daily WGHM-simulated RDP in recovering sub-monthly and geographically localized water mass variations. For this purpose, we model the spatio-temporal characteristics of a water mass anomaly as a time and space Gaussian-type varying function of length of a few days and placed at the centre of a 2-by-2 degree grid. Then, we simulate RDP using Eq. 1 from this error-free water mass variation model but with no use of propagated errors since no error information is really available. The challenge is to detect both the geographical location and the magnitude of such a time and space located anomaly by tuning input parameters, mainly the satellite track density (i.e., number of days of observation) and the a priori uncertainties. Once a set of these parameters are chosen and the error covariance matrices constructed, the process of integration is run (Eqs. 2-7). The errors of recovery are estimated at each step of integration as the differences between the current Kalman-cumulated solutions and the daily reference model maps. Note that the solutions may be biased towards the reference dataset. The sensitiveness of the input a priori uncertainty parameters and the impact of critical satellite track coverage are examined through different values of these a priori parameters. This trial-and-error approach helps us finding the best combination of parameters for the detection of the modelled water mass event, if the minimum RMS value of the recovery error is used as a criterion. resolution. For these limitations of not using accurate RDP data, we deemed that using σ d = 10 -2 -10 -3 m 2 /s 2 , which corresponds to the level of noise of the GRACE-based potential differences, represents a good compromise (see [START_REF] Ramillien | GRACE-derived surface mass anomalies by energy integral approach. Application to continental hydrology[END_REF][START_REF] Ramillien | Constrained regional recovery of continental water mass time-variations from GRACE-based geopotential anomalies over South America[END_REF].

Several simulations have revealed that the error of recovery increases with a priori model uncertainty σ m as well as the duration of the water mass anomaly (Figs. 14,15). This recovery error is ten times more important for water mass amplitude of 1,200 mm than considering a 200-mm water Fig. 16 Geographical variability of the error of recovery versus iteration number to the final estimate. Total integration period is 30 days. The true water mass solution is Gaussian and centred at t = 15 days, with amplitudes of 50 mm (duration: triangles: 1 day, stars: 5 day) and 500 mm (durations: squares: 1 day, circles: 5 days) numerical instabilities in the inversion, is the aliasing of fastmoving water mass events at periods shorter than the total integration period.

As the method cumulates all the hydrological signals, including the very short-term ones that can occur during the integration process, the final average corresponds to a mixing of events. We consider a synthetic Gaussian water mass anomaly located at the center of the region with varying amplitude from 50 to 500 mm of EWH, radius from 200 to 1200 km, and duration lasting from 1 to 5 days. Alongtrack RDP have been simulated from this reference model of water mass anomaly using Eq. 1, and used as input to build iteratively a Kalman filter solution. The error of recovery is then evaluated as the difference between the computed Kalman-based and reference water mass anomalies. These errors versus the number of days of integration are represented in Fig. 16. Persistent errors clearly appear after the water mass anomaly maximum occurring at day #15, and we note that they are not attenuated afterwards, even in the case of a sudden event. This error of aliasing increases drastically with amplitude of the water mass anomaly, by a factor 9 from 50 to 500 mm. In the particular case of water mass anomalies centred over the integration period, the final aliasing error should tend to zero after 30 days of integration. As it is presently built, the scheme keeps any water mass change "in memory" up to the final water mass estimate which is equivalent to an average of cumulated signals over the total period of integration. Due to the poor daily distribution of the GRACE satellite tracks, it is clear that the process cannot update the solution efficiently enough. Fortunately, in the case of long period of integration, high-frequency and mass anomaly. It reaches low values when a priori model uncertainty σ m ranges from 10 to 50 mm. In this domain of small model uncertainty, the Kalman filter strategy provides less error than least-squares integration of one day of sparse RDP data, since its advantage is, by construction, to inherit (or cumulate) information from the previous stage k -1 to build an averaged solution at the following stage k.

Obviously, the error of recovery in our final cumulated solution is lower when the data coverage at each stage is twice, indicating that the Kalman filtering is well adapted to estimate time-constant water mass map by progressive integration of RDP. Unfortunately, it is a particular case as hydrological signals surely vary in time from a day to another. This suggests that these unavoidable variations of water mass during the total period of integration (e.g., ∼30 days) create errors of time aliasing by accumulation of rapid successive hydrological events. While considering constant monthly intervals produces no noticeable aliasing error for averages over large surfaces, quantification of time aliasing errors proposed by [START_REF] Encarnação | Influence of hydrology-related temporal aliasing on the quality of monthly models derived from GRACE satellite gravimetric data[END_REF] shows that GRACE satellite orbit configuration permits an optimal detection of hydrological events of at least 11-15 days, like in the case of the 3,300-to 4,400-km-wide Zambezi river basin.

Aliasing effects polluting the final estimate

According to the previous results, even if GRACE observations contain time-varying hydrological signals, the proposed method provides estimations of water mass anomaly that are considered constant over 10-30 days. Hence, the main source of error, excepting the presence of spurious noise that creates i.e., integration of one day data), but with a twice-denser coverage (i.e., 2-day steps) of the GRACE satellite tracks at each iteration zero-mean noise should cancel out, and thus have a reduced impact on the final cumulated solution.

Reduction of aliasing error

Without including time correlations (or constraints) through hydrology model outputs, decreasing the period of integration is the only way of gaining in temporal resolution. As illustrated in Fig. 17, in the case of integration of daily RDP, the persistency of aliasing error is reduced after the maximum of anomaly is reached at k = 15. However, when the integration is made over 30 days (see Fig. 16), the amplitude of error is slightly less than independent daily integrations, thanks to heritage of useful information from the previous iteration to the next one. Obviously, coverage of the daily GRACE tracks made artificially denser would make the aliasing error decrease significantly, as presented in Fig. 18.

According to the uncertainty principle, it is not possible to benefit from temporal and spatial resolutions at the same time, as previously mentioned by [START_REF] Freeden | Spherical functions if mathematical geosciences, a scalar, vectorial and tensorial setup[END_REF], or equivalently, any representation cannot provide a precise localization of a particular event in both space and time. This is well illustrated by the previous results of [START_REF] Kurtenbach | Deriving daily snapshots of the Earth's gravity field from GRACE L1B data using Kalman filtering[END_REF][START_REF] Kurtenbach | Improved daily GRACE gravity field solutions using a Kalman smoother[END_REF] who have used a similar Kalman filter Fig. 19 Linear correlation between sub-monthly time averages of the regional Kalman filter solutions and the river discharge variations measured at four in situ stations versus the time interval used for averaging these two datasets These long time series of daily-step Kalman filter solutions have been validated by comparing them to GRACE-based and independent datasets (i.e., model outputs and in situ discharge observations). The seasonal amplitudes of TWS and their inter-annual modulations are well restored on the daily sampling solutions compared with regional/global GRACE solutions and model outputs. High correlations (>0.7) were found between daily TWS and in situ discharge data in the four largest drainage basins of South America, but were generally lower than these obtained at 10-day or monthly averages.

While long periods of the hydrological variations are well recovered by the sequential accumulation (e.g., by choosing suitable high values for σ m to allow strong seasonal signals to be easily retrieved), the construction of sub-monthly hydrological events remains problematic because of the poor per-day coverage of the GRACE satellite tracks. According to the simulations we made, the refreshing process is efficient for accumulation of RDP information of 10-15 days, but fails to reach the daily resolution of the Kalman filter solutions and even a resolution of a few days. The detection of sub-monthly events can be slightly improved by tuning a priori error uncertainty parameters; however, the spatial distribution of the satellite tracks is the most limiting factor. approach to estimate the dominant trend, annual and semiannual cycles of the low-degree smooth spherical harmonics of the geopotential from daily GRACE data, with no access to the shorter time periods of the water mass variations. Comparing sub-monthly variations of Kalman filter solutions and in situ river discharge in terms of residuals (i.e., difference between average of TWS and discharge over t varying from 1 to 30 days and monthly respective quantities) reveals that their linear correlation increases with the interval of averaging from 10 to 15 days at a rate of 0.5-0.7 % per day (Fig. 19). This correlation remains greater than 65-70 % for averaging intervals greater than 10 days. Below this interval of time, the correlation rate flattens, suggesting that the regional Kalman filter solutions are not refreshed enough by the sparse coverage of the daily GRACE tracks to be compared to the local discharge measurements.

Conclusion

A new sequential approach for estimating regional 2-by-2 degree variations of water mass has been developed by integrating successive days of GRACE-type potential difference data. This method is based on cumulating satellite information to build progressively stable regional time averages over large continental areas. Recovering piece-wise timeconstant maps of water mass is probably the most adapted case for the proposed method, as the permanent hydrological structures are reinforced during the iterative Kalman filter estimation. This sequential estimation was applied to estimate time-varying regional solutions over South America by daily updates using real GRACE potential differences. Series of regional cumulated solutions were derived at daily intervals by refreshing using real daily (and thus sparse) GRACE RDP, and testing ranges of sensitive a priori error uncertainties on the parameters and observations, before being analysed and confronted to other GRACE-based datasets.

Fig. 1

 1 Fig. 1 Recovery of monthly mass variations over South America for March 2005 from simulated along-track RDP. Units are mm of equivalent-water height. Estimated regional maps obtained by accumulation of 1-day data that converge to a stable solution (top row), and the decreasing of the corresponding a posteriori uncertainties on the

Fig. 2

 2 Fig. 2 Convergence analysis of the recovery using a priori potential anomaly standard deviation of σ d = 10 -9 m 2 /s 2 (i.e., exact observations) and different a priori parameter uncertainties σ m : 1 mm (circles), 10 mm (triangles), 100 mm (squares), 500 mm (stars)

Fig. 4

 4 Fig. 4 Daily 2-by-2 degree regional maps of TWS over South America plotted at monthly intervals [units: mm of equivalent-water height (EWH)]

Fig. 7

 7 Fig. 7 continued

Fig. 10

 10 Fig. 10 Time series of TWS for the largest drainage basins of South America obtained by averaging monthly CSR (blue), GFZ (red), JPL (green) and daily Kalman filter solutions (black)

Fig. 11

 11 Fig. 11 Time series of TWS over the largest drainage basins of South America computed using WGHM model (dashed line) and GLDAS model (dots) outputs, as well as daily Kalman filter solutions (solid line)

Fig. 12

 12 Fig. 12 River discharge variations (grey) measured at different in situ stations: a Obidos (Amazon), b Chapeton (Paranà), c Ciudad Bolivar (Orinoco) and d Tucurui (Tocantins), as well as time series of

Fig. 14 Fig. 15

 1415 Fig. 14 Maximum error of recovery versus a priori parameter σ m after 30 days of Kalman filter integration in the case of a simulated 200-mm amplitude water mass anomaly centred at time k = 15, and assuming an a priori error uncertainty of the observations of 0.001 m 2 /s 2

Fig. 17

 17 Fig. 17Same as Fig.16, but the total integration period is 1 day, instead of 30 days, to gain in temporal resolution and avoid persistent signals

  Fig. 18 Same as Fig. 17(i.e., integration of one day data), but with a twice-denser coverage (i.e., 2-day steps) of the GRACE satellite tracks at each iteration

  • 39 ′ 25 ′′ W; 1 • 55 ′ 23 ′′ S), Manacapuru (60 • 36 ′ 32 ′′ W; 3 • 18 ′ 58 ′′ S), Fazenda Vista Alegre (60 • 01 ′ 34 ′′ W; 4 • 53 ′ 53 ′′ S), Porto Velho (63 • 56 ′ 46 ′′ W; 8 • 47 ′ 59 ′′ S) and Tucurui (49 • 40 ′ 59 ′′ W; 3 • 46 ′ 59 ′′ S) gauges; and (3) the Argentinian water agency (Instituto Nacional del Agua-INA) through the online database Base de Datos Hidrológica Integrada (BDHI-http://www.hidricosargentina.gov.ar/) for the Chapeton station (60 • 16 ′ 59 ′′ W; 31 • 34 ′ 26 ′′ S).

	in situ
	records were downloaded for the period 2003-2010 from:
	(1) the Venezuelian water agency (Instituto Nacional de
	Meteorologia e Hidrologia-INAMEH) for Ciudad Bolivar
	(63 • 36 ′ 29 ′′ W; 8 • 26 ′ 20 ′′ N); (2) the hydrological information
	system Hidroweb (http://hidroweb.ana.gov.br/) of the Brazil-
	ian water agency (Agência Nacional de Aguas-ANA) for
	Obidos (55

Table 1

 1 Root mean square of the absolute differences between daily Kalman filter solutions and different TWS datasets for spatial averages over the main drainage basins of South America [units: mm of equivalent-water height (EWH)]

	Amazon	Paranà	Orinoco	Tocantins

Acknowledgments We thank Dr. Lucia Seoane from GET laboratory in Toulouse for having provided us daily GRACE-type orbit data over continental areas. This research work was partly funded by the TOSCA/CNES "Surcharge & Déformation" project lead by P. Gégout.