On the use of GRACE intersatellite tracking data for improved estimation 1 of soil moisture and groundwater in Australia 2

8 An accurate estimation of soil moisture and groundwater is essential for monitoring the 9 availability of water supply in domestic and agricultural sectors. In order to improve the 10 water storage estimates, previous studies assimilated terrestrial water storage variation 11 ( ) derived from Gravity Recovery and Climate Experiment (GRACE) into land surface 12 models. However, the GRACE-derived was generally computed from the high level 13 products (e.g., land grid). The gridded data products are subjected to several drawbacks such 14 as signal attenuation and/or distortion caused by ad hoc posteriori filters, and a lack of error 15 covariance information. The consequence is undesired alteration of data and its 16 statistical property. To exploit the GRACE information rigorously and negate these 17 limitations, this study uses the fundamental GRACE satellite tracking Level 1B (L1B) data, 18 not the post-processed grid data. The approach combines the GRACE’s least-squares 19 normal equation (full error variance-covariance information) of L1B data with the results 20 from the Community Atmosphere Land Exchange (CABLE) model to improve soil moisture 21 and groundwater estimates. This study demonstrates, for the first time, an importance of 22 using the raw GRACE data. The GRACE-combine (GC) approach is developed for optimal 23 least-squares combination maximizing the strength of the model and observations while 24 suppressing the weaknesses. The approach is applied to estimate the soil moisture and 25 groundwater over 10 Australian river basins and the results are validated against the satellite 26 soil moisture observation and the in-situ groundwater data. We demonstrate the GC approach 27 delivers evident improvement of water storage estimates, consistently from all basins, 28 yielding better agreement at seasonal and inter-annual time scales. Significant improvement 29 is found in groundwater storage while marginal improvement is observed in surface soil 30 moisture estimates likely due to limitation of GRACE’s temporal and spatial resolution. 31


Introduction
The changes of Terrestrial Water Storage ( ) derived from the Gravity Recovery And Climate Experiment (GRACE) data products have been used in the last decade to study global water resources, including groundwater depletion in India and Middle East (Rodell et al., 2009;Voss et al., 2013), water storage accumulation in Canada (Lambert et al., 2013), flood-influenced water storage fluctuation in Cambodia (Tangdamrongsub et al., 2016).The gravity data obtained from GRACE satellites are commonly processed and released in three different product levels (L) that increase in the amount of processing, L1B -satellite tracking data (Wu et al., 2006), L2 -global gravitational Stokes coefficients (Bettadpur, 2012), and L3 -global grids (Landerer and Swenson, 2012).The original (L1B) GRACE information is Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.
inevitably altered or sheered due to data processing and successive post-processing filterings, because the error covariance information is not propagated through each post-processing step.
The GRACE-derived has been computed widely from the higher-level products (e.g., L2 and L3) on which various ad hoc post-processing filters were applied (e.g., Gaussian smoothing filter (Jekeli, 1981), destripe filter (Swenson and Wahr, 2006)).obtained from these filters lacks proper error covariance information and is attenuated and distorted.
To overcome the signal attenuation in GRACE high-level products, empirical approaches have been developed, including the application of scale factors computed from land surface models (Landerer and Swenson, 2012) to the GRACE L3 products.GRACE uncertainty in high level product is usually unknown or assumed.For example, Zaitchik et al. (2008) derived empirically a global average uncertainty that is variable depending on choices of post-processing filters (Sakumura et al., 2014).Furthermore, GRACE error and sensitivity is dependent on latitudes due to the orbit convergence toward poles (Wahr et al., 2006) and any post-processing filters will alter the GRACE data and their error information.Rigorous statistical error information is of equal importance to derivation of for data assimilation and model calibration (Tangdamrongsub et al., 2017).and its uncertainty estimates should be formulated directly from L1B data considering the complete statistical information.
The GRACE information is not fully exploited in many studies.For example, groundwater storage variation ( ) is often computed by subtracting the soil moisture variation ( ) component simulated by the land surface model from GRACE-derived data (Rodell et al., 2009, Famiglietti et al., 2011), assuming the model is error-free.This may result in the inaccurate and the associated error estimate as the uncertainties of GRACE and of the land surface model outputs are neglected in the combination of two noisy data.In data assimilation application, albeit its importance, the GRACE uncertainty is commonly derived empirically not necessarily reflecting the true GRACE error characteristics (e.g., Zaitchik et al., 2008;Tangdamrongsub et al., 2015;Tian et al., 2017).For example, Girotto et al. (2016) used L3 product and showed that it was necessary to adjust GRACE observation and its uncertainty in order to make their water storage estimates more accurate.Similarly, Tian et al. (2017) reported the need of applying a scale factor to GRACE uncertainty (from mascon product) in their GRACE assimilation process.It is apparent that the use of post-processed GRACE products often requires data tuning, leading possibly to an integration of incorrect gravity information into the data assimilation system.Some recent studies began to employ the full variance-covariance information in the data assimilation scheme (Eicker et al., 2014, Schumacher et al., 2016;Tangdamrongsub et al., 2017), however, the GRACE signal used were still affected by the post-processing filters.
This study aims to use the GRACE information of measurement directly from the raw L1B data.The approach optimally combines the GRACE's least-squares normal equations with the model simulation results from the Community Atmosphere Land Exchange (CABLE, Decker, 2015) to improve and estimates.The proposed approach presents three main advantages.Firstly, one can exploit the full GRACE signal and error information by using the normal equation data sets.Secondly, the approach is developed for optimal least-squares combination, which maximizes the model and observation strength while simultaneously supressing their weaknesses.Finally, the method bypasses empirical, multiple-step post-processing filters.The main objective of this study is to present the GRACE-combined (GC) approach to estimate improved and at regional scales.We demonstrate our approach applied to 10 Australian river basins (Fig. 1a).We validate the top layer of estimates against the satellite soil moisture observation (the Advanced Microwave Scanning Radiometer aboard EOS (AMSR-E), Njoku et al., 2003) over all 10 basins and the estimates against the in-situ groundwater data available over Queensland and Victoria (Fig. 1b, 1c).
This paper is outlined as follows: Firstly, the derivation of GC approach is presented in Sect.
2 while the description of GRACE data processing, including the use of GRACE normal equation is given in Sect.3. Secondly, the CABLE modelling is outlined in Sect. 4. This includes the derivation of model uncertainty based on the quality of precipitation data and the model parameter inputs.The processing of validation data is also described in Sect. 4.
Thirdly, Sect. 5 presents the result of and estimates and comparison to in-situ data.The long-term trends in the Australian mass variation over the last 13 years is also investigated in this section.In Sect.6, the purposed approach is discussed in terms of effectiveness, and data assimilation implementation.

