Frédéric Frappart
email: frederic.frappart@lmtg.obs-mip.fr

Guillaume Ramillien
email: guillaume.ramillien@dtp.obs-mip.fr

Marc Leblanc
email: marc.leblanc@jcu.edu.au

Sarah O Tweed
email: sarah.tweed@jcu.edu.au

Marie-Paule Bonnet
email: marie-paule.bonnet@ird.fr

Philippe Maisongrande
email: philippe.maisongrande@cnes.fr

An Independent Component Analysis filtering approach for estimating continental hydrology in the GRACE gravity data

An approach based on Independent Component Analysis (ICA) has been applied on a combination of monthly GRACE satellite solutions computed from official providers (CSR, JPL and GFZ), to separate useful geophysical signals from important striping undulations. We pre-filtered the raw GRACE Level-2 solutions using Gaussian filters of 300, 400, 500-km of radius to verify the non-gaussianity condition which is necessary to apply the ICA. This linear inverse approach ensures to separate components of the observed gravity field which are statistically independent. The most energetic component found by ICA corresponds mainly to the contribution of continental water mass change. Series of ICA-estimated global maps of continental water storage have been produced over 08/2002-07/2009. Our ICA estimates were compared with the solutions obtained using other post-processings of GRACE Level-2 data, such as destriping and Gaussian filtering, at global and basin scales. Besides, they have been validated with in situ measurements in the Murray Darling Basin. Our computed ICA grids are consistent with the different approaches. Moreover, the ICA-derived time-series of water masses showed less north-south spurious gravity signals and improved filtering of unrealistic hydrological features at the basin-scale compared with solutions obtained using other filtering methods.

Introduction

Continental water storage is a key component of global hydrological cycles and plays a major climate system via controls over water, energy and biogeochemical fluxes.

In spite of its importance, the total continental water storage is not well-known at regional and global scales because of the lack of in situ observations and systematic monitoring of the groundwaters [START_REF] Alsdorf | Tracking fresh water from space[END_REF].

The Gravity Recovery and Climate Experiment (GRACE) mission provides a global mapping of the time-variations of the gravity field at an unprecedented resolution of ~400 km and a precision of ~1 cm in terms of geoid height. Tiny variations of gravity are mainly due to redistribution of mass inside the fluid envelops of the Earth (i.e., atmosphere, oceans and continental water storage) from monthly to decade timescales (Tapley et al., 2004).

Pre-processing of GRACE data is made by several providers (University of Texas, Centre for Space Researh -CSR, Jet Propulsion Laboratory -JPL, GeoForschungsZentrum -GFZ and Groupe de Recherche en Géodésie Spatiale -GRGS) which produce residual GRACE spherical harmonic solutions that mainly represent continental hydrology as they are corrected from known mass transfers using ad hoc oceanic models (i.e., Toulouse Unstructured Grid Ocean model 2D -T-UGOm 2D) and atmospheric reanalyses from National Centers for Environmental Prediction (NCEP) and European Centre for Medium Weather Forecasting (ECMWF). Unfortunately these solutions suffer from the presence of important north-south striping due to orbit resonance in spherical harmonics determination and aliasing of short-time phenomena which are geophysically unrealistic.

Since its launch in March 2002, the GRACE terrestrial water storage anomalies have been increasingly used for large-scale hydrological applications (see [START_REF] Ramillien | Detection of continental hydrology and glaciology signals from GRACE: a review[END_REF]Schmitt et al., 2008 for reviews). They demonstrated a great potential to monitor extreme hydrological events [START_REF] Andersen | GRACE-derived terrestrial water storage depletion associated with the 2003 European heat wave[END_REF][START_REF] Seitz | Signals of extreme weather conditions in Central Europe in GRACE 4-D hydrological mass variations[END_REF][START_REF] Chen | 2005 drought event in the Amazon River basin as measured by GRACE and estimated by climate models[END_REF], to estimate water storage variations in the soil [START_REF] Frappart | Interannual variations of river water storage from a multiple satellite approach: A case study for the Rio Negro River basin[END_REF], the aquifers [START_REF] Rodell | Estimating groundwater storage changes in the Mississippi River basin (USA) using GRACE[END_REF][START_REF] Strassberg | Comparison of seasonal terrestrial water storage variations from GRACE with groundwater-level measurements from the High Plains Aquifer (USA)[END_REF][START_REF] Leblanc | Basin scale, integrated observations of the early 21st century multiyear drought in southeast Australia[END_REF] and the snowpack [START_REF] Frappart | Evolution of high-latitude snow mass derived from the GRACE gravimetry mission (2002-2004[END_REF]in press), and hydrological fluxes, such as basin-scale evapotranspiration (Rodell et al., 2004a;Ramillien et al. , 2006a) and discharge (Syed et al. , 2009).

Because of this problem of striping that limits geophysical interpretation, different postprocessing approaches for filtering GRACE geoid solutions have been proposed to extract useful geophysical signals (see [START_REF] Ramillien | Detection of continental hydrology and glaciology signals from GRACE: a review[END_REF]Schmidt et al., 2008 for reviews).

These include the classical isotropic Gaussian filter [START_REF] Jekeli | Alternative methods to s Tech[END_REF], various optimal filtering decorrelation of GRACE errors [START_REF] Han | Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement[END_REF][START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF]Swenson and Wahr, 2006;[START_REF] Sasgen | Wiener optimal filtering of GRACE data[END_REF][START_REF] Kusche | Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models[END_REF][START_REF] Klees | The design of an optimal filter for monthly GRACE gravity models[END_REF], as well as statistical constraints on the time evolution of GRACE coefficients [START_REF] Davis | A statistical filtering approach for Gravity Recovery and Climate Experiment (GRACE) gravity data[END_REF] or from global hydrology models [START_REF] Ramillien | Time variations of land water storage from the inversion of 2-years of GRACE geoids[END_REF]. However, these filtering techniques remain imperfect as they require input non-objective a priori information which are most of the time simply tuned by hand (e.g. choosing the cutting wavelength while using the Gaussian filtering) or based on other rules-of-thumb.

We propose another post-processing approach of the Level-2 GRACE solutions by considering completely objective constraints, so that the gravity component of the observed signals is forced to be uncorrelated numerically using an Independent Component Analysis (ICA) technique. This approach does not require a priori information except the assumption of statistical independence of the elementary signals that compose the observations, i.e., geophysical and spurious noise. The efficiency of ICA to separate gravity signals and noise from combined GRACE solutions has previously been demonstrated on one month of Level-2 solutions [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF]. In this paper, we use this new statistical linear method to derive complete time series of continental water mass change.

The first part of this article presents the datasets used in this study: the monthly GRACE solutions to be inverted by ICA and to be used for comparisons, and the in situ data used for validation of our estimates over the Murray Darling drainage basin (~ 1 million of km²). This region has been selected for validation because of available dense hydrological observations. The second part outlines the three steps of the ICA methodology. Then the third and fourth parts present results and comparisons with other post-processed GRACE solutions at global and regional scales, and in situ measurements for the Murray Darling Basin respectively.

Error balance of the ICA-based solutions is also made by considering the effect of spectrum truncation, leakage and formal uncertainties.

Datasets

The GRACE data

The GRACE mission, sponsored by National Aeronautics and Space Administration (NASA) and Deutsches Zentrum für Luft-und Raumfahrt (DLR), has been collecting data since mid-2002. Monthly gravity models are determined from the analysis of GRACE orbit perturbations in terms of Stokes spherical harmonic coefficients, i.e., geopotential or geoid heights. The geoid is an equipotential surface of the gravity field that coincide with mean sea level. For the very first time, monthly global maps of the gravity time-variations can be derived from GRACE measurements, and hence, to estimate the distribution of the change of mass in the Earth system.

The Level-2 raw solutions

The Level-2 raw data consist of monthly estimates of geo-potential coefficients adjusted for each 30-day period from raw along-track GRACE measurements by different research groups (i.e., CSR, GFZ and JPL). These coefficients are developed up to a degree 60 (or spatial resolution of 333 km) and corrected for oceanic and atmospheric effects [START_REF] Bettadpur | CSR level-2 processing standards document for level-2 product release 0004[END_REF] to obtain residual global grids of ocean and land signals corrupted by a strong noise. These data are available at: ftp://podaac.jpl.nasa.gov/grace/.

The destriped and smoothed solutions

The monthly raw solutions (RL04) from CSR, GFZ, and JPL were destriped and smoothed by [START_REF] Chambers | Evaluation of new GRACE time-variable gravity data over the ocean[END_REF] for hydrological purposes. These three datasets are available for several averaging radii (0, 300 and 500 km on the continents and 300, 500 and 750 km on the oceans) at ftp://podaac.jpl.nasa.gov/tellus/grace/monthly.

In this study, we used the Level-2 RL04 raw data from CSR, GFZ and JPL, that we filtered with a Gaussian filter for radii of 300, 400 and 500 km, and the destriped and smoothed solutions for the averaging radii of 300 and 500 km over land.

The hydrological data for the Murray Darling Basin

In the predominantly semiarid Murray Darling Basin, most of the surface water is regulated using a network of reservoirs, lakes and weirs [START_REF] Kirby | The shared water resources of the Murray-Darling Basin. Murray-Darling Basin Commission[END_REF] and the surface water stored in these systems represent most of the total surface water (SW) present across the basin. A daily time series of the total surface water storage in the network of reservoirs, lakes, weirs and in-channel storage was obtained from the Murray-Darling Basin Commission and the state governments from January 2000 to December 2008.

In the Murray Darling Basin, we derived monthly soil moisture (SM) storage values for the basin from January 2000 to December 2008 from the NOAH land surface model [START_REF] Ek | Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational mesoscale Eta model[END_REF], with the NOAH simulations being driven (parameterization and forcing) by the Global Land Data Assimilation System [START_REF] Rodell | The Global Land Data Assimilation System[END_REF]. The NOAH model 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 NOAH model outputs of soil moisture estimates have a 1° spatial resolution and, using four soil layers, are representative of the top 2 m of the soil.

