3.1 GRACE least-squares normal equations

In this study, the least-squares normal equations are obtained from the ITSG-Grace2016 products between January 2003 and March 2016. All L1B data including K-band range (KBR) inter-satellite tracking data, attitude, accelerometer, GPS-based kinematic orbit data, and 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 N and the vector c (as shown in Eq. 4) 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 (Δx), which is referenced to the a priori time-varying gravity field (x₀), through:

where d and x₀ are given. To obtain a complete gravity field variation between the study period (x term in Eq. 4), the a priori time-varying gravity field, x₀, is firstly restored to Eq. (12), and the mean gravity field (${\stackrel{\mathrm{‾}}{\mathit{x}}}_{\mathrm{0}}$) computed from all x₀ between January 2003 and March 2016 is then removed as follows:

Therefore, in Sect. 2 (e.g. Eq. 5), the matrix N remains unchanged while the vector c can be simply replaced by c = d + N(x₀ − ${\stackrel{\mathrm{‾}}{\mathit{x}}}_{\mathrm{0}})$.

3.2 GRACE-derived ΔTWS products

Three monthly GRACE-derived ΔTWS products are also used: the ITSG-Grace2016 DDK5 solution (ITSG-DDK5 for short, http://icgem.gfz-potsdam.de/series/99_non-iso/ITSG-Grace2016), the CNES/GRGS Release 3 (RL3) (GRGS for short, Lemoine et al., 2015; http://grgs.obs-mip.fr/grace/variable-models-grace-lageos/grace-solutions-release-03), and the JPL RL05M mascon-CRI version 2 product (mascon for short, Watkins et al., 2015; Wiese et al., 2016; http://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons). The ITSG-DDK5 product is the post-processed version of the ITSG L2 solution where the non-isotropic filter DDK5 (Kusche et al., 2009) is applied. The DDK5 solution is empirically selected here to be a good balance between the over-smoothed (e.g. DDK1) and noisy (e.g. DDK8) solutions. The GRGS solution provides ΔTWS at 1^∘ × 1^∘ globally, derived from the Earth's geopotential coefficients up to the maximum degree and order 80, and no filter nor scale factor is applied (L2 data product). Mascon provides ΔTWS at equal-area 3^∘ spherical cap grid globally. In contrast to the ITSG-DDK5 and GRGS solutions, the mascon uses a gain factor derived from the land surface model to restore mitigated signals and reduce leakage errors (L3 data products) (Watkins et al., 2015; Wiese et al., 2016). Additionally, mascon provides the ΔTWS uncertainty together with the solution. The uncertainty is computed based on several geophysical models (see Watkins et al., 2015, and Wiese et al., 2016, for more details). The uncertainty information is not available in the ITSG-DDK5 or GRGS product.

The GRACE data are obtained between January 2003 and March 2016. After retrieval, the long-term mean value between January 2003 and March 2016 is computed and subtracted from the monthly products. To be consistent with CABLE grid spacing (see Sect. 4), the ΔTWS is computed using 0.5^∘ spatial resolution. The coarse-scale datasets (e.g. mascon, GRGS) are resampled to 0.5^∘ × 0.5^∘ using the nearest grid values.

In this study, the independent GRACE solutions are used for two main reasons:

1. To obtain the ΔTWS values outside Australia. As shown in Eq. (7), the ${\widehat{\mathit{h}}}_{\mathrm{o}}$ vector needs to be known, which can be from the GRACE-derived ΔTWS solution. We use the GRGS solutions as the GRGS solution is not subject to the filter choice and it provides ΔTWS at a spatial resolution comparable to the normal equation data.

2. To compare with the ΔTWS estimates from our approaches. All solutions are used to compare and validate our ΔTWS estimates.

4.1 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 to this study. CABLE is a public available land surface model and can be used to estimate soil moisture and groundwater in terms of volumetric water content every 3 h at a 0.5^∘ × 0.5^∘ 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 six 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), ΔSM_top is defined as the soil moisture storage variation at the top 0.022 m thick layer, while ΔSM_rz is the variation accumulated over the second to the bottom soil layers (depth between 0.022 and 4.6 m).