A method of combining GRACE L1B data with land surface model outputs
The statistical information of computed from a model can be written as: where is the "truth" (unknown) model state vector while is the calculated state vector characterized with the model error .The model error is assumed to have zero mean and covariance .
The term is used to represent a vector including global grid, and terms with a subscript (e.g., , ) is used to represent only a regional set of (for example, in Australia).As such, the observation equation over a region can be rewritten as: As soil moisture and groundwater are the major components of in Australia (surface water storage being insignificant), the vector can be defined as: where , , represent the vectors of top (surface) soil moisture, root zone soil moisture, and groundwater storage variations, respectively.
A linearized GRACE satellite-tracking observation equation is formulated as: where is the observation vector containing the L1B inter-satellite ranging data, is the design (partial derivative) matrix relating the data and the Earth gravity field variations, contains the Stokes coefficients of time-varying geopotential fields (e.g., Wahr et al., 1998), and is the L1B data noise, which has zero mean and covariance .Eq. ( 4) can be modified explicitly in terms of soil moisture and groundwater storage variations as: where contains a factor used to convert to geopotential coefficients considering the load Love numbers (e.g., Wahr et al., 1998), converts the gridded data into the corresponding spherical harmonic coefficients, and is the operational matrix converting , , and to .This model is based on the assumption that the GRACE orbital perturbation is a result of variation on the surface, which is very particular in Australia.For convenience, the term = is used in the further derivation.If M is the number of model grid cells, N max is the maximum degree of the geopotential coefficients, and L=(N max +1) 2 -4 is the number of geopotential coefficients, the dimension of , , and are L×M, M×3M, and 3M×1, respectively.
A least-squares solution of Eq. ( 5) is given as: It can be simplified as: where = and = .(The rationales of introducing and are explained in the following section).Note that, the above derivations (Eq.( 5) -Eq.( 7)) are defined with the global grid of .For a regional application, Eq. ( 7) can be modified as: where the subscript indicates the grid only in a region of interest, and for the rest of the globe.If the number of the model grid cells associated with is J and that of the outside cells is M-J.As such, the dimensions of , , , , , are L×J, J×3J, 3J×1, L× (M-J), (M-J)×3(M-J), 3(M-J)×1, respectively.The dimension of and remain unchanged, since they are essentially from the normal equations of the original GRACE L1B data (to be discussed in the following section).
From Eq. ( 8), the normal equations associated with in the region of interest can then be written as where = and = .As seen, Eq. ( 9) is the regional representation of Eq. ( 7) where only the grid cells inside the study region are used, while the contribution from the grid cells outside the region needs to be removed or corrected.Combining the normal equation of Eq. ( 2) and Eq. ( 10), the optimal combined solution of can be resolved as follows: The computation of model covariance matrix will be discussed in Sect.4.2.The posteriori covariance of can be estimated as follows: and the uncertainty estimate of is simply calculated as: where () represents the diagonal element of the given matrix.