In situ estimates of annual changes in the total groundwater storage (GW) across the drainage basin were obtained from an analysis of groundwater levels observed in government monitoring bores from 2000 to 2008. Compared to earlier estimates by [START_REF] Leblanc | Basin scale, integrated observations of the early 21st century multiyear drought in southeast Australia[END_REF], the groundwater estimates presented in this paper provide an update of the in situ water level and a refinement of the distribution of the aquifers storage capacity.

Assuming that (1) the shallow aquifers across the Murray-Darling drainage basin are hydraulically connected and that (2) at a large scale the fractured aquifers can be assimilated to a porous media, changes in groundwater storage across the area can be estimated from observations of groundwater levels (e.g., [START_REF] Rodell | Estimating groundwater storage changes in the Mississippi River basin (USA) using GRACE[END_REF][START_REF] Strassberg | Comparison of seasonal terrestrial water storage variations from GRACE with groundwater-level measurements from the High Plains Aquifer (USA)[END_REF].

Variations in groundwater storage ( GW ) were estimated from in situ measurements as:

GW y S S H (1)
where S y is the aquifer specific yield (%) and H is the groundwater level (L - [START_REF] Macumber | Groundwater flow and resource potential in the Bridgewater and Salisbury West GMAs[END_REF][START_REF] Cresswell | Assessment of salinity management options for Kyeamba Creek, New South Wales: Data analysis and groundwater modelling[END_REF]Hekmeijer and Dawes, 2003a;CSIRO, 2008); from 10 to 15% for the shallow sandy clay unconfined aquifer group [START_REF] Macumber | Groundwater flow and resource potential in the Bridgewater and Salisbury West GMAs[END_REF]Urbano et al., 2004); and from 1 to 10% for the fractured rock aquifer group [START_REF] Cresswell | Assessment of salinity management options for Kyeamba Creek, New South Wales: Data analysis and groundwater modelling[END_REF][START_REF] Hekmeijer | Assessment of salinity management options for Axe Creek, Victoria: Data analysis and groundwater modelling[END_REF][START_REF] Smitt | Assessment of salinity management options for the Brymaroo catchment, South-eastern Queensland[END_REF][START_REF] Petheram | Testing in class variability of groundwater systems: local upland systems[END_REF]. In situ estimates of changes in GW storage are calculated using the spatially averaged change in annual groundwater level across each type of unconfined aquifer group and the mean value of the specific yield for that group using (Eq. 1); while the range of possible values for the specific yield was used to estimate the uncertainty.

Groundwater changes in the deep, confined aquifers (mostly GAB and Renmark aquifers) are either due to: 1) a change in groundwater recharge at the unconfined outcrop; 2) shallow pumping at the unconfined outcrop or 3) deep pumping in confined areas for farming (irrigation and cattle industry). GRACE TWS estimates accounts for all possible sources of influence, while GW in situ estimates only include those occurring across the outcrop. Total pumping from the deep, confined aquifers was estimated to amount to -0.42 km 3 .yr -1 in 2000 (Ife and Skelt, 2004), while groundwater pumping across the basin was -1.6 km 3 in 2002

2003 [START_REF] Kirby | The shared water resources of the Murray-Darling Basin. Murray-Darling Basin Commission[END_REF]. To allow direct comparison between TWS and in-situ GW estimates, pumping from the deep aquifers was added to the in-situ GW time series assuming the -0.42 km 3 .yr -1 pumping rate remained constant during the study period. 3. Methodology

ICA-based filter

ICA is a powerful method for separating a multivariate signal into subcomponents assuming their mutual statistical independence [START_REF] Comon | Independent component analysis: a new concept?[END_REF][START_REF] De Lathauwer | An introduction to independent component analysis[END_REF]. It is commonly used for blind signals separation and has various practical applications [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF], including telecommunications (Ritsaniemi and Joutsensalo, 1999;[START_REF] Cristescu | Fading Channel Estimation by Mutual Information Minimization for Gaussian Stochastic Processes[END_REF], medical signal processing (Vigário, 1997;[START_REF] Van Hateren | Independent component filters of natural images compared with simple cells in primary visual cortex[END_REF], speech signal processing [START_REF] Stone | Independent Component Analysis: A Tutorial Introduction[END_REF], and electrical engineering [START_REF] Gelle | Blind source separation: a tool for rotating machine monitoring by vibration analysis[END_REF][START_REF] Pöyhönen | Independent Component Analysis of vibrations for fault diagnosis of an induction motor[END_REF].

Assuming that an observation vector y collected from N sensors is the combination of P (N P) independent sources represented by the source vector x, the following linear statistical model can be considered:

y Mx (2)
where M is the mixing matrix whose elements m ij (1 i N, 1 j P) indicate to what extent the j th source contribute to the i th observation. The columns {m j } are the mixing vectors.

The goal of ICA is to estimate the mixing matrix M and/or the corresponding realizations of the source vector x, only knowing the realizations of the observation vector y, under the assumptions [START_REF] De Lathauwer | An introduction to independent component analysis[END_REF]:

1) the mixing vectors are linearly independent,

2) the sources are statistically independent.

The original sources x can be simply recovered by multiplying the observed signals y with the inverse of the mixing matrix ICA consists of three numerical steps. The first step of ICA is to centre the observed vector, i.e., to substract the mean vector m E y to make y a zero mean variable. The second step consists in whitening the vector y to remove any correlation between the components of the observed vector. In other words, the components of the white vector ỹ have to be uncorrelated and their variances equal to unity. Letting t C E yy be the correlation matrix of the input data, we define a linear transform B that verifies the two following conditions:

y By (4) and:

t P E yy I (5)
where I P is identity matrix of dimension P x P.

This is easily accomplished by considering:

1 2 B C (6)
The whitening is obtained using an EVD of the covariance matrix C:

t C EDE (7)
where E is the orthogonal matrix of the eigenvectors of C and D is the diagonal matrix of its eigenvalues. D=diag(d 1 d P ) as a reduction of the dimension of the data to the number of independent components (IC) P is performed, discarding the too small eigenvalues.

For the third step, an orthogonal transformation of the whitened signals is used to find the separated sources by rotation of the joint density. The appropriate rotation is obtained by maximizing the non-normality of the marginal densities, since a linear mixture of independent random variables is necessarily more Gaussian than the original components.

Many algorithms of different complexities have been developed for ICA [START_REF] Stone | Independent Component Analysis: A Tutorial Introduction[END_REF]. The FastICA algorithm, a computationally highly efficient method for performing the estimation of ICA [START_REF] Hyvärinen | Independent Component Analysis: Algorithms and Applications[END_REF] has been considered to separate satellite gravity signals. It uses a fixed-point iteration scheme that has been found to be 10 to 100 times faster than conventional gradient methods for ICA [START_REF] Hyvärinen | Fast and Robust Fixed-Point Algorithms for Independent Component Analysis[END_REF].