Table 1Precipitation data from seven 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), 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 is 3 h. Most products are available to present while GSWP3, MERRA, and Princeton terminate earlier.

Download Print Version | Download XLSX

Table 2Model parameters that 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 in Decker (2015) and Ukkola et al. (2016).

Download Print Version | Download XLSX

CABLE is initially forced with the data from the Global Soil Wetness Project Phase 3 (GSWP3), which is currently available until December 2010 (http://hydro.iis.u-tokyo.ac.jp/GSWP3, https://doi.org/10.20783/dias.501). 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

1. all eight forcing data components of GSWP3 and

2. GSWP3 data with the replacement of 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 in other forcing components. We use the GLDAS forcing data in this study and also further test seven different precipitation data products (see more details in Sect. 4.2). The forcing data are up- or downsampled to a 0.5^∘ × 0.5^∘ spatial grid to reconcile with the CABLE spatial resolution.

Figure 2Histograms of the model errors computed from 210 ensemble members (${\mathcal{H}}_{\mathrm{R}}^{\prime }$) without the mean. The basin averaged values (from all 10 Australian basins) of January 2003, for example, are shown.

Download

4.2 Model uncertainty

In this study, the CABLE uncertainty is derived from 210 ensemble estimates associated with different forcing data and model parameters. The seven 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 ΔSM and ΔGWS 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 K = 210 sets of the ensemble water storage estimates (h_e) are obtained:

and the mean value of 𝓗_R is computed as follows:

Note that due to the absence of GSWP3, Princeton, and MERRA data, the number of ensembles reduces to K = 180 after December 2010, K = 150 after December 2012, and K = 120 after February 2016, respectively. The GC approach assumes that model errors are normally distributed with zero mean. Any violation of this assumption will yield a bias in the combined solutions. Therefore, the mean value is removed from each ensemble member, ${\mathcal{H}}_{\mathrm{R}}^{\prime }$ = 𝓗_R − ${\stackrel{\mathrm{̃}}{\mathit{h}}}_{\mathrm{R}}$, and the error covariance matrix of the model is empirically computed as

The ${\stackrel{\mathrm{̃}}{\mathit{h}}}_{\mathrm{R}}$ (Eq. 16) and C_R (Eq. 17) terms can be directly used in Eq. (9). The distribution of model errors is demonstrated in Fig. 2. The figure illustrates the histogram of model errors (${\mathcal{H}}_{\mathrm{R}}^{\prime }$) computed using 210 ensemble members of the model-estimated ΔSM and ΔGWS in January 2003. The histogram indicates that the model error may be approximately described by a normal distribution as introduced in Eq. (1).

Furthermore, in practice, 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 C_R. In this study, the correlation length is determined based on the empirical covariance of model-estimated ΔTWS. The covariance function of ΔTWS is firstly assumed isotropic, and it is computed empirically based on the method given in Tscherning and Rapp (1974). The distance where 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 processes within the LSM is not taken into account in the above description. The omission error may increase the model covariance and introduce a bias as well. We account for the omission error by increasing 20 % of the model covariance. (i.e. multiplying C_R by 1.2). We determine such omission error based on trial and error such that it increases the model error (due to the omission error) but does not exceed the model error value reported by Dumedah and Walker (2014). We acknowledge that this is only a simple practical way of accounting for the omission error into the total model error.

4.3 Validation data

5.1 Model-only performance

We study the model ΔTWS changes under different meteorological forcing and land parameterization. A total of 210 estimates of monthly TWS (sum of SM_top, SM_rz, and GWS) are obtained between January 2003 and March 2016 from the ensemble run based on seven different precipitation inputs. Then, the averaged values of the TWS estimates are computed from the 30 precipitation-associated ensemble members. This results in seven sets of monthly mean TWS estimates from seven different precipitation data. For each set, the monthly ΔTWS is computed by removing the long-term mean computed between January 2003 and March 2016.

The precipitation-based ΔTWS values are then compared with the GRACE mascon solution (see Sect. 3.2) over 10 different Australian basins. The comparison is carried out between January  2003 and March 2016. Due to the availability of the data, the periods used are shorter in cases of GSWP3, Princeton, and MERRA precipitation (see Table 1). The metric used to evaluate a goodness of fit between the CABLE run and GRACE mascon estimates is the Nash–Sutcliffe (NS) coefficient (see Eq. B1) (Fig. 3).

Figure 3 demonstrates CABLE ΔTWS varies noticeably by precipitation as well as locations. The area-weighted average values (see Eq. B2) computed from Princeton, GSWP3, and TRMM yield the model ΔTWS, reasonably agreeing with GRACE by giving an NS coefficient greater than 0.45, while MERRA, PERSIANN, and GLDAS show NS = ∼ 0.3. The lower agreement is mainly due to the quality of rainfall estimates over Australia. The NS of ECMWF is around 0.4.

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

Download

All model ensembles are consistent with the GRACE data over the Timor Sea and inner parts of Australia (e.g. Lake Eyre, LKE; Murray–Darling, MRD; North West Plateau, NWP) where the NS value can reach as high as 0.9 (see e.g. TRMM over TIM, Timor Sea). On the contrary, the lower agreement is found mostly over the coastal basins. Very small or even negative NS values indicate the misfit between the CABLE and GRACE mascon solutions, and they are observed over the Indian Ocean (see GLDAS), North East Coast (NEC; see GSWP3, PERSIANN, TRMM), South East Coast (see MERRA, TRMM), South West Coast (SWC; see GSWP3, GLDAS, MERRA), and South West Plateau (SWP; see MERRA).

By averaging all ΔTWS estimates from seven different precipitation datasets, 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. On average, it can be clearly seen that using the mean value (MN) is a viable option to increase the overall performance of the ΔTWS estimates. Therefore, only the 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. 3.

Figure 4The monthly time series of ΔSM_top, ΔSM_rz, ΔGWS, and ΔTWS estimated from model (blue) and GC (red) solutions over the 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.

Download

Figure 5Similar to Fig. 3, but estimated over the North West Plateau (NWP), South East Coast (SEC), South West Coast (SWC), South West Plateau (SWP), and Timor Sea (TIM).

Download

5.2 Impact of GRACE on storage estimates

5.2.2 Impact on long-term trend estimates

The spatial patterns of the long-term trends of water storage changes over January 2003 and March 2016 are analysed before and after applying the GC approach (Fig. 6). For comparison, the long-term trends of ΔTWS derived from the ITSG-DDK5, mascon, and GRGS solutions are shown in Fig. 7. From Fig. 6b, GRACE effectively changes the long-term trend estimates in most basins in a way the spatial pattern of the ΔTWS trend of the GC solution consistent to the GRACE solutions while satisfying the model processes and keeping the spatial resolution. The trend of ΔSM_top is insignificant (Fig. 6c) and the GC approach does not change (Fig. 6d). The largest adjustment is seen in ΔSM_rz and ΔGWS components, to be consistent with the GRACE data in most basins (Fig. 6 and h).

GRACE shows significant changes in the ΔTWS trend estimates particularly over the northern and western parts of the continent (Fig. 7). 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. 4d). Applying the GC approach clearly improves the trend (Fig. 6a vs. Fig. 6b). The other example is seen over the western part of the continent (see rectangular area in Fig. 6a and b) where the averaged long-term trend of ΔTWS was predicted to be −0.4 cm yr⁻¹ but changed to be −1.2 cm yr⁻¹ (see also Sect. 5.4) by the GC approach. The precipitation over Western Australia is understood to be overestimated after 2012, evidently seen by the fact that the model ΔTWS is always greater than the GC solution (see e.g. Figs. 4h, 5d and p). The GC approach reveals that the water loss over Western Australia is at least 2 times greater than what has been 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. 6) 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 are discussed in more detail in Sect. 5.4.