GRACE least-squares normal equations
In this study, the least-squares normal equations are obtained from the ITSG-Grace2016 AOD1B corrections are reduced in terms of the normal equations.These data products are usually used to compute the Earth's geopotential field to the maximum harmonic degree and order of 90, or at a spatial resolution of ~220 km.The products contain the information of the normal matrix and the vector (as shown in Eq. ( 7)) as well as the a-priori time-varying gravity field coefficients predicted with the GOCO05s solution (Mayer-Gürr et al., 2015).
Note that the solution of the ITSG-Grace2016 normal equation is the anomalous geopotential coefficient vector ( ), which is referenced to the a-priori time-varying gravity field ( ), through: where and are given.To obtain a complete gravity field variation between the study period ( term in in Eq. ( 4)), the a-priori time-varying gravity field, is firstly restored to Eq. ( 14), and the mean gravity field ( ) computed from all between January 2003 and March 2016 is then removed as follows: Therefore, in Sect. 2 (e.g., Eq. ( 7) -( 11)), the matrix remains unchanged while the vector can be simply replaced by = + ( ).

GRACE-derived products
Two monthly GRACE-derived products are also used, the CNES/GRGS Release 3 (RL3) (GRGS for short, Lemoine et al., 2015)  In this study, these two independent GRACE solutions are used for two main reasons: 1. To obtain the values outside Australia.As shown in Eq. ( 9), the vector needs to be known, which can be from the GRACE-derived solution.We use the GRGS solutions as the GRGS solution provides at a spatial resolution comparable to the normal equation data.
2. To compare with the estimates from our approaches.Both GRGS and mascon solutions are used to compare and validate our estimates.

Model setup
The extensive description of the CABLE model is given in Decker (2015) and Ukkola et al. (2016).This section describes the model setup and specific changes applied for this study.
CABLE can be used to estimate soil moisture and groundwater in terms of volumetric water content every 3 hours at a 0.5 o ×0.5 o spatial resolution.The soil moisture and groundwater storage can be simply computed by multiplying the estimates with thicknesses of various layers.For soil moisture, the thickness of 6 soil layers is 0.022, 0.058, 0.154, 0.409, 1.085, and 2.872 m, from top to bottom, respectively.The thickness of the groundwater layer is modelled to be 20 m uniformly.Recalling Eq. ( 3), is defined as the soil moisture storage variation at the top 0.022 m thick layer, while is the variation accumulated over the second to the bottom soil layers (depth between 0.022 cm and 4.6 m).
CABLE is initially forced with the data from the Global Soil Wetness Project Phase 3 (GSWP3) (Dirmeyer et al., 2011; http://hydro.iis.u-tokyo.ac.jp/GSWP3), which is currently available until December 2010.We replace GSWP3 forcing data with GLDAS data (Rodell et al., 2004) to compute the water storage changes to 2016.The forcing data used in CABLE are precipitation, air temperature, snowfall rate, wind speed, humidity, surface pressure, and short-wave and long-wave downward radiations.To investigate the impact of different forcing data, the offline sensitivity study is conducted by comparing the water storage estimates computed using:

Model uncertainty
In this study, the CABLE uncertainty is derived from 210 ensemble estimates associated with different forcing data and model parameters.The 7 different precipitation products (see Table 1) are used to run the model independently.Most products are available to present day while GSWP3, Princeton, and MERRA are only available until December 2010, December 2012, and February 2016, respectively.For each precipitation forcing, 30 ensembles are generated by perturbing the model parameters within +/-10% of the nominal values.The perturbed size of 10% is similar to Dumedah and Walker (2014).Based on the CABLE structure, the and estimates are most sensitive to the model parameters listed in Table 2.For example, the fractions of clay, sand, and silt (f clay , f sand , f silt ) are used to compute soil parameters including field capacity, hydraulic conductivity, and soil saturation which mainly affect soil moisture storage.Similarly, the drainage parameters (e.g., q sub , f p ) control the amount of subsurface runoff, which has a direct impact on root zone soil moisture and groundwater storages.
From ensemble generations, total = 210 sets of the ensemble water storage estimates ( ) are obtained: and the mean value of is computed as follows: Note that due to the absence of GSWP3, Princeton, and MERRA data, the number of ensembles reduce to = 180 after December 2010, = 150 after December 2012, and = 120 after February 2016, respectively.The mean value is removed from each ensemble member, = , and the error covariance matrix of the model is empirically computed as: The (Eq. ( 18)) and (Eq.( 19)) terms can be directly used in Eq. ( 11).
Note that the sampling error caused by finite sample size might lead to spurious correlations in the model covariance matrix (Hamill et al., 2001).The effect can be reduced by applying an exponential decay with a particular spatial correlation length to .In this study, the correlation length is determined based on the empirical covariance of model estimated .The covariance function of is firstly assumed isotropic, and it is computed empirically based on the method given in Tscherning and Rapp (1974).The distance where Since there is no satellite observed or ground measured root zone soil moisture data for meaningful comparison with our results, particularly at continental scale.Validation of at regional and continental scales is currently unachievable due to a complete lack of observations at this spatial scale.

In-situ groundwater
The in-situ groundwater level from bore measurements are obtained from 2 different ground observation networks (see Fig. 1).The data in Queensland are obtained from Department of Natural Resources and Mines (DNRM) while the data in Victoria is from Department of Environment and Primary Industries (DWPI).More than 10,000 measurements are available from each network, but the data gap and outliers are present.Therefore, the bore measurement is firstly filtered by removing the sites that present no data or data gap longer than 30 months during the study period.
To obtain the monthly mean value, the hourly or daily data are averaged in a particular month.The outliers are detected and fixed using the Hampel filter (Pearson, 2005) where the remaining data gaps are filled using the cubic spline interpolation.To obtain the groundwater level variation, the long-term mean groundwater level computed between the study period is removed from the monthly values.The groundwater level variation ( ) is then converted to using = , where is specific yield.Based on Chen et al. ( 2016), = 0.1 is used for the Victoria network.Specific yields of Queensland's network have been found ranging from 0.045 (Rassam et al., 2013) to 0.06 (Welsh 2008), and an averaged =

Model-only performance
We  1).The metric used to evaluate a goodness of fit between CABLE run and GRACE mascon estimates is the Nash-Sutcliff (NS) coefficient (see Eq. ( A1)) (Fig. 2).This section investigates the impact of the GC approach on the estimates of various water storage components.The estimate obtained from the GC approach is demonstrated in Sect.5.1, by comparing with the independent GRACE mascon solution.Figure 2 shows the GC result yields the highest NS values in all basins, outperforming all other CABLE runs.In average (AVG), the NS value increases by ~35% (0.55 to 0.74) from the MN case.The similar behaviour is also seen when compared with the GRGS GRACE solution (not shown); the average NS value increases from 0.50 to 0.74.This is not surprising as the GC approach uses the fundamental GRACE tracking data as GRACE mascon and GRGS solutions do.
Improvement of NS coefficient indicates merely the successfulness of integrating GRACE data and the model estimates.The impact of GRACE varies across the individual storage as well as across the geographical location (climate regime).In general, the major contributors to are and .
Due to a small store size (only ~2 cm thick), contributes only ~2 % to .As such, , and have greater variations, which commonly lead to greater uncertainty compared to , and therefore, the stores anticipate greater shares from the GRACE update.This behaviour is seen over all basins where the differences between CABLEsimulated and GC , and estimates are greater (compared to ).
Furthermore, the impact of GRACE on , and is different across the continent.
For example, over central and southern Australia (see e.g., LKE, MRD, NWP, SWP), the dry climate is responsible for a small amount of groundwater recharge and most of the infiltration is stored in soil compartments.In this climate condition, amplitude is significantly larger than and it plays a greater role in , and consequently, the GRACE contribution is mostly seen in component.Different behaviour is seen over the northern Australia (GOC, NEC, TIM) where amplitude are greater (~40 % of ) compared to other basins (only ~17 % of ).This is due to the sufficient amount of rainfall over the wet climate region, replenishing groundwater recharges and resulting in greater variability in .Therefore, compared to the dry climate basin, the GRACE contributes to over these basins by the larger amount.March 2016 are analysed before and after applying the GC approach (Fig. 5).For comparison, the long-term trends of derived from the mascon and GRGS solutions are also shown (Fig. 5a, 5b).From Fig. 5d, GRACE effectively changes the long-term trend estimates in most basins in a way the spatial pattern of the trend of the GC solution consistent to the mascon and GRGS solutions, while satisfying the model processes and keeping the spatial resolution.The trend of is insignificant (Fig. 5e) and the GC approach does not change (Fig. 5f).The largest adjustment is seen in and components, to be consistent with the GRACE data in most basins (Fig. 5h, 5j).