We used the FastICA algorithm (available at http://www.cis.hut.fi/projects/ica/fastica/) to unravel the IC of the monthly gravity field anomaly in the Level-2 GRACE products. We previously demonstrated, on a synthetic case, that land and ocean mass anomalies are statistically independent from the north-south stripes using information from land and ocean models and simulated noise [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF]. Considering that the GRACE Level-2 products from CSR, GFZ and JPL are different observations of the same monthly gravity anomaly and, that the land hydrology and the north-south stripes are the independent sources, we applied this methodology to the complete 2002-2009 time series. The raw Level-2 GRACE solutions present Gaussian histograms which prevent the successful application of the ICA method. To ensure the non-Gaussianity of the observations, the raw data have been preprocessed using Gaussian filters with averaging radii of 300, 400 and 500 km as in [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF].

Time-series of basin-scale total water storage average

For a given month t, the regional average of land water volume V(t) (or height h(t)) over a

given river basin of area A is simply computed from the water height h j , with (expressed in terms of mm of equivalent-water height) inside A, and the elementary surface

R e 2 sin j : 2 ( ) ( , , )sin j e j j j j A V t R h t (8) 2 ( ) ( , , )sin e j j j j j A R h t h t A (9)
where j and j are co-latitude and longitude of the j th point, and are the grid steps in longitude and latitude respectively (generally = ). In practice, all points of A used in (Eq. 8 and 9) are extracted for the eleven drainage basins masks at a 0.5° resolution provided by [START_REF] Oki | Design of Total Runoff Integrating Pathways (TRIP) -A global river channel network[END_REF], except for the Murray Darling Basin where we used basin limits from [START_REF] Leblanc | Basin scale, integrated observations of the early 21st century multiyear drought in southeast Australia[END_REF].

Regional estimates of formal error

As ICA provides separated solutions which have Gaussian distributions, the variance of the regional average for a given basin is:

2 2 1 2 L k k formal L ( 10
)
where formal is the regional formal error, k is the formal error at a grid point number k, and L is the number of points used in the regional averaging.

If the points inside the considered basin are independent, this relation is slightly simplified:

k formal L (11)

Frequency cut-off error estimates

Error in frequency cut-off represents the loss of energy in the short spatial wavelength due to the low-pass harmonic decomposition of the signals that is stopped at the maximum degree

N 1 .
For the GRACE solution separated by ICA; N 1 =60, thus the spatial resolution is limited and stopped at ~330 km by construction. This error is simply evaluated by considering the difference of reconstructing the remaining spectrum between two cutting harmonic degrees N 1

and N 2 , where N 2 > N 1 and N 2 should be large enough compared to N 1 (e.g., N 2 =300 in study):

2 1 2 1 0 0 N N N truncation n n n n n n N (12)
using the scalar product

0 n n nm nm nm nm m C A S B ( 13
)
where A nm and B nm are the harmonic coefficients of the considered geographical mask, and C nm and S nm are the harmonic coefficients of the water masses.

Leakage error estimates

We define « leakage » as the portion of signals from outside the considered geographical region that pollutes the estimates. By construction, this effect can be seen as the limitation of the geoid signals degree in the spherical harmonics representation. For each basin and at each period of time, leakage is simply computed as the average of outside values by using an « inverse » mask, which is 0 and 1 in and out of the region respectively, developed in spherical harmonics and then truncated at degree 60. This method of computing leakage of continental water mass has been previously proposed for the entire continent of Antarctica (Ramillien et al., 2006b), which revealed that the seasonal amplitude of this type of error can be quite important (e.g. up to 10% of the geophysical signals). In case of no leakage, this average should be zero (at least, it decreases with the maximum degree of decomposition). However, the maximum leakage of continental hydrology remains in the order of the signals magnitude itself.

Results and discussion

ICA-filtered land water solutions

The methodology presented in [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF] has been applied to the Level-2 RL04 raw monthly GRACE solutions from CSR, GFZ and JPL, preprocessed using a Gaussian filter with a radius of 300, 400 and 500 km, over the period July 2002 to July 2009. The results of this filtering method is presented in Fig. 1 for four different time periods (March and September 2006, March 2007and March 2008) using the GFZ solutions Gaussian-filtered with a radius of 400 km. Only the ICA-based GFZ solution is presented since, for a specific radius, the ICA-based CSR, GFZ, and JPL solutions only differ from a scaling factor for each specific component. The ICA-filtered CSR, GFZ, and JPL solutions are obtained by multiplying the j th IC with the j th mixing vector (2). As the last two modes correspond to the north south stripes, we present their sum in Fig. 1.

The first component is clearly ascribed to terrestrial water storage with variations in the range of 450 mm of Equivalent Water Thickness (EWT) for an averaging radius of 400 km. The larger water mass anomalies are observed in the tropical regions, i.e., the Amazon, the Congo, the Ganges and the Mekong Basins, and at high latitudes in the northern hemisphere. The components 2 and 3 correspond to the north-south stripes due to orbits. They are smaller than the first component by a factor of 3 or 4 as previously found [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF].

The FastICA algorithm was unable to retrieve realistic patterns and/or amplitudes of TWSderived from GRACE data preprocessed using a Gaussian filter with a radius 300 km for Some of these dates, such as the period between June and November 2004, correspond to deep resonance between the satellites caused by an almost exact repeat of the orbit, responsible for a significantly poorer accuracy of the monthly solutions [START_REF] Chambers | Evaluation of new GRACE time-variable gravity data over the ocean[END_REF]. As ICA is based on the assumption of independence of the sources, if the sources exhibit similar statistical distribution, the algorithm is unable to separate them.

A classical measure of the peakiness of the probability distribution is given by the kurtosis.

The kurtosis K y is dimensionless fourth moment of a variable y and classically defined as:

4 2 2 y E y K E y (14)
If the probability density function of y is purely Gaussian, its kurtosis has the numerical value of 3. In the followings, we will consider the excess of kurtosis (K y -3) and refer to the kurtosis as it is commonly done. So a variable y will be Gaussian if its kurtosis remains close to 0.

The time series of the kurtosis of the sources separated using ICA are presented in Fig. 2 We also observed that the number of time steps with only one IC (the outputs are identical to the inputs, i.e., no independent sources are identified and hence no filtering was performed)

increases with the radius of the Gaussian filter (none at 300 km, 2 at 400 km, 7 at 500 km).

In the following, as the ICA-derived TWS with a Gaussian prefiltering of 300 km, exhibits an important gap of 6 months in 2004, we will only consider the solutions obtained after a preprocessing with a Gaussian filter for radii of 400 and 500 km (ICA400 and ICA500).

Global scale comparisons

Global scale comparisons have been achieved with commonly-used GRACE hydrology preprocessings: the Gaussian filter [START_REF] Jekeli | Alternative methods to s Tech[END_REF] and the destriping method (Swenson and Wahr, 2006) for several smoothing radii.

ICA versus Gaussian-filtered solutions

Advantages of extracting continental hydrology using ICA after a simple Gaussian filtering have to be demonstrated for the complete period of availability of the GRACE Level-2 dataset, as it was for one period of GRACE Level-2 data in [START_REF] Frappart | Denoising satellite gravity signals by Independent Component Analysis[END_REF]. Numerical tests of comparisons before and after ICA have been made to show full utility of considering a post-treatment by ICA for signals separation. For the period 2003-2008, we present correlation (Fig. 3) and RMS (Fig. 4) maps computed between the ICA400 (respectively ICA500) and Gaussian-filtered solutions with a radius of 400 km (500 km), named in the followings G400 (G500). High correlation coefficients are generally observed over land (greater than 0.9), increasing with the smoothing radius, especially over areas with large hydrological signals, i.e., Amazon, Congo and Ganges basins, boreal regions. On the contrary, low correlation values, structured as stripes, are located over arid and semi-arid regions (southwest of the US, Sahara, Saudi Arabia, Gobi desert, centre of Australia), especially for GFZ and JPL solutions. These important RMS differences between Gaussian and ICA-based solutions reveal that the GRACE signals still contains remaining stripes after the Gaussian filtering. This justifies that extracting the useful continental hydrology signals requires a further processing. For this purpose, ICA succeeds in isolating this noise in its second and third components (as illustrated in Figure 1).

In Fig. 4, we observe that the spatial distribution of the RMS between ICA solutions and Gaussian solutions presents north-south stripes with values generally lower than 30 mm, except for some spots between (20 30)° of latitudes on the GFZ (Fig. 4c andd) and JPL (Fig. 4e ande) solutions. The ICA approach allows the filtering of remnant stripes present in the Gaussian solutions, especially for GFZ and JPL (Fig. 4c to f). These unrealistic structures (stripes and spots), which correspond to resonances in the orbit of the satellites, are clearly filtered out using the ICA approach (compare with Fig. 1). A more important smoothing due to a larger radius caused a decrease of the RMS between ICA and Gaussian solutions (Fig 4a, c and e). The RMS can reach 100 mm between ICA400 and G400 and only 65 mm between ICA500 and G500 for the GFZ solutions.

ICA versus destriped and smoothed solutions

Similar to Fig. 3 and4, we present correlation (Fig. 5) and RMS (Fig. 6) maps over the period [2003][2004][2005][2006][2007][2008]; between ICA400 (ICA500 respectively) and destriped and smoothed solutions with radii of 300 km (DS300) and 500 km (DS500), made available by [START_REF] Chambers | Evaluation of new GRACE time-variable gravity data over the ocean[END_REF]. The correlation maps between ICA400 (respectively ICA500) and DS300 (DS500) exhibit very similar patterns compared with those presented in Fig. 3, except in the McKenzie Basin and Nunavut (northern Canada) where low correlation or negative correlations for GFZ and JPL solutions (Fig. 5c to f) are found.

The RMS differences between ICA and destriped and smoothed solutions present low values above 30° (or below -30°) of latitude (< 30 mm of EWT), with the exceptions of Nunavut in northern Canada and extreme values (up to 100 mm) in the tropics. This is clearly due to the low performance of the destriping method in areas close to Equator as previously mentioned by Swenson and Wahr (2006) and [START_REF] Klees | The design of an optimal filter for monthly GRACE gravity models[END_REF]. These extrema are especially present in the DS300 GFZ and JPL solutions. Some secondary maxima (up to 80 mm) are also noticed for the GFZ and the JPL solutions. These important differences correspond to north-south stripes that still appear in the destriped and smoothed solutions despite the filtering process (and that can be filtered out by applying an ICA approach see the results for DS300 GFZ solution of March 2006 in Fig. 7). For an averaging radius of 500 km, the RMS between ICA500 and DS500 is lower than 20 to 30 mm, except for the Nunavut (Fig. 6b,d,f), along the Parana stream (40 to 50 mm in the CSR solutions Fig. 6b), and along the Amazon and Parana streams (60 to 80 mm in the GFZ solutions Fig. 6d).

These low spatial correlations suggest Gaussian-ICA provides at the least equivalent results on the continents to the smoothing-destriping method. Besides, it is interesting that both approaches are based on a pre-Gaussian filtering. Short wavelength differences between the maps obtained separetely using ICA and destriping reveals the limitation of the destriping which generates artefacts in the tropics.

Trend comparisons

The effects of the stripes on the trends estimated using GRACE-based TWS is supposed to be significant. From the time-series of TWS anomaly grids derived from GRACE (using equation 1), the temporal trend, seasonal and semi-annual amplitudes were simultaneously fitted by least-square adjustment at each grid point over the period 2003-2008 (see Frappart et al., in press, for details) (see Frappart et al., in press, for details). We present in Fig. 8 trends of TWS for ICA400 (Fig. 8a) and ICA500 (Fig. 8b), for G400 (Fig. 8c) and G500 (Fig. 8d), and for DS300 (Fig. 8e) and DS500 (Fig. 8f).

The trend estimates exhibit large differences in spatial patterns and the amplitude of the signal, especially between ICA and other processing methods. The most significant differences over land (except Antarctica) are located at high latitudes, over Scandinavia, and the Laurentide region, in the northeast of Canada. The ICA400 and ICA500 solutions present negative trends of TWS (Fig. 8a andb), whereas the G400 and G500 and the DS300 and DS500 present large positive trends (Fig. 8c to f). These two zones are strongly affected by the post-glacial rebound (PGR) which has a specific signature in the observed gravity field.

This effect accounts for positive trends in these regions. According to the PGR models developed by [START_REF] Peltier | Global Glacial Isostasy and the Surface of the Ice-Age Earth: The ICE-5G(VM2) model and GRACE[END_REF] and [START_REF] Paulson | Inference of mantle viscosity from GRACE and relative sea level data[END_REF], its intensity can be greater in these to the redistribution of water masses. Unfortunately, trends computed on the 2 nd and 3 rd ICs, and their sum does not exclusively correspond to GIA, but to a mixture of geophysical remaining signals and noise. It is also worth noticing that the gravity signature of the Sumatra event in December 2004, which is clearly apparent on the Gaussian-filtered solutions (Fig. 8c and d), is not visible in the 1 st mode of the ICA solutions (Fig. 8a andb), but is present in the sum of the 2 nd and 3 rd ICs. This confirms that the ICA is able to isolate a pure hydrological mode in the GRACE products.

A second important difference between ICA solutions and Gaussian-filtered and destriped and smoothed solutions concerns the impact of the filtering radius on the trends estimate. An increase of the filtering radius causes a smoothing of the solutions, and, consequently, a decrease of the intensity of the trends on Gaussian-filtered and destriped and smoothed solutions (Fig. 8c to f). On the contrary, the intensity of the trends increases with radius of prefiltering on the ICA solutions. The location of the extrema is also shifted. This change in the location is a side-effect of the prefiltering with the Gaussian filter. A better location of the trends is observed when a Gaussian-filter of 400 km of radius is used instead of 500 km (see for instance the trends pattern in the Amazon and Orinoco Basins in Fig. 8a andb).

Basin scale comparisons

Changes in total water volume were estimated for 27 drainage basins whose locations are shown in Fig. 9 (the corresponding areas are given in Table 1). ICA400 and ICA500 were used to compute regional TWS averages versus time (Eq. 8 and 9). These times-series were compared to the G400 and G500, and the DS300 and DS500 respectively. Examples of the Amazon, the Ob and the Mekong basins for ICA400, G400 and DS300 (GFZ solutions) are presented in Fig. 10. TWS exhibit very similar temporal patterns for all the types of filtering radii and basins. The correlation coefficients between the ICA and the Gaussian-filtered, or the ICA and the destriped and smoothed time series of TWS are greater than 0.9 for 21 out of 27 basins. For five other basins (Amur, Colorado, Hwang Ho, Parana, St Lawrence), most of the correlation coefficients are greater than 0.8 or 0.9, and the others (generally the correlation coefficients between ICA and destriped and smoothed GFZ solutions) greater than 0.65 or 0.7. The only exception is the McKenzie Basin where all the correlation coefficients between ICA and destriped and smoothed are lower than 0.75 (r (ICA400,DS300) JPL =0.45 and r (ICA500,DS500) JPL =0.51). For the GFZ solutions, some unrealistic peaks are present in the destriped and smoothed TWS time series for some periods (Fig. 10a, 10b and 10c). These peaks only appear on the Gaussian-filtered solutions for the smallest basins, such as the one of the Mekong river (Fig. 10c), but are not present in the ICA-filtered solutions. RMS difference between ICAfiltered solutions and the other type of solutions are generally lower than 30 mm of equivalent water height and logically decrease with the radius of filtering. Differences with ICA solutions are generally larger for the GFZ-based destriped solutions, especially in tropical regions where the performances of the destriping are not the best and the hydrological signal the largest (Fig. 11).

We present an analysis of possible sources of error on the computation of regional averages versus time using the 300, 400, 500-km pre-filtered ICA solutions. This task is made on the longest available period of time for each center (CSR, JPL, GFZ), and for the 27 drainage basins (Table 1).

The formal error decreases with the number of points (see Eq. 10 and 11), the surface of the considered region, and the value of error k at each grid point. In the case of the Amazon Basin (~6 millions of km²), the formal error is only of 4.5 mm when k = 100 mm of equivalent-water height. For the Dniepr river basin that represents the smallest surface of our chosen basins (~0.52 million of km²), the formal error reaches 12.3 mm. For the series of basins, the maximum values of formal errors on the regional average are in the range of 10-15 km 3 of water volume per hundred of mm of error on the gridded points.

To estimate the frequency cut-off error , we made statistics of the numerical tests were performed to see what maximum error can be reached using Eq. 12. We computed this residual quantity for 300, 400 and 500 km-filtered ICA solutions and for each hydrological basin for N 1 =60 and N 2 =300. The maximum error is always less than 1 km 3 , as shown previously (Ramillien et al., 2006a), and it decreases with the filtering wavelength of the pre-processing. In other words, this error simply increases with the level of noise in the data. For the 300 km and 400 km pre-filtered solutions, the maximum values are found for the Amur River: 0.3 km 3 (CSR), 0.1 km 3 (JPL), and 0.8 km 3 (GFZ), and 0.04 km 3 (CSR), 0.09 km 3 (JPL), 0.08 km 3 (GFZ) respectively. While using the 500-km filtered solutions, the error of truncation is less than 0.005 km 3 .

The leakage error on the ICA solutions was computed per drainage basin (see section 3.5).

The results are presented in Fig. 12 for the different centres and radii of Gaussian prefiltering of 300, 400 and 500 km. We observed that the leakage decreases with radius of filtering and is generally lower in the JPL solutions, which presents lowest peak to peak amplitudes. This leakage error is logically greater in areas where several basins with large hydrological signal are close, i.e., South America with the Amazon, the Parana, the Orinoco and the Tocantins, or tropical Asia with the Ganges, the Brahmaputra and the Mekong rivers. The leakage error is also greater for basins with small hydrological signals close to a predominant basin having important hydrological variations, i.e., the Okavango and the Zambezi with the Congo. and NOAH outputs (SM). In the Murray-Darling River basin, GLDAS-NOAH simulations of SM range from 5 to 29% (in volumetric water content) across the basin for the study period, and are within typical values for monthly means at 1° resolution [START_REF] Lawrence | Soil moisture variability across climate zones[END_REF]. Besides, [START_REF] Leblanc | Basin scale, integrated observations of the early 21st century multiyear drought in southeast Australia[END_REF] show that, between 2003 and 2007, the linear rate of water changes for the GRACE TWS time series is similar to that observed for the annual total water storage from in situ observations and modeling. The combined annual anomalies of surface water, groundwater and soil moisture are highly correlated with the annual GRACE TWS (R = 0.94 and mean absolute difference =13 km 3 for the 2003 2007 period). These latter interannual field-based data are considered as the reference in the following analysis.

The results are presented on Fig. 14 for the GFZ solutions. The ICA solutions generally present a temporal pattern closer to the so-called reference (i.e., the sum of in situ data for SW and GW and model outputs for SM). This temporal pattern and the associated dynamics do not change with the radius of filtering, which is not the case with the Gaussian filtered, and the destriped and smoothed solutions. None of the solutions are able to retrieve the minimum observed in 2007 in the reference dataset. It is important to notice that the largest part of the interannual variations of TWS comes from the GW reservoir. For this hydrological component, the water storage is derived from in situ measurements through Eq. ( 1), and is highly dependent on the specific yield coefficient; where averages of local measurements are used to determine ranges of regional-scale estimates. For the Murray Darling Basin, composed of several aquifers, this leads to a broad range of spatial variability for the GW estimates. This variability is around or greater than 20 km 3 for 2002, 2007 and 2008 and around 10 km 3 for 2003, 5 km 3 in 2004 and 2006, 2.5 km 3 in 2005. We present yearly deviations to the reference for the different considered solutions (Fig. 15) over 2003-2006, where the range of variability of the GW is lower. Generally, the absolute deviation of ICAderived TWS to the reference is lower than 5 km 3 , especially for solutions filtered at 400 km.

At 500 km, the smoothing due to the preprocessing using a Gaussian filter is more important and explains the slight increase of the absolute deviation as the spatial resolution is degraded.

The destriping method and the Gaussian filtering also exhibit good performances, even if important deviations are sometimes observed (particularly with the Gaussian filter for a radius of 400 km).

Conclusion

The ICA-based approach is a very efficient method for successfully separating TWS from noise in the GRACE Level-2 data. We demonstrated that this method is more robust than classical filtering methods, such as the Gaussian filtering or the destriping. Comparisons at a global-scale showed that the ICA-based solutions present less north-south stripes than

Gaussian and destriped solutions on the land, and more realistic hydrological structures than the destriped solutions in the tropics. Trend maps over 2003-2008 have also been computed.

The corresponding trend maps present more realistic trend patterns than those obtained with other types of solutions (see for example over the Amazon and Orinoco Basins with the 400 km radius of prefiltering). ICA filtering seems to allow the separation of the GIA from the TWS as negative trends were found over the Laurentides and Scandinavia. Unfortunately, this important geophysical parameter does not appear clearly in an ICA mode yet. The major drawback of this approach is that it can not directly be applied to the GRACE Level-2 raw data, as a first step of prefiltering is required. In this study, we applied a Gaussian filtering which deteriorates the location of the important water mass patterns. This aspect of the pretreatment to be improved and highlights the necessity to replace the Gaussian filter used for preprocessing the GRACE Level-2 raw data by one more suited to improve the quality of the GRACE-derived TWS and thus obtain trustworthy estimate of the trends.

At the basin-scale, the ICA-based solutions allowed us to filter out the unrealistic peaks present in the time-series of TWS obtained using classical filtering for basins with areas lower than one million km². Among the ICA-based solutions, the JPL solutions are less affected by leakage compared with other solutions. JPL solutions also exhibit the lowest peak-to-peak amplitudes. The error balance of the GRACE-derived TWS is dominated by the effect of the leakage. Geofluids, 4, 187 196. Vigário, R. (1997). Extraction of ocular artifacts from EEG using independent component analysis, Encephalography and Clinical Neurophysiology, 103(3), 395-404. Tables Table 1: The 29 drainage basins considered in this study sorted by decreasing area. Figure 9: Location of the 27 drainage basins chosen in this study. See Table 1 for the correspondence between basins and numbers.

Figures

Figure 10: Time series of TWS (mm) derived from ICA400 (black), G400 (blue), DS300 (red) for GFZ solutions over the Amazon (a), Ob (b) and Mekong (c) basins.

Figure 11: RMS between ICA400 and G400 for CSR (dark blue), GFZ (red), and JPL (dark green) solutions and between ICA400 and DS300 for CSR (light blue), GFZ (orange), and JPL (green) solutions per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009. Tables Table 1: The 29 drainage basins considered in this study sorted by decreasing area. Figure 9: Location of the 27 drainage basins chosen in this study. See Table 1 for the correspondence between basins and numbers.

  original source signals, at least N observations are necessary if N sources are present. ICA remains applicable for square or over-determined problems. ICA proceeds by maximizing the statistical independence of the estimated components. As a condition of applicability of the method, non-Gaussianity of the input signals has to be checked. The central limit theorem is then used for measuring the statistical independence of the components. Classical algorithms for ICA use centering and whitening based on eigenvalue decomposition (EVD) and reduction of dimension as main processing steps. Whitening ensures that the input observations are equally treated before dimension reduction.

  several months (02/2003, 06 to 11/2004, 02/2005, 07/2005, 01/2006, 01/2007, 02/2009).

  for different radii of Gaussian filtering (300, 400 and 500 km) of GRACE mass anomalies. The kurtosis of the sum of the 2 nd and 3 rd ICs, corresponding to the north-south stripes, is most of the time, close to 0; that is to say that the meridian oriented spurious signals is almost Gaussian. Almost equal values of the kurtosis for the 1 st IC and the sum of the 2 nd and 3 rd ICs can be observed for several months. Most of the time, they correspond to time steps where the algorithm is unable to retrieve realistic TWS(02/2003, 08/2004, 11/2004, 02/2005, 01/2006, 01/2007, 02/2009).

  regions than the trends measured by GRACE. For example, in the Nelson Basin,Frappart et al. (in press) found a trend of TWS from GRACE of (4.5 0.2) mm/yr over 2003-2006, while model-based estimates of PGR represents 18.8 mm/yr in this region. As the space and time characteristics of the glacial isostatic adjustement (GIA) is different than the one from continental hydrology, the ICA approach may have separated it from the signals only related

  4.4. Basin scale validation GRACE observations were used to estimate variations in TWS over the Murray Darling Basin (see Fig. 13 for the location of the basin) at an interannual scale from January 2003 to December 2008. The annual variations in TWS from GRACE for different types of solutions were compared to the annual TWS computed as the sum of in situ observations (SW+GW)

Figure 1 :

 1 GRACE water storage from GFZ filtered with a Gaussian filter of 400 km of radius. (Top) First ICA component corresponding to land hydrology and ocean mass. (Bottom) Sum of the second and third components corresponding to the north south stripes. (a) March 2006, (b) September 2006, (c) March 2007, (d) March 2008. Units are millimeters of EWT.Figure 2: Time series of the kurtosis of the mass anomalies detected by GRACE after Gaussian filtering for radii of a) 300 km, b) 400 km, c) 500 km.

