Regional GRACE-based estimates of water mass variations over Australia : validation and interpretation

Time series of regional 2 ◦ × 2 Gravity Recovery and Climate Experiment (GRACE) solutions have been computed from 2003 to 2011 with a 10-day resolution by using an energy integral method over Australia (112 ◦ E–156 E; 44 S–10 S). This approach uses the dynamical orbit analysis of GRACE Level 1 measurements, and specially accurate along-track K-band range rate (KBRR) residuals with a 1 μm s−1 level of errors, to estimate the total water mass over continental regions. The advantages of regional solutions are a significant reduction of GRACE aliasing errors (i.e. north– south stripes) providing a more accurate estimation of water mass balance for hydrological applications. In this paper, the validation of these regional solutions over Australia is presented, as well as their ability to describe water mass change as a response of climate forcings such as El Niño. Principal component analysis of GRACE-derived total water storage (TWS) maps shows spatial and temporal patterns that are consistent with independent data sets (e.g. rainfall, climate index and in situ observations). Regional TWS maps show higher spatial correlations with in situ water table measurements over Murray–Darling drainage basin (80–90 %), and they offer a better localization of hydrological structures than classical GRACE global solutions (i.e. Level 2 Groupe de Recherche en Géodésie Spatiale (GRGS)) products and 400 km independent component analysis solutions as a linear combination of GRACE solutions provided by different centers.


Introduction
Since its launch in 2002, the Gravity Recovery and Climate Experiment (GRACE) space mission, has provided first-time global estimates of changes in total water storage (surface water, soil moisture and groundwater) with unprecedented accuracy at global scales.GRACE data have already demonstrated a strong potential for estimating hydrological information such as discharges (Syed et al., 2008), evapotranspiration (Rodell et al., 2004;Ramillien et al., 2006a), groundwater (Rodell et al., 2007;Leblanc et al., 2009;Taylor et al., 2013) and extreme climate events as floods and droughts (Andersen et al., 2005;Frappart et al., 2012Frappart et al., , 2013a;;Houborg et al., 2012).Several analysis centres, including the Centre for Space Research, University of Texas at Austin (CSR), Jet Propulsion Laboratory (JPL), Ge-oForschungsZentrum (GFZ) and Groupe de Recherche en Géodésie Spatiale (GRGS), use GRACE Level 1B observations to produce time series of global gravity solutions developed in spherical harmonics.Static gravity field and its timevarying variations (i.e.atmospheric and oceans mass changes including the periodic tides) are removed during the GRACE orbits pre-treatment.However, these correcting models remain imperfect by their incompleteness in the description of natural phenomena by omission and/or by the lack of resolution.As GRACE solutions are averages over constant time intervals, errors in the correcting models with periods from hours to days contaminate 10-day and monthly GRACE solutions by aliasing, and thus degrades their accuracy (Duan et al., 2012).This effect of signal distortion deteriorates the Published by Copernicus Publications on behalf of the European Geosciences Union.
L. Seoane et al.: Regional GRACE-derived water mass variations over Australia detection quality of true water mass signals into other time frequencies, and makes these signals indistinguishable once they are sampled.In the case of GRACE, aliasing is complicated since mass variations occur both in time and space, but they are sampled mainly along satellite tracks in the nearly latitudinal direction, and represented through global spherical harmonic (SH) functions afterwards.The particular distribution of the GRACE tracks due to the single polar orbit geometry translates into north-south "stripes" on the 10-day and monthly gravity field solutions, by the fact that the ability of such an orbit configuration for recovering short-term mass and small variations is limited.The determination of the GRACE-Stokes coefficients leads to underdetermined systems of normal equations to be solved, by creating correlations between SH of high degrees (i.e.> 10-15) (see Swenson and Wahr, 2006), and amplification of the orbit error and the data noise (Himanshu et al., 2012).Another serious drawback of using SH functions for representing the gravity field is "leakage", as energetic signals propagate on the entire terrestrial sphere, since SH are defined as global undulations, and thus they come to pollute the water mass estimates in the geographical region of interest.This is particularly the case for a small region (i.e.typically about 1 millions of square km), whereas averaging over a large area, such as the Australian continent, surely cancels the spurious SH undulations from outside.Filtering techniques can partly remove some of these effects (Swenson et al., 2006;Klees et al., 2008;Frappart et al., 2011).New regional GRACE solutions using an energy integral method (Ramillien et al., 2011(Ramillien et al., , 2012) ) yield to significant gains in the accuracy and resolution of total water storage estimates compared to classical strategies using the spherical harmonics.The present study aims to validate these new regional solutions over Australia, the driest inhabited part of the world.On average, the seasonal amplitude of total water storage (TWS) in Australia is 3 times lower than the one observed in the Amazon.Thus, the challenge for us is to detect realistic structures while noise amplitude in GRACE data is closer to the amplitude of the true water mass variations over this continent.We propose that water storage variations from regional solutions over Australia are compared with independent information given by the in situ observations (rainfall, surface water storage groundwater levels).GRACE products have already been used to study time variations of water mass storages and transfers over Australia.In addition, the GRACE-based water storage solutions have been used to evaluate the accuracy of model estimates (Van Dijk et al., 2011;Forootan et al., 2012).Syed et al. (2008) used GRACE data to estimate a depletion of 1.3 mm month −1 of total water storage in Australia.Leblanc et al. (2009) used GRACE global solutions to analyze the impact of the 2007 Millennium Drought on the water resources in the Murray-Darling Basin.Awange et al. (2009) noted the typically small hydrological signals in Australia that are not detectable because of errors in the standard GRACE data processing (e.g.oceanic non-tidal model deficiencies) and min Sy mean 0.70, 13.17 0.67, 14.20 0.48, 16.30 Sy max 0.68, 14.63 0.65, 15.53 0.43, 18.44  isotropic filtering methods.Tregoning et al. (2008) showed the propagation of error of non-tidal oceanic modelling in the Gulf of Carpentaria onto the continent in the northern pat of Queensland and Northern Territory states.Frappart et al. (2011) proposed a filtering technique based on independent component analysis (ICA) that reproduced better the in situ observations over the Murray-Darling Basin than the Level-2 GRACE GFZ, CSR and JPL solutions when these data sets are used individually.
Our regional approach, based on the energy integral method, can potentially improve further the accuracy of GRACE TWS estimations.Over South America, Ramillien et al. (2012) have shown that this regional approach offers a reduction of both north-south striping due to the distribution of GRACE satellite tracks, and temporal aliasing of correcting models.Besides, regional solutions present more realistic spatial and temporal patterns than the global ones when compared to in situ data sets (Frappart et al., 2013b).
In the following section, all data sets are presented.To make the comparison between different data sets, we separate the GRACE and models signals into main modes of TWS variability through a PCA (principal component analysis;Preisendoerfer, 1988) so that we can recognize similar water mass dynamics from the several sources of information, and quantify statistically the differences between them.We interpret each of them in terms of individual hydrological structures by the distinction of ground and surface water systems and empirical climatic indexes.Analysis of GRACE signals over Australia using PCA have already been proposed in previous studies (Rieser et al., 2010;Awange et al., 2011).For further validation, we compared the regional solutions to three in situ and modelled hydrology data sets: (1) rainfall across the continent; (2) surface and subsurface storage in the Murray-Darling Basin (see Fig. 1); (3) soil moiture derived from australian hydrologic model.The conclusion of our analysis is given in Sect. 4.