Figure 6Long-term trends of ΔTWS (a, b), ΔSM_top (c, d), ΔSM_rz (e, f), and ΔGWS (g, h) estimated from the model-only (left) and the GC solutions (right). The eastern part of North West Plateau basin is shown as a rectangle polygon in (a) and (b).

Download

Figure 7Long-term trends of GRACE-derived ΔTWS from ITSG-DDK5 (a), mascon (b), and GRGS solution (c).

Download

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

Download

5.3 Comparison with independent data

5.3.2 Groundwater

The ΔGWS estimates from the model and the GC method are compared with the in situ data obtained from two different ground networks in Queensland and Victoria. For each network, all ΔGWS data inside the groundwater network boundary (see polygons in Fig. 1) are used to compute the average ΔGWS time series. From the comparison given in Fig. 9, it is found that the GC solutions of ΔGWS follow the overall inter-annual pattern of CABLE but with a considerably larger amplitude. This results in a better agreement with the in situ ΔGWS data seen from both networks. The NS coefficient of ΔGWS between the estimates and the in situ data are given in Table 4. The CABLE ΔGWS 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 increases from 0.5 to 0.7 (∼ 32 %) and 0.4 to 0.8 (∼ 93 %) in Queensland and Victoria, respectively.

Table 3NS coefficients between top soil moisture estimates and the satellite soil moisture observations from AMSR-E products over 10 different Australian basins. The area-weighted average value (AVG) is also shown.