Figure 3 :

 3 Correlation maps over the period 2003Correlation maps over the period -2008 between the ICA-filtered TWS and the Gaussian-filtered TWS. Left column: ICA400-G400 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-G500 (b: CSR, d: GFZ, f: JPL).

Figure 4 :

 4 RMS maps over the period 2003RMS maps over the period -2008 between the ICA-filtered TWS and the Gaussian-filtered TWS. Left column: ICA400-G400 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-G500 (b: CSR, d: GFZ, f: JPL). Figure 5: Correlation maps over the period 2003-2008 between the ICA-filtered TWS and the destriped and smoothed TWS. Left column: ICA400-DS300 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-DS500 (b: CSR, d: GFZ, f: JPL).

Figure 6 :

 6 Figure6: RMS maps over the period[2003][2004][2005][2006][2007][2008] between the ICA-filtered TWS and the destriped and smoothed TWS. Left column: ICA400-DS300 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-DS500 (b: CSR, d: GFZ, f: JPL). Figure7: GRACE water storage from GFZ destriped and smoothed with a Gaussian filter of 300 km of radius for March 2006. (Top) First ICA component corresponding to land hydrology. (Bottom) Sum of the second and third components corresponding to the north south stripes. Figure8: Trend maps over the period 2003-2008 of TWS using the GFZ solutions a) ICA400, b) ICA500, c) G400, d) G500, e) DS300 and f) DS500. Figure9: Location of the 27 drainage basins chosen in this study. See Table1for the correspondence between basins and numbers. Figure10: Time series of TWS (mm) derived from ICA400 (black), G400 (blue), DS300 (red) for GFZ solutions over the Amazon (a), Ob (b) and Mekong (c) basins. Figure11: RMS between ICA400 and G400 for CSR (dark blue), GFZ (red), and JPL (dark green) solutions and between ICA400 and DS300 for CSR (light blue), GFZ (orange), and JPL (green) solutions per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009. Figure 12: Standard deviation of the leakage error (mm) per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009 for the ICA solutions.Figure 13: Map of the Murray Darling drainage basin in Australia. Cumulative rainfall deficit across the Murray Darling Basin for the 2001 2006 period and location of the shallow groundwater monitoring bores.Figure14: Comparison of GRACE TWS annual anomalies with hydrological estimates from in situ measurements (SW and GW) and modeling (SM) for the period 2003 2008. The grey curves correspond to in situ + model TWS, the blue to ICA-filtered solutions, the green to Gaussian filtered, and the red to destriped GFZ solutions at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. Figure15: Yearly deviation over 2003-2006 from in situ + model TWS of GRACE-derived TWS filtered with different approaches at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. In blue, the ICA solutions (dark blue: CSR, blue: GFZ, light blue: JPL), in green, the Gaussian solutions (dark green: CSR, green: GFZ, light green: JPL), and the destriped and smoothed solutions (red: CSR, orange: GFZ, yellow: JPL).