Impact on long-term trend estimates
GRACE shows significant changes in the trend estimates particularly over the northern and western parts of the continent.The model estimates around the Gulf of Carpentaria basin show a strong negative trend that is inconsistent from the GRACE data.It is found that underestimated precipitation after 2012 is likely the cause of such an incompatible negative trend (see Fig. 3d).Applying the GC approach clearly improves the trend (Fig. 5c vs. 5d).The other example is seen over the western part of the continent (see rectangular area in Fig. 5c, 5d) where the averaged long-term trend of was predicted to be -0.4 cm/year but changed to be -1.2 cm/year (see also Sect.5.4) by the GC approach.
The precipitation over the western Australia is understood to be overestimated after 2012, evidently seen by that the model is always greater than the GC solution (see e.g., Fig. 3h, 4d, 4p).The GC approach reveals that the water loss over the western Australia is at least twice greater than what has predicted by the CABLE model.
In addition, the shortage of water storage in the south-eastern part of the continent from the millennium drought (McGrath et al., 2012) has been recovered (seen as a positive water storage trend in Fig. 5) after the rainfall between 2009 and 2012, while the western part is still drying out (seen as negative trends).The trend estimates in terms of mass change is discussed in more detail in Sect.5.4.

Reduction of uncertainty
Influenced by climate pattern, the uncertainty of water storage estimates significantly varies across Australia.The uncertainty of the model estimate is computed from the variability induced by different precipitation and model parameters while the uncertainty of GC solution is computed using Eq. ( 13).As expected, larger uncertainties are observed in and than in (an order of magnitude smaller) since is smaller than others (Fig. 6).Over the wet basins, larger amplitude of the water storage leads to larger uncertainty, seen over Gulf of Carpentaria, North East Coast, South East Coast, and Timor Sea where the CABLE-simulated uncertainty is approximately 28 % larger than other basins.The smaller uncertainty is found over the dry regions (e.g., LKE, SWP).In most basins, the uncertainty of is larger than the , except the wet basins (e.g., GOC, NEC, TIM) where the greater groundwater recharge leads to a larger uncertainty of .
Figure 6 demonstrates how much the formal error of each of storage components is reduced by the GC approach.Overall, the estimated CABLE uncertainties averaged over all basins (AVG) are 0.2, 4.0, 4.0, and 5.7 cm for , , , and , respectively.With the GC approach, the uncertainties of , , , and decrease by approximately 26%, 35%, 39%, and 37%, respectively.
It is worth mentioning that the model uncertainty is mainly influenced by the meteorological forcing data.The uncertainty of precipitation derived from 7 different precipitation products is shown in Fig. 6e.The spatial pattern of the precipitation uncertainty is correlated with the uncertainty of water storage estimates.The larger water storage uncertainty is deduced from the larger precipitation uncertainty.The quality of precipitation forcing data is found to be an important factor to determine the accuracy of water storage computation.