Download Print Version | Download XLSX

Table 4NS coefficient and long-term trend of ΔGWS estimated from the model-only and GC solutions in Queensland and Victoria groundwater network. The long-term trend of the in situ data is also shown.

Download Print Version | Download XLSX

The comparison of the long-term trend of ΔGWS is also evaluated. The estimated trends in Queensland and Victoria are given in Table 4. Beneficially from the GC approach, the ΔGWS trend is improved by approximately 20 % (from 0.4 to 0.6, compared to 1.6 cm yr⁻¹) in Queensland. The increase in ΔGWS is mainly influenced by the large amount of rainfall during the 2009–2012 La Niña episodes (see Fig. 9a). In Victoria, a significant improvement of the ΔGWS trend by about 76 % (from 0.1 to −0.2, compared to −0.3 cm yr⁻¹) is observed. Similar improvement of long-term trend estimates is seen in deseasonalized time series (improves by ∼ 15 % in Queensland and by ∼ 74 % in Victoria). The decrease in ΔGWS 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 decrease in the ΔGWS estimate cannot be properly captured in the model simulation. Applying GC approach effectively reduces the model deficiency and improves the quality of the groundwater estimations.

Figure 9The monthly time series of ΔGWS 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.

Download

Figure 10Mass changes (Gt, Gigatonne) of ΔTWS, ΔSM_top, ΔSM_rz, and ΔGWS estimated from GC solutions over 10 Australian basins in three different periods: Big Dry (January 2003–December 2009), Big Wet (January 2010–December 2012), and the entire period (January 2003–March 2016).

Download

5.4 Assessment of mass variation in the past 13 years

Australia experiences significant climate variability; for example, the Millennium Drought starting from late 1990 (Van Dijk et al., 2013) and extremely wet conditions during several La Niña episodes (Trenberth, 2012; Han, 2017). These periods are referred as “Big Dry” and “Big Wet” (Ummenhofer et al., 2011; Xie et al., 2016). To understand the total water storage (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 it is defined as the Big Wet period following Xie et al. (2016). In each period, the long-term trends of GC estimates of ΔTWS, ΔSM_top, ΔSM_rz, and ΔGWS are firstly calculated. Then, the total water storage variation (in m) is simply obtained by multiplying the long-term trend (in m yr⁻¹) 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⁻³). The estimated mass variations during Big Dry and Big Wet are displayed in Fig. 10. The long-term mass variation in the entire period (January 2003–March 2016) is also shown.