Figure 7 :

 7 Figure6: RMS maps over the period[2003][2004][2005][2006][2007][2008] between the ICA-filtered TWS and the destriped and smoothed TWS. Left column: ICA400-DS300 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-DS500 (b: CSR, d: GFZ, f: JPL). Figure7: GRACE water storage from GFZ destriped and smoothed with a Gaussian filter of 300 km of radius for March 2006. (Top) First ICA component corresponding to land hydrology. (Bottom) Sum of the second and third components corresponding to the north south stripes. Figure8: Trend maps over the period 2003-2008 of TWS using the GFZ solutions a) ICA400, b) ICA500, c) G400, d) G500, e) DS300 and f) DS500. Figure9: Location of the 27 drainage basins chosen in this study. See Table1for the correspondence between basins and numbers. Figure10: Time series of TWS (mm) derived from ICA400 (black), G400 (blue), DS300 (red) for GFZ solutions over the Amazon (a), Ob (b) and Mekong (c) basins. Figure11: RMS between ICA400 and G400 for CSR (dark blue), GFZ (red), and JPL (dark green) solutions and between ICA400 and DS300 for CSR (light blue), GFZ (orange), and JPL (green) solutions per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009. Figure 12: Standard deviation of the leakage error (mm) per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009 for the ICA solutions.Figure 13: Map of the Murray Darling drainage basin in Australia. Cumulative rainfall deficit across the Murray Darling Basin for the 2001 2006 period and location of the shallow groundwater monitoring bores.Figure14: Comparison of GRACE TWS annual anomalies with hydrological estimates from in situ measurements (SW and GW) and modeling (SM) for the period 2003 2008. The grey curves correspond to in situ + model TWS, the blue to ICA-filtered solutions, the green to Gaussian filtered, and the red to destriped GFZ solutions at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. Figure15: Yearly deviation over 2003-2006 from in situ + model TWS of GRACE-derived TWS filtered with different approaches at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. In blue, the ICA solutions (dark blue: CSR, blue: GFZ, light blue: JPL), in green, the Gaussian solutions (dark green: CSR, green: GFZ, light green: JPL), and the destriped and smoothed solutions (red: CSR, orange: GFZ, yellow: JPL).

Figure 12 :

 12 Figure6: RMS maps over the period[2003][2004][2005][2006][2007][2008] between the ICA-filtered TWS and the destriped and smoothed TWS. Left column: ICA400-DS300 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-DS500 (b: CSR, d: GFZ, f: JPL). Figure7: GRACE water storage from GFZ destriped and smoothed with a Gaussian filter of 300 km of radius for March 2006. (Top) First ICA component corresponding to land hydrology. (Bottom) Sum of the second and third components corresponding to the north south stripes. Figure8: Trend maps over the period 2003-2008 of TWS using the GFZ solutions a) ICA400, b) ICA500, c) G400, d) G500, e) DS300 and f) DS500. Figure9: Location of the 27 drainage basins chosen in this study. See Table1for the correspondence between basins and numbers. Figure10: Time series of TWS (mm) derived from ICA400 (black), G400 (blue), DS300 (red) for GFZ solutions over the Amazon (a), Ob (b) and Mekong (c) basins. Figure11: RMS between ICA400 and G400 for CSR (dark blue), GFZ (red), and JPL (dark green) solutions and between ICA400 and DS300 for CSR (light blue), GFZ (orange), and JPL (green) solutions per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009. Figure 12: Standard deviation of the leakage error (mm) per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009 for the ICA solutions.Figure 13: Map of the Murray Darling drainage basin in Australia. Cumulative rainfall deficit across the Murray Darling Basin for the 2001 2006 period and location of the shallow groundwater monitoring bores.Figure14: Comparison of GRACE TWS annual anomalies with hydrological estimates from in situ measurements (SW and GW) and modeling (SM) for the period 2003 2008. The grey curves correspond to in situ + model TWS, the blue to ICA-filtered solutions, the green to Gaussian filtered, and the red to destriped GFZ solutions at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. Figure15: Yearly deviation over 2003-2006 from in situ + model TWS of GRACE-derived TWS filtered with different approaches at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. In blue, the ICA solutions (dark blue: CSR, blue: GFZ, light blue: JPL), in green, the Gaussian solutions (dark green: CSR, green: GFZ, light green: JPL), and the destriped and smoothed solutions (red: CSR, orange: GFZ, yellow: JPL).