Soil moisture
The estimates are compared with the AMSR-E derived soil moisture.The processing of AMSR-E data is described in Sect 4.3.1.The performance is assessed using Nash-Sutcliff coefficients, given in Table 3.In general, CABLE (MN) shows a good performance in the top soil moisture simulation showing NS value of >0.4 for most of the basins.The top soil moisture estimate shows slightly better agreement with the C-band measurement of the AMSR-E product.This is likely caused by the greater emitting depth of the C-band measurement (~1 cm), which is closer to the depth of the top soil layer (~2 cm) used in this study (Njoku et al., 2003).
The GC approach leads to a small bit of improvement of the top soil estimate consistently from C-and X-band measurements and from all basins.No degradation of the NS value is observed in the GC solutions.The largest improvement is seen over LKE and NEC, where NS increases by 10 -15%.For other regions, the change in the NS coefficient may be incremental.

Groundwater
The estimates from the model and the GC method are compared with the in situ data obtained from 2 different ground networks in Queensland and Vitoria.For each network, all data inside the groundwater network boundary (see polygons in Fig. 1) are used to compute the average time series.From the comparison given in Fig. 7, it is found that the GC solutions of follows the overall inter-annual pattern of CABLE but with a considerably larger amplitude.This results in a better agreement with the in situ data seen from both networks.The NS coefficient of between the estimates and the in situ data are given in Table 4.The CABLE performs significantly better in Queensland (NS = ~0.5)than Victoria (NS = ~0.3).Significant improvement is found from the GC solutions in both networks, where the NS value increases from 0.5 to 0.6 (~ 22 %) in Queensland and from 0.3 to 0.6 (~85 %) in Victoria.Even greater improvement is seen when the inter-annual patterns are compared.The NS value increase from 0.5 to 0.7 (~ 32 %), and 0.4 to 0.8 (~93 %) in Queensland and Victoria, respectively.
The comparison of the long-term trend of is also evaluated.The estimated trends in Queensland and Victoria are given in Table 4. Beneficially from the GC approach, the trend is improved by approximately 20 % (from 0.4 to 0.6, compared to 1.6 cm/year) in Victoria is mainly due to the highly-demanded groundwater consumption by agriculture and domestic activities (van Dijk et al., 2007;Chen et al., 2016).As the groundwater consumption is not parameterized in CABLE, the decreasing of estimate cannot properly captured in the model simulation.Applying GC approach effectively reduces the model deficiency and improves the quality of the groundwater estimations.