During Big Dry (2003–2009), a significant loss of total storage (40–60 Gt over 7 years) is observed over the LKE, MRD, NWP, and SWP basins. The largest groundwater loss of > 20 Gt 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 ΔSM_rz and ΔGWS explains more than 90 % of the ΔTWS value. The mass loss is also observed in ΔSM_top but is > 10 times smaller than ΔSM_rz and ΔGWS.

During Big Wet (2010–2012), the basins such as LKE, MRD, and TIM exhibit the significant total storage gain of > 100 Gt. The gain is particularly larger in ΔSM_rz over the basins that experienced the significant loss during Big Dry. For example, over LKE and MRD, the gain of ΔSM_rz is approximately 2–3 times greater than ΔGWS. 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 ΔSM_rz 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 ΔGWS compared to ΔSM_rz. In NEC, ΔGWS gain is ∼ 8 times larger than ΔSM_rz 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 ΔSM_rz loss observed over GOC is simply caused by the timing selection of the Big Wet period, which ends earlier (∼ 2011) in GOC than in other basins. The ΔSM_rz gain becomes ∼ 26 Gt if the Big Wet period is defined as 2008–2011. During the post-Big-Wet period (2012 and afterwards), the decreasing trend of water storage is observed from all basins (see Figs. 4 and 5). 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 ΔTWS 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 Gt during the 13-year period (see Fig. 10c). 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 ΔTWS loss of ∼ 79 Gt is observed over NWP, which is ∼ 25 Gt greater than the loss during Big Dry (see Fig. 10a and c). 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 ΔGWS, compared to the Millennium Drought period (see Fig. 10j and l). The continual decrease in water storage over western basins is likely caused by the interaction of complex climate patterns such as the El Niño–Southern Oscillation, Indian Ocean Dipole, and Southern Annular Mode cycles (Australian Bureau of Meteorology, 2012; Xie et al., 2016).

5.5 Comparison of the GC approach with alternatives

The simplest approach to estimate ΔGWS is to subtract the model soil moisture component from GRACE ΔTWS 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 in model outputs and GRACE data, the water storage states are updated through a Kalman filter:

where ${\stackrel{\mathrm{̃}}{\mathit{h}}}_{\mathrm{R}}$, H, C_R are described in Sect. 2; b is an observation vector containing GRACE-derived ΔTWS; and R is an error variance–covariance matrix of the observation. The GRACE-derived ΔTWS and its error information is obtained from the mascon solution. The matrix R is a (diagonal) error variance matrix since no covariance information is given in the mascon product. Note that the model uncertainty remains the same as in GC approach (Sect. 4.2). The different results from App 1 and App 2 are mainly attributed to the different estimates of the uncertainty.

Figure 11ΔGWS estimated from Approach 1 (App 1) and Approach 2 (App 2) in Queensland (a) and Victoria (b). The in situ groundwater network data and the GC solutions are also shown. Deseasonalized time series are shown in thick lines.

Download

The ΔGWS estimates from App 1, App 2, and GC in Queensland and Victoria are shown in Fig. 11. It is clearly seen that ΔGWS values from App 1 are overestimated while the one from App 2 fits the ground data significantly better. This behaviour was also seen in Tangdamrongsub et al. (2017), where the water storage estimates tend to be overestimated when error components such as spatial correlation error were neglected as in App 1. ΔGWS from App 2 shows clear improvements in terms of NS coefficients in both networks. Considering the deseasonalized ΔGWS estimates, in Queensland, the trend increases from 0.39 ± 0.03 to 0.42 ± 0.03 cm yr⁻¹ (improves by 1.5 %) and the NS value increases from 0.46 to 0.53. In Victoria, the trend decreases from 0.73 ± 0.10 to 0.46 ± 0.05 cm yr⁻¹ (improves by 27 %) and the NS value increases from −0.89 to 0.30. Although App 2 is not yet as good as the GC solution based on the most comprehensive error propagation, this simple test demonstrates the importance of considering the uncertainty. The reason for App 2 being less accurate than GC is likely the too-simplified error information implemented in App 2.