Figure 13 :

 13 Figure6: RMS maps over the period[2003][2004][2005][2006][2007][2008] between the ICA-filtered TWS and the destriped and smoothed TWS. Left column: ICA400-DS300 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-DS500 (b: CSR, d: GFZ, f: JPL). Figure7: GRACE water storage from GFZ destriped and smoothed with a Gaussian filter of 300 km of radius for March 2006. (Top) First ICA component corresponding to land hydrology. (Bottom) Sum of the second and third components corresponding to the north south stripes. Figure8: Trend maps over the period 2003-2008 of TWS using the GFZ solutions a) ICA400, b) ICA500, c) G400, d) G500, e) DS300 and f) DS500. Figure9: Location of the 27 drainage basins chosen in this study. See Table1for the correspondence between basins and numbers. Figure10: Time series of TWS (mm) derived from ICA400 (black), G400 (blue), DS300 (red) for GFZ solutions over the Amazon (a), Ob (b) and Mekong (c) basins. Figure11: RMS between ICA400 and G400 for CSR (dark blue), GFZ (red), and JPL (dark green) solutions and between ICA400 and DS300 for CSR (light blue), GFZ (orange), and JPL (green) solutions per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009. Figure 12: Standard deviation of the leakage error (mm) per basin (sorted by decreasing area of drainage basin) over the period October 2002 July 2009 for the ICA solutions.Figure 13: Map of the Murray Darling drainage basin in Australia. Cumulative rainfall deficit across the Murray Darling Basin for the 2001 2006 period and location of the shallow groundwater monitoring bores.Figure14: Comparison of GRACE TWS annual anomalies with hydrological estimates from in situ measurements (SW and GW) and modeling (SM) for the period 2003 2008. The grey curves correspond to in situ + model TWS, the blue to ICA-filtered solutions, the green to Gaussian filtered, and the red to destriped GFZ solutions at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. Figure15: Yearly deviation over 2003-2006 from in situ + model TWS of GRACE-derived TWS filtered with different approaches at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km. In blue, the ICA solutions (dark blue: CSR, blue: GFZ, light blue: JPL), in green, the Gaussian solutions (dark green: CSR, green: GFZ, light green: JPL), and the destriped and smoothed solutions (red: CSR, orange: GFZ, yellow: JPL).

  Figures Figure 1: GRACE water storage from GFZ filtered with a Gaussian filter of 400 km of radius for March 2006. (Top) First ICA component corresponding to land hydrology and ocean mass. (Bottom) Sum of the second and third components corresponding to the north south stripes. (a) March 2006, (b) September 2006, (c) March 2007, (d) March 2008. Units are millimeters of EWT. March 2006 September 2006

Figure 2 :

 2 Figures Figure 1: GRACE water storage from GFZ filtered with a Gaussian filter of 400 km of radius for March 2006. (Top) First ICA component corresponding to land hydrology and ocean mass. (Bottom) Sum of the second and third components corresponding to the north south stripes. (a) March 2006, (b) September 2006, (c) March 2007, (d) March 2008. Units are millimeters of EWT. March 2006 September 2006

Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 34567 Figure 3: Correlation maps over the period 2003-2008 between the ICA-filtered TWS and the Gaussian-filtered TWS. Left column: ICA400-G400 (a: CSR, c: GFZ, e: JPL). Right column: ICA500-G500 (b: CSR, d: GFZ, f: JPL).

Figure 8 :

 8 Figure 8: Trend maps over the period 2003-2008 of TWS using the GFZ solutions a) ICA400, b) ICA500, c) G400, d) G500, e) DS300 and f) DS500.

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: Time series of TWS (mm) derived from ICA400 (black), G400 (blue), DS300 (red) for GFZ solutions over the Amazon (a), Ob (b) and Mekong (c) basins.

Figure 13 :

 13 Figure 13: Map of the Murray Darling drainage basin in Australia. Cumulative rainfall deficit across the Murray Darling Basin for the 2001 2006 period and location of the shallow groundwater monitoring bores.

Figure 14 :Figure 15 :

 1415 Figure14: Comparison of GRACE TWS annual anomalies with hydrological estimates from in situ measurements (SW and GW) and modeling (SM) for the period 2003 2008. The grey curves correspond to in situ + model TWS, the blue to ICA-filtered solutions, the green to Gaussian filtered, and the red to destriped GFZ solutions at a) 400 km of filtering (300 km for the destriped solutions) and b) 500 km.

  Validation with in situ measurements was performed in the Murray Darling Basin where broad networks of in situ measurements of SW and GW are available. The ICA-based solutions are in better agreement with in situ data compared with the other types of solutions, especially where prefiltering at a 400 km of radius. The maximum deviations are lower by a factor two or three compared with the other filtering methods. Acknowlegdements[START_REF] Strassberg | Comparison of seasonal terrestrial water storage variations from GRACE with groundwater-level measurements from the High Plains Aquifer (USA)[END_REF] Comparison of seasonal terrestrial water storage variations from GRACE with groundwater-level measurements from the High Plains Aquifer (USA). Geophysical ResearchLetters, 34, L14402, doi:10.1029/2007GL030139.

	Swenson, S., & Wahr, J. (2006). Post-processing removal of correlated errors in GRACE
	data, Geophysycal Research Letters, 33, L08402, doi:10.1029/2005GL025285.
	Syed, T.H., Famiglietti, J.S., Chambers, D. (2009). GRACE-based estimates of terrestrial
	freshwater discharge from basin to continental scales, Journal of Hydrometeorology, 10(1),
	doi: 10.1175/2008JHM993.1.
	Tapley, B.D., Bettadpur, S., Ries, J.C., Thompson, P.F. & Watkins M. (2004). GRACE
	measurements of mass variability in the Earth system. Science, 305, 503-505.
	Tregonning, P., Ramillien, G., McQueen, H. & Zwartz, D. (2009). Glacial isostatic
	adjustement observed by GRACE. Journal of Geophysical Research, 114, B06406,
	doi:10.1029/2008JB006161.
	Urbano, L. D., Person, M., Kelts, K. & Hanor, J. S. (2004). Transient groundwater impacts on
	the development of paleoclimatic lake records in semi-arid environments.

This work was partly supported by the fondation Sciences et Techniques pour (STAE) in the framework of the CYMENT Project (post-doctoral grant for Frédéric Frappart) and the Australian Research Council grant DP0666300 (Marc Leblanc and Sarah Tweed). The authors also wish to thank Dr. Matthew Rodell for his help during the revision process.


Evrot Icssea 2010 -Pétin

Jean-François Pétin

Dominique Evrot

Gérard Morel

Pascal Lamy

Combining SysML and formal models for safety requirements verification

Keywords: System specification, SysML, formal verification, Discrete Event Systems, Safety requirements

Industrial safety-related standards strongly recommend the use of formal methods to control the complexity of software-intensive automation. This paper deals with the verification of safety requirements for the design of an industrial machinery control system embedding safety-critical software applications. It combines SysML semi-formal modelling approach to capture and structure safety requirements and model-checking techniques for the formal verification purposes.

INTRODUCTION

This paper deals with the verification of safety requirements for the design of complex systems involving software, mechanical, electrical, or pneumatic components for industrial safety-critical machinery. The control of such complex systems results from an interaction network between all the system components that can introduce undesired behaviour with critical impact on safety. Consequently, safety properties of the control system cannot be verified by only proving local properties of each component, but need to be studied through the emerging behaviour that issues from the interaction network. The main difficulties of this last verification are caused by the heterogeneity of the control components and technology domains that use proper design formalisms and tools.

In the development of industrial automation, standards, such as the safety-related IEC 61508 1 or its application in the machinery domain (IEC 62061 2 ), strongly recommend the use of formal methods to control the complexity of softwareintensive applications. In this way, the role of INRS (Institut National de Recherche et de Sécurité), acting as a French independent institute, is to anticipate future risk by acquiring new knowledge and converting current knowledge into practical know-how available to professional in charge of prevention. This paper presents a combination of SysML semi-formal modelling approach, to identify, specify and refine properties and architectures of an INRS case study, with formal models for the design of its control system and for the verification of its dynamic properties. Section 2 introduces the problem statement by showing complementarities between systemoriented approaches and formal models. This section ends by the description of INRS case study which is developed using our proposed methodology in section 3. The final section provides some conclusions and open issues.

PROBLEM STATEMENT

Automation engineering has to ensure that the system under control to be developed, noted S, is compliant with the enduser requirements, noted R. Taking into account the fact that:

• requirements R are usually captured in multiple levels of abstraction and broken down into sub-requirements at different levels of abstraction, noted R = {R 1 , R 2 , … R n }, • the system under control is often composed of heterogeneous components (mechanical, software, electronic subsystems) that cooperate to achieve the systems goal, noted S = {S 1 , S 2 , …, S m }, the initial problem becomes a non bijective relationship between the set of requirements and the set of system components. Fusaoka [START_REF] Fusaoka | A description and reasoning of plant controllers in temporal logic[END_REF] and Supervisory Control Theory [START_REF] Ramadge | Supervisory control of a class of discrete event processes[END_REF] have formalised this problem focusing on Discrete Event System dynamics point of view. This formulation can be extended [START_REF] Pétin | Formal specification method for production systems automation[END_REF] by taking into account system features (dynamics but also functions, structure, information, etc) according to the predicate S S Goal (where S c means Control system and S p means physical system) that has to be preserved through the automation engineering processes. More generally speaking, it means that:

• a given system component S i may satisfy a sub-set of requirements such that S R , , • a given requirement R j may be satisfied by a subset of system components such that S , R . Identifying and proving the preservation of these relationships for the safety requirements, all along the system engineering process, is our main objective. Two ways of thinking can be found in the scientific literature to partially cover this issue: formal methods & models and system engineering approaches.

Formal models

Formal modelling process relies on mathematical foundations that may ensure that a given system design and/or implementation fulfils a given model of requirements. Conceptual and practical approaches have been widely explored by academic bodies related to computer sciences and automatic control: examples are software verification by theorem proving [START_REF] Abrial | The B-book: Assigning programs to meanings[END_REF], model checking [START_REF] Clarke | Model checking[END_REF], or automatic synthesis [START_REF] Cassandras | Introduction to Discrete Event Systems[END_REF]. More precisely in the area of automation, these approaches are applied:

• during the design phase, to prove control model properties, using analysis techniques and model checking based on the building of the reachable states space of the control software, • during the implementation phase, to prove PLC (Programmable Logic Controller) programs properties by applying formal techniques to check the properties satisfied by controllers implemented using IEC 61131-3 programming languages [START_REF] Roussel | An algebraic approach for PLC programs verification[END_REF].

Due to the combinatory explosion phenomena induced by the exploration of the whole state space, these techniques are efficient for small-sized components. For larger scale applications, system global requirements must be transformed into multi-component local properties that can be processed by a model-checker. Bridging the gap from global system properties to a composition of local properties remains difficult [START_REF] Barragan | Obtaining temporal and timed properties of logic controllers from fault tree analysis[END_REF].

Second difficulty is linked to the underlying mathematical representation which is limited to the modelling of system dynamics, and hardly covers the description of other system features. More precisely, the decomposition of a complex system in several sub-systems, involving software, hardware and man-machine interface components, does not facilitate the traceability and the refinement of global system properties through the sub-system decomposition, design, verification and validation. This aspect requires method that enables model refinement to capture very abstract requirements including safety properties and to project them into all system components [START_REF] Evrot | Formal specification of safe manufacturing machines using the B method: application to a mechanical press[END_REF]. Following the consensus that early phases of a system definition are the most important in ensuring that the target system will satisfy the user's requirements, many systems and automation engineering practitioners [START_REF] Pétin | Formal specification method for production systems automation[END_REF] [10] [START_REF] Johnson | Improving automation software dependability: a role for formal methods ?[END_REF] consider that time is ripe to formalize the earlier phases of system specification. This challenge justifies the use of system-oriented approaches [START_REF] Sheard | Complex Systems Science and its Effects on Systems Engineering[END_REF] in automation engineering.

System-oriented models

Standards define system engineering as an "interdisciplinary collaborative approach to derive, evolve and verify a life cycle balanced system solution which satisfies customer expectations and meets public acceptability" (IEEE 1220 3 , ISO 15288 4 ). System engineering is usually involving two complementary and inter-correlated processes related to system definition that aims to identify and analyse the system needs and requirements and to system development that focus on the system design in accordance with identified requirements. This definition meets our goal that is to identify safety system requirements issued from an interdisciplinary approach and to design and verify a system control that contributes to satisfy the identified requirements. Models and tools promoted by system engineering are most of the time based on a wide set of graphical formalism that help to capture the various system features as well for its definition as for its development (requirements engineering, functional, structural and behavioural modelling, etc). Some of them, based on the Unified Modelling Language (UML 5 ) and its extensions, cover the various aspect of system modelling, while other ones focus on a specific part such as requirements expression (DOORS 6 ) or documentation traceability (Reqtify7 ). UML and its extension UML2.0, is a system modelling toolbox that should address use cases identification and architectural, behavioural design of the system. These tools allow a multi-trades approach due to their high level of abstraction and its graphic notation. However, it seems that its origin and primary application being software-oriented, this language cannot be adequate for all aspects of systems engineering (hardware, software, man-machine interface…). The new profile SysML has been developed to meet system engineering community expectations. SysML toolbox contains new diagrams, like the requirement diagram, that allows a most efficient system analysis and design. Using UML/SysML based models for automation engineering has been identified as an efficient way to handle the complexity of current control applications [START_REF] Bonfé | Object-Oriented modeling of logic control systems for industrial applications[END_REF]. OOONEIDA [START_REF] Auinger | Requirements and solutions to Software encapsulation and engineering in next generation manufacturing systems: OOONEIDA approach[END_REF], TORERO [START_REF] Ferrarini | Control functions development for distributed automation systems using the Torero approach[END_REF] and CORFU projects [START_REF] Tranoris | A tool supported engineering process for developing control applications[END_REF] as well as the definition of UML profile for process automation [START_REF] Katzke | UML-PA as an Engineering Model for Distributed Process Automation[END_REF] can be mentioned in this way. However, the wide notation toolbox provided by UML/SysML suffers from a lack of formal semantics for validation and verification issues. Efforts have been made towards UML formalization but they are often focused on static equivalence between diagrams in a Model Driven Engineering perspective: executive semantics of the dynamics diagrams are formalized by providing equivalence to DES models such as Petri Nets.

Our objective is to bridge the gap between a necessary semi-formal approach, such as SysML and formals model that are efficient for validation and verification of safety properties. In this way, an integrated method has been proposed to identify the safety requirements of a complex control system using the SysML requirements diagram, to allocate them on a given set of system components, to define the impact of this allocation in terms of safety properties to be satisfied by each system component, and then to provide the model-checkers with this correct model of component properties.

3/10

Case study

We consider a mechanical press that aims at stamping metal products with a tool in vertical movement. A crankshaft equipped with two sensors for top and bottom dead centres gives this vertical movement. A clutch, which is actuated by a pneumatic valve, ensures the transmission between motor and crankshaft. Operator protection is supplied thanks to a two hands command. An emergency stop enable or not the stamping action. We focus on the press function that is in charge of clutching and braking the press motor in the normal operating mode (motor is considered as always working).

In this mode, the press is initially stopped in up position and waits for a two hands command. Then, the tool is falling until the bottom dead centre is reached, then is rising back, and finally stopped when in up position. Moreover, if the operator releases the hand command during the falling phase, or if emergency stop occurs, the press is stopped.

Fig. 1. Mechanical press case study

Usually, control of such safety control system is performed by wired technology. Due to economic constraints (control application reusability, functions complexity), control of industrial machinery ensures more and more safety functions using Safe Programmable Logic Controller (Safe-PLC). Integration of safety functions into the software parts of programmable systems can be envisaged only if the safety level of programmable technology set is at least equal to the wired technology. One role of INRS is to study methods for the development of safety machines according to recommendations listed in the annex IV of the machinery directive [18] and IEC 61508. In particular, the higher level SIL 4 recommends the use of formal methods to control the quality of software intensive applications.

In this way, INRS has studied the B method [START_REF] Abrial | The B-book: Assigning programs to meanings[END_REF] and ordered a case study development using this method [START_REF] Abrial | Case study of a complete reactive system in Event-B: a mechanical press controller, tutorial 3 of ZB[END_REF]. This study has clearly shown the benefit of using B for proved code generation. Nevertheless, modelling methodology is not explicit and relies on the expert practice. In the context of formal system engineering and taking into account the need for traceability of safety properties, this approach must be enrich by modelling guidelines:

• to better detail how the safety requirements are taken into account when modelling,

• to enlarge software engineering to automation and system engineering including software, hardware, man-machine interface and physical systems. In other words, starting the study of a complex system using formal models appears to be hardly achievable in practice. It justifies the use of less formal model such as SysML [START_REF] Evrot | Using SysML for identification and refinement of machinery safety properties[END_REF] to cover the first phases of specification using a graphical language that promotes the knowledge sharing between the involved engineers.

SYSML-BASED METHODOLOGY TOWARDS REQUIREMENT VERIFICATION

Requirements are the formalization of the user needs for one part and the formalization of provider limits or standards prescriptions for other ones. Requirements can be expressed at a high level of abstraction when dealing with system global features but can also be expressed at a very low level of abstraction when depending on technical local choices for system component (for example material and immaterial barriers to be used). Structuring this set of requirements is based on an incremental modelling, in the same idea than the B formal refinement, which consists in substituting an abstract requirement by a set of more concrete ones. Using SysML Requirements diagram, requirements are modelled with a class stereotype including an open list of attributes such as text definition, source, id, and so on. Structuring the requirements is done using the composition link provided by SysML, which means that the composed requirement is realized if and only if all the components requirements are realized. Figure 2 illustrates the decomposition of one abstract requirement into several more concrete ones. Benefit of such a diagram is to facilitate the elicitation and structuring of requirements that are usually expressed by the end-users and/or the standards in natural language [START_REF] Blaise | Formalisation of normative knowledge for safe design[END_REF]. This set of requirements must then be allocated (or mapped) onto system functions and/or components. Functions are modelled using classical Activity and/or Use Case diagrams. System components are modelled using SysML Block defined as class stereotypes. Structure of the control system is modelled using Block Definition Diagram (BDD) while detailed definition of control components can be modelled using Internal Block Diagram (IBD). IBD is an internal description of a component, similar to synchronous DES models, and which can be exported to widely used tools for control design (Simulink 8 , Scade 9 , ControlBuild 10 , etc). 8 Simulink is a MathWorks, Inc. product Allocation is made using the following SysML links: satisfy for the mapping between requirements and blocks (figure 3) and refines for the mapping between requirements and functions (figure 4). Block and activity relationship is modelled with allocate link: it means that an activity having a refines link with a requirement must be implemented by the block that satisfy this given requirement (figure 4). Note that requirement allocation is preserved by the block composition link. It means that if a requirement is satisfied by a component, then it is also satisfied by all its sub-components. Property is a requirement stereotype that formalizes requirements textual description in terms of logic predicates and/or numerical parameters (figure 4a). It proposes a list of attributes involving language type, regular expression, kind, and state machine describing some sequential properties. This extension enables:

• to make the link between system modelling and trade-oriented design with the objective to prove that behaviour of a given component is compliant with local expected properties; in this way, concrete properties Pi of each components can be reused as logic or temporal predicates by model checking tools, such as UPPAAL 11 or can be reused as post-conditions by simulation and testing tools, • to formally demonstrate that a logic combination of component's concrete properties Pc(j) establishes the system abstract properties Pa(i), using a theorem prover, such as COQ 12 (in figure 4a, P2 is composed of P2.1 and P2.2). 9 SCADE is a Lustre-based Esterel Technologies product 10 ControlBuild is IEC61131-3-based Geensys product 11 UPPAAL is developed by the Information Technology Department of the Uppsala Univ. (Sweden) and the Computer Science Department of the Aalborg Univ. (Denmark). 12 Huet G., Khan G. and Paulin-Mohring C, 1995. The Coq proof assistant, Technical Report 178, Inria. Integration of tools for system specification with SysML, for control system design & simulation using trade-oriented languages (ControlBuild) and for proving safety properties (UPPAAL model checker), have been impleneted using XML and XSLT transformations based on a meta-model of requirements, properties and blocks shared concepts (Figure 4a). Following sub-sections illustrates all the steps of the proposed methodology using the INRS case study.