Assessment of mass variation in the past 13 years
Australia experiences significant climate variability; for example the millennium drought starting from late '90 (Van Dijk et al., 2013) and extremely wet condition during several La Niña episodes (Trenberth 2012;Han 2017).These periods are referred as "Big Dry" and "Big Wet" (Ummenhofer et al., 2009;Xie et al., 2016).To understand the total water storage During Big Wet (2010 -2012), the basins like LKE, MRD and TIM exhibit the significant total storage gain of >100 Gton.The gain is particularly larger in over the basins that experienced the significant loss during Big Dry.For example, over LKE and MRD, the gain of is approximately 2 -3 times greater than .It implies that most of the infiltration (from the 2009 -2012 La Niña rainfall) is stored as soil moisture through the long drought period, and that the groundwater recharge is secondary to the increase.
The opposite behaviour is observed over the basins (such as NEC and GOC) that experienced mass gain during Big Dry.The water storage gain is greater in compared to .In During the post-Big Wet period (2012 and afterwards), the decreasing trend of water storage is observed from all basins (see Fig. 3, 4).This is mainly caused by the decrease in precipitation after 2012 and by gradual water loss through evapotranspiration (Fasullo et al., 2013).
The overall water storage change in the last 13 years demonstrates that the severe water loss from most basins during Big Dry (the millennium drought) is balanced with the gain during Big Wet (the La Niña).The negative estimated during Big Dry becomes positive in LKE, MRD, and SEC and less negative in TIM, and the greatest gain is observed from NEC by ~50 Gton during 13 year-period (see Fig. 8c).However, the water mass loss is still detected over the western basins (e.g., IND, NWP, SWP, SWC), and their magnitudes are even larger than the mass loss during Big Dry.For example, the greatest loss of ~79 Gton is observed over NWP, which is ~25 Gton greater than the loss during Big Dry (see Fig.

Comparison of GC approach with alternatives
The simplest approach to estimate is to subtract the model soil moisture component from GRACE data, without considering uncertainty in the model output, as used in Rodell et al. (2009) and Famiglietti et al. (2011).This method is called Approach 1 (App 1).
In Approach 2 (App 2) as in Tangdamrongsub et al. (2017), by accounting for the uncertainty of model outputs and GRACE data, the water storage states are updated through a Kalman filter: where , , are described in Sect.2, is an observation vector containing GRACEderived , and is an error variance-covariance matrix of the observation.The GRACE-derived and its error information is obtained from the mascon solution.The matrix is a (diagonal) error variance matrix since no covariance information is given in the mascon product.

GC data assimilation approach
We so far discussed the GC approach to update the water storage estimates independently every month.The approach can be easily expanded to sequentially update the model initial states whenever the GRACE observation is available (for example, every day) as in data assimilation (DA) like ensemble Kalman filter (Evensen, 2003) and particle filter (Weerts and El Serafy, 2006).We briefly describe a way of modifying the GC approach suitable for DA.The ensemble of simulated monthly water storage estimates is predicted based on the set of ensemble forcing data and model parameters.This is simply running CABLE for K (number of ensemble) times.When GRACE observation is available, the updated state is computed: where the subscript represents the ensemble or perturbed version of the original vector or matrix (see e.g., Eq. ( 11)).The dimension of , , is 3J×K.The estimated can be directly used as in the initial state for the next time step for CABLE run (Eicker et al., 2014;Tangdamrongsub et al., 2015;Tian et al., 2017), or used in the repeated run to avoid any spurious jump of the water storage estimates between the each step (Forman et al., 2012;Tangdamrongsub et al., 2017).This sequential update process can be carried out as long as desired.The feasibility of GRACE DA has been demonstrated with "devised" uncertainty (covariance) information.As a future work, we will develop new DA approach on the basis of full error information of GRACE data by using the least-squares normal equation and thus carrying the error information from the fundamental (satellite tracking) data level.

Conclusion
This study presents an approach of combining the raw GRACE observation with model simulation to improve water storage estimates over Australia.Distinct from other methods, we exploit the fundamental GRACE satellite tracking data and the full data error variancecovariance information to avoid alteration of signal and measurement error information present in higher level data products.
We compare groundwater storage estimates from GC approach and two other approaches, subject to inclusion of GRACE uncertainty in calculation.Validating three results of against the in situ groundwater data, we find that the GC approach delivers the most The land surface model we used is deficient in anthropogenic groundwater consumption.The model calibration will never help and the groundwater consumption must be brought in by external sources.On the contrary, the statistical approach like our GC approach may be useful to fill in the missing component and lead to a more comprehensive water storage inventory.
However, it is difficult to constrain different water storage components by only using total storage observation like GRACE.In addition, it is challenging to improve surface soil moisture varying rapidly in time, using a monthly mean GRACE observation.Tian et al.
(2017) utilized the satellite soil moisture observation from the Soil Moisture and Ocean Salinity (SMOS, Kerr et al., 2001) in addition to GRACE data for their data assimilation and showed a clear improvement in the top soil moisture estimate.The GC approach with complementary observations at higher temporal resolution should be considered particularly to enhance the surface soil moisture computation.
Finally, the GC approach can be simply extended for GRACE data assimilation.Assimilating the raw GRACE data into land surface models like CABLE enables the model state and parameter to be adjusted with the realistic error information, allowing reliable storage computation.The GC data assimilation will be developed in our future study.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.= + ; ~ ( , ), Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.1.All 8 forcing data components of GSWP3, 2. GSWP3 data with replacing one component obtained from GLDAS forcing data.It is found that the water storage estimate is most sensitive to the replacement of precipitation data, as expected, and relatively less sensitive to the change of other forcing components.We use the GLDAS forcing data in this study and also further test 7 different precipitation data products (see more details in Sect.4.2).The forcing data are up/down sampled to a 0.5 o ×0.5 o spatial grid to reconcile with the CABLE spatial resolution.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License. the maximum value of the variance decreases to half is defined as the correlation length.The obtained values vary month-to-month, and the mean value of 250 km is used in this study.It is emphasized that the model omission error caused by imperfect modelling of hydrological process within the LSM is not taken into account in the above description.We assume for such model error by increasing 20% of the model covariance.(i.e., multiplying by 1.2).4.3 Validation data4.3.1 Satellite soil moisture observationThe satellite observed surface soil moisture data is obtained from the Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) cooperating the Land ParameterRetrieval Model(Njoku et al., 2003).The observation is used to validate our estimates of top soil moisture changes ( ).The AMSR-E product provides volumetric water content in the top layer derived from a passive microwave data (from NASA EOS Aqua satellite) and forward radiative transfer model.In this study, the level 3 product, available daily between June 2002 and June 2011 at 0.25 o ×0.25 o spatial resolution is used (Owe et al., 2008).The measurements from ascending and descending overpasses are averaged for each frequency band (C and X).Then, the monthly mean value is computed by averaging the daily data within a month.To obtain the variation of the surface soil moisture, the long-term mean between June 2002 and June 2011 is removed from the monthly data.Regarding the different depth measured in CABLE and AMSR-E, the CDF-matching technique (Reichle and Koster, 2004) is used to reduce the bias between the top soil moisture model and the observation.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.0.05 is used in this study.Finally, the mean value computed from all data (in each network) is used to represent the in-situ data of the network.