Global Level 2 solutions from GRGS
We compare the solutions derived from integral energy approach to the classical global solutions computed by the GRGS in Toulouse.The Level 2 GRGS-EIGEN-GL04 10day models are derived from Level 1 GRACE measurements including KBRR and from LAGEOS-1/2 SLR data for enhancement of lower harmonic degrees (Bruinsma et al., 2010).These gravity fields are expressed in terms of normalized spherical harmonic coefficients of the geopotential (or Stokes coefficients) from degree 2 up to degree 50 using an empirical stabilization approach without any post-processing smoothing or filtering.Spherical harmonics up to degree 50 and the corresponding 10-day global equivalent water height (EWH) maps of 1 • spatial resolution which describe surface water mass change for 2002-2012, are available for the last release (RL02) at http://grgs.obs-mip.fr.
In the present study spatial means of TWS over a specific drainage basin δψ( t) for the 10-day period t, such as the Murray-Darling, is computed as the scalar product of the water mass coefficients δC nm , δS nm with A nm and B nm that are the normalized harmonic coefficients of the considered geo-graphical mask (Wahr et al., 1998;Ramillien et al., 2006b): where R is the Earth's mean radius (≈ 6371 km).TWS variations can be expressed in terms of equivalent water height changes if δψ is divided by the area of the drainage basin.
2.1.2Global Level 2 solutions from CSR, GFZ and JPL filtered using an ICA approach For further comparison, a post-processing method based on independent component analysis (ICA) was applied to the Level 2 GRACE solutions from different official providers (i.e.UTCSR, JPL and GFZ) pre-filtered with 400 km radius Gaussian filters.Before applying an ICA, the Level-2 GRACE solutions need to be somehow low-pass filtered not to have Gaussian distribution.When they are not filtered enough, they are still dominated by striping and their distribution keeps being Gaussian.If they are too low-pass filtered, they correspond the long wavelengths of the continental hydrology and they also exhibit Gaussianity.Thus, the good compromise of 350-400 km to ensure non-Gaussianity has been found by Frappart et al. (2011) (Frappart et al., 2010(Frappart et al., , 2011)).For a given month, the ICA 400 km-filtered solutions only differ from a scaling factor, so that only the GFZderived ICA 400 km-filtered solutions are presented.