APPLICATION ON THE INRS CASE STUDY

Safety requirements modelling

Safety requirements for the mechanical press are issued from the ISO 12100-1 risk analysis including:

• a hazard identification and risk assesment based on risk analysis methods like fault tree analysis or Failures Modes Effects Criticality Analysis (FMECA). • a risk reduction based on prevention measures, integration of safety devices and information notices. Safety analysis shows that the most important risk is a crash risk of operator's hands during the stamping movement. It identifies two access types to the dangerous area, a front access and two lateral accesses. Because of the crankshaft position in the press architecture, lateral access are potentially damageable not only during the downing phase (when the tool stamp the iron sheet) but also during the rising phase (when the tool return to its initial position). So lateral access must be avoided during all the press movement (requirement SR3) and front access must only be avoided during the downing phase (requirement SR1). According to the ISO 12100-2 protector choice algorithm, we choose removable protectors for the lateral access (SR3) and a two hands command (THC) for the front access (SR1). THC role is to require that the operator use his two hands to press on the THC two devices (kind of pushing buttons) if he wants to activate the stamping movement. All protection devices are interlocked with the press stamping movement. It means that if a protection is not activated during the stamping movement then the movement is stopped (Requirement SR2). These first set of abstract requirements are modelled in Figure 5a. Our study focuses onto the two hands command. Requirement SR1 is split according to EN954-1 general standard and to EN 574 two hands commands dedicated standard into a set of five detailed requirements related to:

• ergonomics of the two hands command (requirement SR1.2 prevent that THC can be activated with only one hand), • self fault detection (requirement SR1.4 says that THC must be desactivated if a fault is detected),

• two hands command activation and deactivation (requirements SR1.1 means that THC is reactivated only if both devices have been released, SR1.3. means that THC is activated only if both devices have been actuated in a 0,5 s interval and SR1.5 means that THC is deactivated if one device is released. These new sub-requirements are modelled by enriching the previous SysML Requirement Diagram and by associating semi-formal predicates associated to each requirement (figure 5b) using property stereotype. It means that satisfying property P1 is equivalent to fully satisfy properties P1.1. to P1.5. Proving this equivalence between first order logic predicates can be done manually or using a theorem prover such as COQ. Note that a same approach must be used to refine requirement SR2 and SR3. Next phase consist identifying the system architecture including software and hardware components in order to allocate each requirement on a supporting device.

Safety requirements allocation

Requirements allocation using "satisfy" link enable to link one or more blocks (system component) to one or more requirements. Two hands command requirements allocation is given in table 1. Description of structural and functional views of the press-system is done using respectively:

Requirements

• block diagrams that define the system component architecture. For example classical THC has been architecture involve devices (pushing buttons), a cover, two switches for faults detections, and software component (Figure 6a). • activity diagrams that define operation calls and data exchange. For example, Figure 6b Translation of sequential properties is more difficult. It requires the use of an observer; i.e. a state machine that registers the system evolution and that introduces a faulty state in which the system must never go. For example, P1.4a is represented by such an observer which involves 3 states (figure 10): the first one with THC active, the second one with THC inactive and both devices released and a third one, halfway, with THC inactive and only one device released. The property P1.4a formalizes that activating THC is only possible if both devices have been released since the previous activation. So when state 3 is active (THC have just been deactivated), it is possible to reach state 1 only by activating state 2. Consequently, the faulty state, corresponding to the evolution that must be avoided, can be reached from state 3 if THC is activated. Verification of the properties is done by proving using UPPAAL that the faulty state can never be reached.

Discussion

In this study, only safety properties have been focused on, the functional properties have not been considered. This is justified by the chosen implementation that follows the recommendations provided by the standards related to manufacturing machines, i.e. to separate safety related components from nominal components. Consequently, the chosen implementation in Safe-PLC is based on a safety filter that is placed between control system and physical system and that THC_react THC_desact THC_fault THC 9/10 enables or disables the control system outputs [START_REF] Marangé | Proposition of plant model for the verification of system safety[END_REF]. In other words, it means that this filter ensures that the safety requirements are always satisfied whatever the outputs of the control system are (figure 11).

Fig 11 . Safety filter [START_REF] Marangé | Proposition of plant model for the verification of system safety[END_REF] Benefit is that functional components of the machine control can be developed according classical approaches while more costly methods needed to reach level 2, 3 or 4 of SIL (Safety Integrity Level) are limited to safety components that are involved in safety filter. Note that splitting functional and safety requirements is not always an obvious task: if a component realizes both safety and nominal functions then it must be consider entirely as safety related.

CONCLUSION

There is a growing interest in formal methods and tools that facilitate the validation of software-intensive automation systems. This interest becomes a legal requirement when dealing with safety-critical systems; the IEC 61508 safetyrelated standard strongly recommends the use of formal verification methods to be applied in the certification process by the suppliers, integrators, or independent external authorities, but without defining how they can be applied. On the other way, system engineering semi-formal approaches appears to be essential for capturing and structuring the requirements of a complex system and for proposing high level system functional architectures that cope with the identified requirements.

Trying to match theses apparently opposite way of thinking, we have proposed a system modelling approach that combines non-formal methods based on SysML requirements and block diagrams, and formal methods such as modelchecking to prove that the local behaviour of each system component contributes to satisfy system requirements. This work is based on the mapping between structured models of the requirements to be satisfied, the functions to be realized and the component to be implemented. One limit must be mentioned: composition of properties is assumed to be described by a simple logical expression (i.e. satisfying a property is assumed to be equivalent to satisfy a logical combination of sub-properties) but it is not always true especially when manipulating temporal expression involving operators such as before, after or until. Further work is currently developed in this way. Benefit relies on the requirements traceability all along the system engineering process and on the associated verification capabilities. This point is very helpful for the development of COTS-based (commercial off-the-shelf) control, where subcontractors' requirements must be clearly identified.

Although these interdisciplinary exchanges between computer science and system engineering approaches demonstrate that they contribute to verify the safety integrity highest levels, common experiments on laboratory-scale and industrialscale case-studies emphasize that effort must still be employed to make the proposed engineering framework effective in practice.

  Fig 2. Requirements refinement

Fig 3 .

 3 Fig 3. Requirements allocationIn order to provide proofs about requirements composition and allocation, we need a more formal description of the different safety properties, associated to each abstract and concrete requirement, respectively noted Pa(i) and Pc(j). In this way, we have extended the SysML meta-model to introduce property concept.

Fig 9 .

 9 Fig 7. THC requirements allocation

Fig 10 .

 10 Fig 10. Observer for property P1.4aAll these translations have been implemented using model transformation as recommended by Model Driven Engineering. It is based on the definition of a XML common format using XSLT transformations to enable the exchange of data and models between different commercial tools, such as those used for the INRS case study: SysML system specification (Magic Draw), control system design and simulation using trade-oriented languages (ControlBuild), model-checker for proving the safety properties (UPPAAL).

Table 1 .

 1 Requirements allocation

		Components
	SR 1	All
	SR 1.1	Software / Device
	SR 1.2	Cover / Switch
	SR 1.3	Software / Device
	SR 1.4	Software / Device
	SR 1.5	Software / Device

IEC 61508, Functional safety of electrical/electronic/ programmable electronic (E/E/PE) safety-related systems

IEC 62061, Safety of machinery-Functional safety of electrical, electronic and programmable electronic control systems

Standard for application and management of systems engineering

Systems and software engineering, System life cycle processes

OMG, Object Management Group, www.omg.org

Telelogic, www.telelogic.com

Geensys, www.geensys.com

Acknowledgment

This work has been partially funded by INRS in the framework of Dominique Evrot PhD thesis (Contribution à la vérification d'exigences de sécurité: application au domaine de la machine industrielle, Thèse de l'Université H. Poincaré, juillet 2008, in French). We would address our acknowledgement to our INRS partners who participated to this project.

Properties associated to requirements become the properties that each component must satisfy. This operation can lead:

• to decompose a property into several sub-properties when a requirement and its associated property is allocated to several component; this is the case for properties 1.1, 1.3, 1.5 which are projected on software component (properties 1.1.a, 1.3.a, 1.5.a) and on device component (1.1.b, 1.3.b) and for 1.4 property which is split into software (1.4a), device (1.3b/1.4c) and switch blocks (1.4b). • to possibly merge properties when a same component satisfies several properties; this is the case of the device block which satisfies part of 1.3 and 1.4 properties that can be merged into 1.3b/1.4c property. Merging operation between predicates results from a logical AND between elementary predicates. Correctness of these decomposition and merging operations can be established and demonstrated using manual or a theorem prover. It leads to the final THC requirements allocation given by the figure 11 where LD and RD represent respectively the activation of left and right devices or pushing buttons, C_LD and C_RD represent respectively the activation of fault switch associated to left and right pushing devices while HTC represent the activation of two hand command. Next phase consists in a detailed design of the identified components (and more particularly the behavioural description of the system), and in proving that the system design is compliant with the identified properties.

Safety requirements verification

The behavioural design of the THC software control part has to satisfy the properties 1.1a, 1.3a, 1.4a and 1.5a identified in figure 8. This design must also be compliant with the functional architecture shown in figure 8. It can be done using internal block diagrams (IBD) and state diagrams provided by SysML. Note that the synchronous data flow structure provided by IBD facilitates its transcription into dedicated control design tool based on synchronous DES models such as Simulink, Scade or ControlBuild (figure 9).