Figure 2
Figure 2 demonstrates CABLE varies noticeably by precipitation as well as locations.The area-weighted average values (see Eq. (A2)) computed from Princeton, GSWP3, andTRMM yields the model reasonably agreeing with GRACE by giving the NS coefficient greater than 0.45, while MERRA, PERSIANN, and GLDAS show NS = ~0.3.The less agreement is mainly due to the quality of rainfall estimates over Australia.The NS of ECMWF is around 0.4.All model ensembles are consistent with the GRACE data over Timor Sea and inner parts of Australia (e.g., LKE, MRD, NWP) where the NS value can reach as high as 0.9 (see e.g., TRMM over TIM).On the contrary, the less agreement is found mostly over the coastal basins.Very small or even negative NS values indicate the misfit between CABLE and GRACE mascon solutions, and they are observed over Indian Ocean (see GLDAS), North East Coast (see GSWP3, PERSIANN, TRMM), South East Coast (see MERRA, TRMM), South West Coast (see GSWP3, GLDAS, MERRA), and South West Plateau (see MERRA).By averaging all estimates from 7 different precipitation, the mean-ensemble estimate (MN) delivers the best agreement with GRACE as seen by the highest average NS value (MN of AVG = 0.55) among all ensembles.Particularly, NS values are greater than 0.4 in all basins and no negative NS values are presented in MN.In average, it can be clearly seen that using the mean value (MN) is a viable option to increase the overall performance of the estimates.Therefore, only CABLE MN result will be used in further analyses.The comparison with the GRGS GRACE solution was also evaluated (not shown here) and the overall results are similar to Fig.2.5.2 Impact of GRACE on storage estimates5.2.1 Contribution of GRACE

Figure
Figure 3 shows the GC results of as well as , , and in different basins.The monthly time-series and the de-seasonalized time-series are shown.In general, GRACE tends to increase when the model (MN) is predicted to be underestimated (see e.g., LKE, MRD, NWP, SWP, TIM between 2011 and 2012) and by decrease when determined to be overestimated (see all basins between 2008 and 2010).A clear example is seen over Gulf of Carpentaria (Fig. 3d), where CABLE overestimates and produces phase delay between 2008 and 2010.The over estimated amplitude and phase delay seen in CABLE during this above period (Fig. 3c) is caused by an overestimation of soil and groundwater storage.The positively biased soil and groundwater storage causes a phase delay by increasing the amount of time required for the subsurface drainage (baseflow) to reduce to soil and groundwater stores.The overestimation of water storage is the result of overestimated precipitation or underestimated evapotranspiration.The amplitude and phase of the water storage estimate are adjusted toward GRACE observation in the GC approach.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.Queensland.Increasing of is mainly influenced by the large amount of rainfall during the 2009 -2012 La Niña episodes (see Fig. 7a).In Victoria, significant improvement of trend by about 76 % (from 0.1 to -0.2, compared to -0.3 cm/year) is observed.Similar improvement of long-term trend estimates is seen in de-seasonalized time series (improves by ~15 % in Queensland and by ~74 % in Victoria).Decreasing of in