Regional solutions
An alternative regional approach to the classical global one (i.e.Stokes coefficients determination) has been recently proposed to improve geographical localization of hydrological structures and to reduce leakage and aliasing.This energy integral method consists in recovering equivalent-water thickness of juxtaposed 2 • or 4 • surface tiles by inversion of differences of potential anomaly (DPA) between GRACE vehicles A and B that are related to the GRACE intersatellite KBR range (KBRR) residuals (Ramillien et al., 2011(Ramillien et al., , 2012)).In terms of computation, our regional method differs from the NASA/GFSC "mascons" (Rowlands et al., 2002(Rowlands et al., , 2005;;Lemoine et al., 2007) since, instead of having a good spatial frequency description through classical band-limited spherical harmonics, it takes the benefit of the best spatial localization of surface hydrology structures by construction (Freeden and Schreiner, 2008).This gain of resolution by using a linear Newtonian operator has been lately shown using 2 • regional solutions over South America, where the main modes of variability coincide with the geographical limits   of known hydrological units such as individual groundwater layers (Frappart et al., 2013b).These KBRR residuals were obtained by correcting the raw observations from the a priori gravitational accelerations of known large-scale mass variations (i.e.atmosphere and oceanic mass variations, polar movements, solid and oceanic tides) as well as the static gravity field of the Earth, during the iterative least-squares orbit adjustement made by the GINS (Géodésie par Intégrations Numériques Simultanées) software (Bruinsma et al., 2010).The effects of non-conservative forces measured by on-board GRACE accelerometers are also removed from the along-track observations, in order to get the contributions of unmodelled phenomena, thus mainly water mass changes in continental hydrology.In order to reduce unrealistic orbit error at fractions of the satellite revolution period, and thus to avoid instabilities in the inversion of regional solutions, the GRACE-based DPA tracks passing over Australia are linearly detrended first.This operation locally absorbs orbit error and keeps a subset of short and medium wave-  lengths of GRACE KBRR residuals; it has a unique way of dealing with short arc regionalization of a global-scale problem.Since classical gravimetric inversion does not provide a unique solution and to reduce the spurious effects of the noise in the residuals, regularization strategies have been applied to find numerically stable solutions, either based on singular value decomposition (SVD) and L-curve analysis (Ramillien et al., 2011), or by introducing averaging radius as spatial constraints (Ramillien et al., 2012).Accordingly to this last type of regularization, before least square inversion we added a spatial constraint matrix block C to the Newtonian matrix, this latter relating the equivalent-water heights (i.e.unknown parameters) to the along-track DPA (i.e.observations).The coefficients of this extra matrix block C are obtained by imposing each equivalent-water height to be a linear combination of its neighbours weighted by the inverse of their angular distances, and inside a maximum geographical radius r (e.g. in the case of spatial "averaging", the coefficients inside a given radius r should equal 1/P , as P is the number of equivalent-water heights located at a distance  lesser than r and zero elsewhere).Introduction of spatial constraints enables the ill-conditioned Newtonian matrix to be inverted, especially when the maximum radius remains large (> 800 km), but leads to smoothing and thus the risk of losing details in the regional solutions.After several tests on these filtering effects on the regional estimate (see Ramillien et al., 2012), we chose a compromise by keeping realistic hydrological details and limiting the increase of numerical noise, therefore a radius of r = 600 km is considered over continental areas, and a longer radius of r = 50 000 km over oceans in order to damp unrealistic oscillations efficiently.Following this latter regularization, time series of successive 10-day regional solutions of water mass have been produced over a region including the whole Australian continent (112 • E-156 • E; 44 • S-10 • S) for 2003-2011.The idea for cancelling the lack of long wavelengths in the final 10-day regional solutions is to complete these solutions by adding the large-scale undulations of the GRGS solutions to them (Fig. 2).For this purpose, we computed the time series of the difference between regional and GRGS solutions, both spatially averaged over Australia.This residual time series represents the largescale undulations that are unfortunately not described by the regional solutions because of the linear detrending.Of several maximum harmonic degrees n = 2, 3, 4,. . ., 12 are then used to make the synthesis of GRGS Stokes coefficients onto smooth 10-day global grids.These latter grids are then spatially averaged over the whole of Australia to produce time series for each chosen degree n.These profiles are compared to the residual time series, and the optimal maximum degree to retain is the one showing the least RMS difference with the residuals.In the particular case of our study, considering GRGS Stokes coefficients of degrees less or equal to n = 6 explains the best the residual discrepancy over Australia.This optimal harmonic degree corresponds roughly to the latitudinal dimension of the geographical region considered in the inversion (about 6700 km, here).The long wavelengths from GRGS Stokes coefficients represent about 47 and 45 % of the complete hydrological water mass signals over Australia and the Murray-Darling Basin, respectively.As presented in Fig. 3, from December 2008 to November 2009 these seasonal amplitudes of the augmented 10-day regional solutions are typically of ±250 mm of EWH in the tropical northern region of Australia, where the monsoon occurs from January to April, yielding an alternation of wet and dry seasons.As expected, these strong seasonal amplitudes are in the northern tropical band, for latitudes greater than −20 • S.Moreover, the series of the regional TWS maps reveals the constant loss of groundwater in the southeastern part of Australia, which is affected by the so-called Millennium Drought until 2010.Additionally, we show in Fig. 4

In situ and modelled hydrological data sets
Regional TWS estimates were also compared to in situ and modelled hydrological data sets.To be consistent with GRACE solutions we applied the so-called "boxcar" convolution filter of the Generic Mapping Tools (GMT) software (Wessel and Smith, 1998) which consists of spatially averaging all points inside a filter radius of 220 km in the neighbourhood of the considered grid point (or nearly 2 • , consistently with GRACE resolution).

BoM rainfall product
Rainfall over Australia, from the Bureau of Meteorology (BoM) of the Australian Government, were used for comparison against GRACE regional solutions following Rieser et al. (2010) who found a direct relationship between rainfall and GRACE surface mass changes.Rainfall observations are converted to monthly values and interpolated to geographic grids with the spatial resolution of 15 × 15 based on the Barnes successive-correction method (Jeffrey et al., 2001;Jones et al., 2009).Rainfall grids are available since 1967 using the site http://www.bom.gov.au.

Groundwater
In situ estimates of groundwater storage (GWS) in the Murray-Darling Basin were obtained from an analysis of groundwater levels observed in government monitoring bores from 2002 to 2009.Variations in groundwater storage ( S GWS ) were estimated from in situ measurements as  a Except for regional solutions: here mode 3 is used to estimate correlation values.b Except for regional solutions: here mode 2 is used to estimate correlation values.
where S y is the aquifer specific yield and H is the groundwater-level (L −1 ) observed in monitoring bores.
Groundwater level data (H ) were sourced from government departments of the states covered by the Murray Darling Basin (QLD, Natural Resources and Mines; NSW, Department of Water and Energy; VIC, Department of Sustainability and Environment; and SA, Department of Water Land and Biodiversity Conservation).Only observation bores (production bores excluded) with an average saturated zone < 30 m from the bottom of the screened interval were selected, as deeper bores can reflect processes occurring on longer timescales (Fetter, 2001).In the Murray and Darling drainage basins, changes in groundwater levels were computed at three-and six-month time steps, respectively.Only bores with at least 80 % of the time periods populated with groundwater levels were selected for the analysis (1470 bores in the Murray and 958 bores in the Darling, Fig. 1).The median of the groundwater level for each time period was first calculated at each bore and any change at a bore was computed as the difference of median groundwater level between two consecutive periods.A spatial interpolation of the groundwater level change between two consecutive periods was performed across the basin using a kriging technique (Isaaks and Srivastava, 1989).Spatial averages of groundwater level change were computed for each aquifer group.The Murray-Darling drainage basin comprises sev-eral unconfined aquifers that can be regrouped into 3 groups according to their lithology.The specific yield S y in Eq. ( 2) is estimated for each group: a range from 5 to 10 % for the clay sand unconfined aquifer group (Macumber, 1999;Cresswell et al., 2003;Hekmeijer and Dawes, 2003a;CSIRO, 2008); from 10 to 15 % for the shallow sandy clay unconfined aquifer group (Macumber, 1999;Urbano et al., 2004); and from 1 to 10 % for the fractured rock aquifer group (Cresswell et al., 2003;Hekmeijer and Dawes, 2003b;Smitt et al., 2003;Petheram et al., 2003).In situ estimates of changes in groundwater storage are calculated using the spatially averaged change in groundwater level for an aquifer group and the mean value of the specific yield for that group.The range of possible values for the specific yield was used to estimate the uncertainty.Groundwater changes in the deep, confined aquifers (mostly Great Artesian Basin and Renmark aquifers) are either due to: (1) a change in groundwater recharge at unconfined outcrop; (2) shallow pumping at unconfined outcrop; or (3) deep pumping in confined areas for farming (irrigation and cattle industry).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 (Kirby et al., 2006).To allow direct comparison between TWS and in situ groundwater estimates, pumping from the deep aquifers was added to the in situ groundwater time series, assuming the −0.42 km 3 yr −1 pumping rate remained constant during the study period.

Surface water
In the predominantly semiarid Murray-Darling Basin, most of the surface water (SW) was stored in a network of reservoirs, lakes and weirs (Kirby et al., 2006).A daily time series of the surface water storage in the network of reservoirs were obtained from the Murray-Darling Basin Commission and the state governments from 2002 to 2010.In our study we used these series as a component of total water storage estimated from the in situ measurements.

Soil moisture derived from AWRA hydrologic model
Soil moisture (SM) was derived from the Australian Water Resources Assessment (AWRA) hydrologic model (Van Dijk and Renzullo, 2011;Van Dijk et al., 2011).The system com-

PCA of the GRACE data sets
PCA was applied to the TWS time series of the 10-day GRACE-based regional, 10-day GRGS RL02 and monthly ICA 400 km-filtered solutions, as well as the BoM rainfall grids after having removed the dominant seasonal cycle using  an averaging window of 13 months.Figures 6, 8 and 10 present the spatial components of the three first modes that correspond to the most significant variability.The corresponding temporal components are presented in Figs. 7, 9 and 11.For the regional solutions, the first three modes of PCA explain variances of 60, 15 and 13 % of variability, respectively.Correlation values between the first modes of the different data sets are given in Table 1 for the spatial components and Table 2 for the temporal components, and their structures are discussed in the very next sections.

The first mode of TWS variability
The first PCA mode explains 60, 67 and 52 % of the variability of the regional solutions, GRGS and ICA 400 km-filtered TWS solutions, respectively.The first spatial component of GRGS global solutions still contains important north-south striping (Fig. 6) that degrades the amplitudes and the localization of the hydrological structures.Note that by construction the regional approach reduces this polluting noise due to orbit resonance, since its first step consists of removing the long-term satellite gravity undulations before inversion for equivalent-water heights, as explained in Sect.2.1.3.For each PCA mode, the correlations between the spatial components of GRACE-based data sets and rainfall BoM product are shown in Table 1.Cross-correlation values for temporal components are also estimated and presented in Table 2.
For the first mode, with a confidence level of 95 %, temporal pattern correlations between different GRACE solutions are very consistent (70-80 %).However, by comparing GRACE data with rainfall product only around 60 % of the signals are explained.In fact, GRACE satellites observe the mass variations due to rainfall but also other terms of the water balance equation that includes surface runoff.On the other hand, the correlation between the spatial patterns of GRACE products and rainfall drop due to the loss of mass only observed in the gravity field solutions at the Great Sandy Desert (in the northwest region of Australia) and the Murray Basin.In the Great Sandy Desert, rainfall is nearly ungauged (Van Dijk et al., 2011) and discrepancies are likely due to the missing information in the BoM products.The signal observed over Murray Basin could be associated to aquifer variations that are not well described by hydrological models such as AWRA (for more details see Van Dijk et al., 2011).The loss of mass over Murray-Darling Basin presented in the spatial component of the first mode (Fig. 6) is due to the recent Millennium Drought affecting this region (Leblanc et al., 2009(Leblanc et al., , 2011)).The corresponding negative trend on the temporal component can be related to the cumulative effect of the prolonged drought from 2003 to 2009.The positive trend of water mass changes after 2010 indicates a gain of water mass in this southeastern part of Australia.
Climate indexes are also used to achieve the validation.In fact, precedent studies show that Australian climate remains forced by these climate oscillations (Westra and Sharma, 2006).Strong correlations between oceanic indexes and GRACE TWS over Australia have been reported (e.g.García-García et al., 2011;Foorotan et al., 2012).The South Oscillation Index (SOI) and the Pacific Decadal Oscillation (PDO) interannual variations were smoothed out using the identical 13-month running average window, and then compared to the temporal components of the first modes (Fig. 7a  and b).As illustrated in Table 3, high correlations (≈ 70 %) were found between these climatic indexes and the temporal components of the first TWS modes, suggesting the important role of these climate oscillations on the Australian wa-ter mass storage variability at interannual timescales.Mantua and Hare (2002) compared the PDO index with rainfall data and suggested that warm phases of the PDO (positive values) coincide with anomalously dry periods in eastern Australia.The signature of this phenomenon can be seen in the spatial and temporal components of the first PCA mode over the Murray-Darling Basin (Figs. 6 and 7), and it is probably the reason why the large negative correlations values are obtained (see Table 3).

The second and third modes
In the case of our regional solutions, the second and third modes explain a very similar percentage of total variance (i.e. 15 and 13 %, respectively), so we decided to invert their order for consistency with the modes from other solutions.
The PCA spatial patterns of second-mode GRGS, ICA 400 km-filtered solutions, rainfall data and third-mode regional solutions are presented in Fig. 8.They all show an important and continuous deficit of mass occurring over the  Murray-Darling Basin and Lake Eyre in the centre of Australia.Looking at the temporal components in Fig. 9a, we notice that the strongest deficit happens in April 2006 and it can be interpreted as the maximum of severity of the Millennium Drought that has been observed (Leblanc et al., 2011).However, the low correlation values (spatial and temporal, see Tables 1 and 2) between GRACE products over the whole of Australia seem to indicate that the main amplitudes seen in the Murray-Darling Basin (see Fig. 8 spatial pattern and Fig. 9a temporal pattern) are probably related to the interannual groundwater level variations.Van Dijk et al. (2011) attributed the differences between GRACE data and AWRA model to deficiencies in the description of the groundwater dynamics.For our purpose, we derived interannual groundwater changes from in situ observations for comparison (see Sect. 2.2).They are obtained by computing the mean of observations over the Murray-Darling Basin, after the seasonal variations are removed using a running average window.The water mass rate over the observation period 2003-2010 are essentially presented in the first mode, as it is shown in Sect.3.1.1.Then, the linear trend is also removed from groundwater observations in order to compare the residuals to the signals of second temporal mode.In Fig. 9a and b, we notice long-term consistency between the interannual variations observed in mode 2 and the interannual groundwater level variations, especially over the Murray Basin: there is a maximum of water mass around 2006 for both, followed by a decrease until 2008 to finally reach a quite constant level.There are simple reasons that can explain the presence of the time shift of the 2006 peak: (1) the geographical area where the analysis is made is not the same, the comparison is made between groundwater interannual variations over Murray Basin and the temporal mode 2 of GRACE which includes the whole continent signal, so more signals are spatially averaged in the last case; (2) the localized distribution of bores is clearly not homogeneous since they do not cover the complete region, the shortcoming being to miss a large portion of signal energy in uncovered areas, whereas the GRACE data have lower spacial resolution but their cov- erage made by the satellites along their tracks is surely more homogeneous; (3) and the lower resolution of GWS measurements due to gaps in the starting raw observations before time averaging (6 months for the whole Murray-Darling Basin and 3 months for Murray Basin).The PCA spatial and temporal modes presented in Figs. 10 and 11 show positive variations located in the northwest part of Australia, which correspond to a gain of water mass related to an increase of rainfall (Li et al., 2009).As previously mentioned by Munier et al. (2012), an interannual cycle of TWS appears over the Canning Basin.This long-term variation in the GRACE modes, shown in Figs. 10 and 11, is consistent with the third mode of rainfall (see correlation values in Tables 1 and 2).Cross-correlation values are 64, 50 and 74 % for regional, GRGS and ICA 400 km-filtered solutions, respectively and with a time lag of 4 months with respect to rainfall data.

Detection of localized hydrological events
To demonstrate the improvement in detection brought by our regional solutions, we confronted them with known meteorological events sharply localized in both space and time.

An exceptional rainfall in 2010
Anomaly maps of TWS were generated as the difference between the monthly TWS grids corresponding to the exceptionally high rainfall in July and September 2010, and the mean grid averaged over the previous dry period (i.e.mean of the six months from June to November 2009).GRACEbased TWS and BoM rainfall anomaly grids are presented in Fig. 12.When the BoM rainfall anomaly grids are taken as reference, the TWS anomaly grids derived from regional solutions offer a better spatial localization of the September 2010 rainfall structures than the ones computed with global GRGS solutions.

Case of cyclone Charlotte
The heavy rainfall and subsequent floods that were caused by cyclone Charlotte between 9 and 12 January 2009 in the Gulf of Carpentaria basins are clearly visible in Fig. 13.Tregoning et al. (2012) already noticed the signature of this climate event in GRACE global solutions.The first flood warnings for Queensland were issued in the first week of January from moderate to major flooding in the Gulf rivers, and in Georgina and Diamantina rivers of Western Queensland (see Fig. 1).Since important rainfall continued, these warnings were maintained for nine weeks (see BoM report of 2009).Instantaneous TWS anomalies, defined as the difference of two consecutive TWS 10-day regional and GRGS RL02 solutions are presented in Fig. 13a and b, respectively.In Fig. 13c, TWS anomalies derived from ICA 400 km-filtered solutions are defined as the difference between two consecutive TWS monthly solutions.We notice clearly the strong signature of the floods caused by cyclone Charlotte, since the regional solutions reveal important positive water heights of ≈ 600 mm of EWH, higher than the levels described by the solutions of the previous months.For our regional solutions, this corresponding TWS anomaly is localized in the Carpentaria drainage basins area (see Fig. 1).This sudden meteorological event agrees quite well with the BoM rainfall database (Fig. 13d) and GRGS RL02 estimates.However, ICA 400 km-filtered solutions show a very low amplitudes of this event, as these combined solutions remain smooth.Leblanc et al. (2009) have shown the effects of drought over Murray-Darling Basin from 2002 to 2009, and they compared GRGS-RL02 global solutions to the sum of in situ data for surface and ground waters plus model output for soil moisture.Time series of GRACE TWS averaged over the whole Murray-Darling using regional, GRGS and ICA 400 km-filtered solutions are presented in Fig. 14.We notice TWS variations are consistent between all GRACE solutions (correlations values: regional/GRGS = 90 %, ICA/regional = 76 %, ICA/GRGS = 79 %).The average of rainfall in 2010, caused by a switch to a strong La Niña event, brings an end to this lasting drought and logically leads to the increase of TWS.TWS estimates from reconstructed in situ (GWS + SW) and modelled hydrologic information (SM) were interpolated onto 2 • sampling grids.According to the temporal resolution of the GWS data for the period 2003-mid-2009, the GRACE data sets were smoothed versus time: 3-month and 6-month samplings for Murray and Darling basins, respectively.The corresponding time series for the separate Murray and Darling basins are presented in Fig. 15.Correlations of 70-90% and RMS values of 15-30 mm EWH listed in Table 4 suggest a good agreement between re-constructed TWS (i.e.SM + SW + GWS) variations and the ones derived from GRACE data sets, our regional GRACE solutions giving the best agreement.As seen in Fig. 15, we used a coefficient of proportionality (i.e.Specific yield) ranging as described in Sect.2.2.2 to convert in situ measurements of water Gridded correlations and RMS differences in the Murray-Darling Basin are even better for our regional GRACE solutions (Fig. 16), especially in the Murray Basin, where the correlations are greater than 75 %, and where GWS measurements are the densest (see data coverage in Fig. 1), whereas GRGS and ICA-filtered solutions may drop down to 25 %, and even less.

Conclusions
In this paper, we validated our 10-day GRACE regional solutions of water mass variations over Australia for the period 2003-2011 by comparing them with (i) global GRGS RL02 and ICA 400 km-filtered solutions based on global spherical harmonics representation, and (ii) independent data sets, i.e. estimates of rainfall, groundwater and surface water derived from in situ observations and modelled soil moisture.
Once applied to the GRACE data sets, PCA provided the main modes of variability of these data sets before comparison.Temporal and spatial patterns are consistent for the regional and global solutions.However, the regional solutions offer better geographical localization of hydrological structures, while global GRACE solutions remain affected by aliasing errors (i.e.north-south striping).There are high correlations between SOI and PDO climate indexes during the Millennium Drought and our regional solutions, in particular for the first PCA mode at the scale of the whole of Australia.Besides, second and third modes are also related to the drought that occurred in the southeastern part of Australia up to late 2006.While all the GRACE solutions remain consistent over Australia, our regional solutions yield the best agreement with in situ water table measurements in the Murray-Darling drainage basin.This is probably due to the benefit of the gain of details dealing with regional short tracks of GRACE, instead of using band-limited harmonics functions defined on the entire terrestrial sphere.

Fig. 1 .Fig. 1 .
Fig. 1.Australian drainage regions used in our study.Geographical locations of the bores in the Murray-Darling are marked as black cercles.They represented, with at least 60 % of data period populated (because of their availability during the period 2002-2009), averages of the first 30 m of the water column height.

. 3 .Fig. 3 .
Fig. 3. Time series of regional maps of water mass from December 2008 to November 2009.We notice the annual signal with amplitude of ±250 mm of EWH representing the alternating wet and dry seasons.
Fig. 5. Time series of TWS means over the estimated desert region (Central Australia) of 18 tiles, the largest signal amplitudes are lesser than 30 mm and the RMS value is of 19 mm.

Fig. 5 .
Fig. 5. Time series of TWS means over the estimated desert region (central Australia) of 18 tiles, the largest signal amplitudes are lesser than 30 mm and the RMS value is of 19 mm.
the maps of linear trends and annual amplitudes computed over the period 2003-2009.We have chosen this period to separate El Niño and La Niña cycles (mid-2009).Even if the long-term variations over Australia are not strictly linear, the trends reveal the regions www.hydrol-earth-syst-sci.net/17/4925
Fig. 9. (a) Temporal components of PCA 3rd orthogonal mode of regional solutions and 2nd mode of GRGS, ICA 400 km-filtered solutions and Rainfall product.(b) Groundwater levels from in situ measurements over Murray-Darling Basin.

Fig. 11 .
Fig. 11.Temporal components of PCA 2nd orthogonal mode of regional solutions and 3rd mode of GRGS, ICA 400 km-filtered solutions and rainfall product.

Fig. 12 .Fig. 12 .
Fig. 12. Regional and GRGS RL02 of ∆TWS monthly maps estimated with respect to the mean map for the previous driest season (from June to November 2009), and rainfall maps.For these two periods we can notice that the regional solutions show a better localization of better the exceptional rainfall maximum ocurring in the central region of Australia.

Fig. 13 .
Fig. 13.TWS time series map derived from regional solutions (a), from GRGS RL02 (b), from ICA-filtered solutions (c) and rainfall anomalies for January 2009 (d).We notice in January solutions the strong amplitudes in the Gulf of Carpentaria due to the floods caused by cyclone Charlotte (9 January-12 January 2009).

Fig. 14 .
Fig. 14. 10-day and monthly time series of averaged total water storage over the whole Murray-Darling Basin from 2003 to 2011 computed from 10-day regional solutions (blue line), 10-day GRGS global solutions (red line) and monthly ICA 400 km-filtered solutions (green line).

Fig. 15 .
Fig. 15.(a) 3 month averaged time series of TWS over Murray Basin and (b) 6 month averaged time series of TWS over Darling Basin for 2003-2010.Series are computed from GRACE-based regional solutions (blue line), GRGS RL02 (red line), ICA 400 km-filtered solutions (green line) and TWS estimated as groundwater storage + surface water observations + modeled soil moisture (grey lines).Discountionous grey line: Maximum and minimum values of the specific yield are used (values of 0.06 and 0.14); Continiuous grey line: mean value of specific yield is used (0.01).

Table 1 .
Correlation values between PCA spatial components (mode 1, 2 and 3) of different GRACE-based data sets and rainfall BoM product.Except for regional solutions: here mode 3 is used to estimate correlation values.b Except for regional solutions: here mode 2 is used to estimate correlation values. a

Table 2 .
Cross-correlation maxima values and corresponding temporal lag in number of months (second value in italic) between PCA temporal components (mode 1, 2 and 3) of different GRACE-based data sets regional solutions, GRGS and ICA 400 km-filtered solution, as well as rainfall BoM product.Data derived from regional and GRGS are filtered and interpolated to monthly resolution.The time lag represents the time delay between water mass variations derived from GRACE and the observed rainfall.

Table 3 .
Correlation values between climate index and PCA first mode temporal patterns of regional, GRGS and ICA solutions.