(
mass) variation influenced by these two distinct climate variabilities, the water storage change obtained from the GC approach during Big Dry and Big Wet is separately investigated over 10 basins.The time window between January 2003 and December 2009 is defined as the Big Dry period while between January 2010 and December 2012 is defined as the Big Wet period following Xie et al. (2016).In each period, the long-term trends of GC estimates of , , , and are firstly calculated.Then, the total water storage variation (in meter) is simply obtained by multiplying the long-term trend (in m/year) with the number of years in the specific period, 7 years for Big Dry and 3 years for Big Wet.To obtain the mass variation, the water storage variation is multiplied by the area of the basin and the density of water (1000 kg/m 3 ).The estimated mass variations during Big Dry and Big Wet are displayed in Fig. 8.The long-term mass variation of the entire period (January 2003 -March 2016) is also shown.During Big Dry (2003 -2009), a significant loss of total storage (40 -60 Gton over 7 years) is observed over LKE, MRD, NWP, and SWP basins.The largest groundwater loss of >20 Gton is found from LKE and MRD.No significant change is observed over the tropical climate regions (e.g., GOC, NEC).The mass loss mostly occurs in the root zone and groundwater compartments where the sum of and explains more than 90% of the value.The mass loss is also observed in but >10 times smaller than and .
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.NEC, gain is ~8 times larger than during Big Wet.The soil compartment may be saturated during Big Dry and additional infiltration from the Big Wet precipitation leads to an increased groundwater recharge.The loss observed over GOC is simply caused by the timing selection of Big Wet period, which ends earlier (~2011) in GOC than in other basins.The gain becomes ~26 Gton if the Big Wet period is defined as 2008 -2011.
8a and 8c).The basin is less affected by the La Niña, and the rainfall during Big Wet is clearly inadequate to support the water storage recovery in the basin.Rainfall deficiency also reduces the groundwater recharge, resulting in even more decreasing of , compared to the Millennium Drought period (see Fig. 8j and 8l).The continual decrease in water storage over western basins is likely caused by the interaction of complex climate patterns like El Niño Southern Oscillation, Indian Ocean Dipole, and Southern Annular Mode cycles (Australian Bureau of Meteorology, 2012; Xie et al., 2016).
Fig. 9.It is clearly seen that from App1 are overestimated while the one from App2fits the ground data significantly better.This behaviour was also seen in Tangdamrongsub et Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.accurate groundwater estimate, followed by the approach based on incomplete information of GRACE's data error.The poorest estimate of groundwater storage is seen when the GRACE uncertainty is completely ignored.This confirms the critical value of using the complete GRACE signal and error information at the raw data level.The analysis of water storage change between 2003 and 2016 reveals that half of the continent (5 out of 10 basins) is still not fully recovered from the Millennium Drought.The TWS decrease in Western Australia has been most characteristic and the GC approach finds that the water loss mainly occur in groundwater layer.Rainfall inadequacy is attributed to the continual dry condition, leading to a greater decreasing of groundwater recharge and storage over Western Australia.

Figure 1 .
Figure 1.(a) Geographical location of 10 Australian river basins.Red and blue polygons indicate the boundaries of groundwater networks in Queensland (b) and Victoria (c), respectively.Triangles (in b and c) represent the selected bore locations used in this study.

Figure 2 .
Figure 2. NS coefficients between the model and GRACE-mascon over 10 Australian basins (in ordinate).The NS values were computed based on CABLE computed with 7 different precipitation data (in abscissa), GSWP3 (GS), GLDAS (GL), ECMWF (EC), MERRA (MR), PERSIANN (PR), TRMM (TR).The NS value of the mean estimates (the average of 7 variants) is also shown (MN).The area-weighted average NS value over all basins is also shown (AVG).The NS value of from the GRACE-combined (GC) approach is shown in the last column.The full name of the basins can be found in Fig. 1.

Figure 3 .
Figure 3.The monthly time series of , , , and estimated from model (blue) and GC (red) solutions over Gulf of Carpentaria (GOC), Indian Ocean (IND), Lake Eyre (LKE), Murray-Darling (MRD), and North East Coast (NEC).The deseasonalized time series is also shown.

Figure
Figure 6.Uncertainties of , , , and estimated from the model (blue) and the GC solutions (red) in 10 different Australian basins.The uncertainty of the precipitation is shown in (e).The area-weighted average value (AVG) is also shown.

Figure 7 .
Figure 7.The monthly time series of estimated from the model, GC solutions, and measured from the in situ groundwater network in Queensland (a) and Victoria (b).Deseasonalized time series are shown in thick lines.

Figure 9 .
Figure 9. estimated from Approach 1 (App1) and Approach 2 (App2) in Queensland (a) and Victoria (b).The in-situ groundwater network data and the GC solutions are also shown.De-seasonalized time series are shown in thick lines.
o using the nearest grid values.

Table 1 .
This work is funded by University of Newcastle to support NASA's GRACE and GRACE Follow-On missions.Mark Decker was supported by ARC Centre of Excellence for Climate Systems Science.We thank Torsten Mayer-Gürr for GRACE data products in the form of the least-squares normal equations.Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.Precipitation data from 7 different products used in this study, the Global Soil Wetness Project Phase 3 (GSWP3), the Global Land Data Assimilation System (GLDAS), the Tropical Rainfall Measuring Mission (TRMM), the Modern-Era Retrospective Analysis for Research and Applications (MERRA), the European Centre for Medium-Range Weather Forecasts (ECMWF), the Princeton's Global Meteorological Forcing Dataset (Princeton), and the Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks (PERSIANN).The temporal resolution of all products are 3 hours.Most products are available to present while GSWP3, MERRA, and Princeton terminate earlier.

Table 2 .
Ukkola et al. (2016)t are sensitive to SM and GWS estimates.The following parameters were perturbed using the additive noise with the boundary conditions given in the last column.The further parameter description can be found inDecker (2015)andUkkola et al. (2016).Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-318Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.