Assessing Changes in Terrestrial Water Storage Components over the Great Artesian Basin Using Satellite Observations

: The inﬂuence of climate change and anthropogenic activities (e.g., water withdrawals) on groundwater basins has gained attention recently across the globe. However, the understanding of hydrological stores (e.g., groundwater storage) in one of the largest and deepest artesian basins, the Great Artesian Basin (GAB) is limited due to the poor distribution of groundwater monitoring bores. In this study, Gravity Recovery and Climate Experiment (GRACE) satellite and ancillary data from observations and models (soil moisture, rainfall, and evapotranspiration (ET)) were used to assess changes in terrestrial water storage and groundwater storage (GWS) variations across the GAB and its sub-basins (Carpentaria, Surat, Western Eromanga, and Central Eromanga). Results show that there is strong relationship of GWS variation with rainfall (r = 0.9) and ET (r = 0.9 to 1) in the Surat and some parts of the Carpentaria sub-basin in the GAB (2002–2017). Using multi-variate methods, we found that variation in GWS is primarily driven by rainfall in the Carpentaria sub-basin. While changes in rainfall account for much of the observed spatio-temporal distribution of water storage changes in Carpentaria and some parts of the Surat sub-basin (r = 0.90 at 0–2 months lag), the relationship of GWS with rainfall and ET in Central Eromanga sub-basin (r = 0.10–0.30 at more than 12 months lag) suggest the effects of human water extraction in the GAB. Factors groundwater-dependent


Introduction
The Great Artesian Basin (GAB) is one of the world's most extensive artesian aquifer systems, underlying approximately 25% of Australia and containing approximately 65,000 km 3 of groundwater. It is a substantial water source for human needs, agriculture, and mining industries [1]. Groundwater discharges from the GAB sustain numerous spring wetlands, which have substantial ecological, scientific, and socio-economic significance [2]. However, the GAB has seen an overall decline in groundwater levels during the past century, exacerbated by human activity (e.g., mining), changing climate conditions [3][4][5], and extraction (e.g., through bore wells), with massive demand from the pastoral industry [3]. In a recent review of monitored groundwater flow and its underground vertical leakage in the GAB, Habermehl [6] observed that some artesian springs have dried up in highly developed regions as a result of up to 100 m reductions in artesian groundwater pressure. In addition, groundwater extraction across the GAB has resulted in decreasing groundwater levels and the drying up many springs [7][8][9].
The GAB spans a range of climates, from tropical, semi-arid and arid, and surface water bodies are largely non-perennial [10]. The scarcity of surface water in the GAB makes predominantly arid and semi-arid regions [3]. The groundwater temperature in the basin varies between 30 • to 100 • C [42]. Groundwater recharge occurs through precipitation in the GAB. Rainfall enters primarily along the uplifted eastern edge of the Great Dividing Range [43,44], with the western edge being a lower recharge zone due to negligible rainfall [45]. It follows that groundwater recharge zones are mostly located in the southeast and central part of the GAB (Figure 1a). Groundwater recharges occur through rainfall and discharge occur through evapotranspiration and extraction from the GAB. Groundwater within the GAB flows in a generally southwest direction, with some of the groundwater moving upward to the surface through various faults and fractures, creating ecologically important springs [46]. Since 1900, the groundwater level in the GAB has decreased [8] due to established mining companies in this area [47]. In addition, climate change (affecting rainfall and aquifer recharge rates) and human activities (e.g., groundwater extraction, mining) over time have an extensive influence on regional aquifer storage, groundwater levels, and surface discharge rates across the GAB [4]. The spatial distribution of average rainfall varies over the GAB between 2002 and 2017 indicates the complex hydrological processes in the GAB (Figure 1b).
The spatial distribution of ET loss (Figure 1c) gives spatial distribution of averaged ET in the GAB system and GLDAS Noah data was used to assess ET between 2002 and 2017. The ET loss is greater in the north and some parts of the southeast region of the GAB.

Datasets Used
We used four key datasets in this study: (1) terrestrial water storage (TWS) from GRACE satellite, (2) soil moisture storage from GLDAS Noah, (3) rainfall from gridded silos and (4) evapotranspiration from GLDAS Noah. The details of the datasets are outlined in Table 1 and Figure 2.

Gravity Recovery and Climate Experiment (GRACE) Data
In this study, we used terrestrial water storage derived from GRACE-mascons data and covered the period between the years 2002-2017. The data was retrieved from the Center for Space Research (CSR) (http://www2.csr.utexas.edu/grace/RL06_mascons.html, accessed on 15 January 2020) [48] and down sampled to a spatial resolution (0.5 • × 0.5 • ). Hydrological patterns (majorly TWS changes) have been obtained using especially GRACE data in other parts of the world [19,[49][50][51].

Global Land Data Assimilation System (GLDAS) Data
The Global Land Data Assimilation System (GLDAS) generates soil moisture time series and other parameters, such as evapotranspiration, canopy water evaporation, and surface runoff amount, which can improve the understanding of surface water extent [52]. The soil moisture component is an integral part of terrestrial water storage and aids in computing groundwater storage [53]. The GLDAS dataset is a terrestrial modelling system that integrates satellite and ground-based observations. The novelty of this dataset is that it separates the satellite and ground-based observations from the atmosphere [52]. In addition, GLDAS-based soil-moisture data provides an adjustable model time step and output interval and spatial resolution of 1 • × 1 • [54]. GLDAS soil-moisture dataset (with a spatial resolution of 0.25 • × 0.25 • ) from National Aeronautics and Space Administration (NASA's) Goddard Space Flight Centre (GSFC) portal (https://disc.gsfc.nasa.gov/datasets/ GLDAS_NOAH025_M_2.1/summary?keywords=GLDAS, accessed on 1 April 2020) [55] was used in this study to assess the soil moisture storage component over the GAB from 2002 to 2017.

Rainfall
We used the silos gridded rainfall data to characterize the spatio-temporal patterns of rainfall in the GAB. The silos dataset (Queensland Government) provides monthly gridded estimates of precipitation in high spatial resolution (0.05 • × 0.05 • ). The monthly rainfall datasets downloaded from the Queensland Government website portal (https: //www.longpaddock.qld.gov.au/silo/gridded-data/, accessed on 18 April 2020) [56] were used in this study. In addition, the rainfall data was used to determine its impact on GWS variation.

Evapotranspiration (ET)
Evapotranspiration is an essential component in water balance estimation over an area [57]. There are various techniques available to assess the ET component for a region [58], and in this study, GLDAS Noah satellite data from the NASA GSFC portal (https://disc.gsfc.nasa.gov/ datasets/GLDAS_NOAH025_M_2.1/summary?keywords=Evapotranspiration, accessed on 1 April 2020) [55] was used to extract monthly ET data for the GAB between 2002 and 2017. Further, ET data was used to understand the water loss across the basin and how it affects variation in GWS.

Groundwater Storage Variation Estimate from GRACE and GLDAS Satellite Data
GWS variation can be determined from TWS variation (GRACE satellite data) by removing other water storage components (e.g., soil moisture storage and surface water storage) from TWS. For semi-arid regions such as the GAB, the surface water storage component observed in the GRACE TWS data is likely to be negligible, as indicated by the lack of strong gravimetric fluctuation in surface water. Yan et al. [10] concluded that surface water and canopy storage variations contribute only 3% to observed terrestrial water storage over the GAB. If these water storage components are ignored following Yan et al. [10], groundwater storage can be estimated from terrestrial water storage as, where ∆WS ground represents changes in GWS, ∆TWS shows changes in TWS and ∆WS soil represents changes in soil moisture storage. Further, the surface water storage variations are negligible compared to soil moisture and terrestrial water storage variation in Australia [27].

Spatial-Temporal Patterns of Water Storage Components Using Principal Component Analysis
This study implemented the Principal Component Analysis (PCA) technique on rainfall, TWS and GWS datasets to summarize spatio-temporal variations in rainfall, TWS and GWS. PCA is a statistical decomposition method that decomposes multi-dimensional data and reduces its dimensionality and interpretability [59,60]. The usefulness of this analysis technique has gained popularity in atmospheric science and hydrological science for its dimensionality minimization and simple interpretation nature [61][62][63]. PCA transforms the dataset (e.g., TWS, GWS and rainfall) linearly and obtains a set of orthogonal vectors encompassing the very same region [60,64]. Mathematically, the eigenvalues and eigenvectors of a covariance matrix determine the principal components (PCs) of a given dataset [65]. This method helped in determining principal components (i.e., temporal variations) and empirical orthogonal functions (EOFs) (i.e., spatial maps). A scree plot analysis was employed to ensure that only significant orthogonal modes of variability were interpreted in all the hydrological units such as TWS, GWS and rainfall over the GAB [61]. The following equation was used to decompose variations in rainfall, TWS and GWS, where a (k) (t) represents temporal variations (also known as standardized scores) and p k are the corresponding spatial patterns (known as the empirical orthogonal functions [EOF]loadings). The standardized score is part of the total variation proportional to the total covariance in the time described by the eigenvector (EOF). EOFs have been normalized using the standard deviation of their corresponding principal components. For instance, while the EOF represents the spatial distribution of TWS, GWS or rainfall, the EOF/PC pairs are called PCA modes. In our study, PCA was employed to statistically decompose GRACE and rainfall datasets into PCs (temporal) and EOFs (spatial) to help in identifying the dominant patterns of GWS, TWS and rainfall in the GAB. Across the entire space-time dataset, 20 out of 183 months (10.9%) of total observations were missing over the 2002-2017 study period. These missing values occurred as random gaps in between years and were filled using linear interpolation, which is a common method to reconstruct or predict missing hydrological time series of this nature [27,59]. This interpolation did not impact on the overall data quality. With a consecutive monthly time-series of GRACE observations (183 time-steps starting from April 2002-June 2017) following the linear interpolation, we then implemented the PCA.

Time Series Analyses of Water Storage Components
Time series analysis of monthly averaged water storage components (TWS, GWS, ET and rainfall) was performed to determine the changes in these hydrological fluxes in time. In addition, time-series analyses were also executed to understand the variation and connectivity in different water storage components at each sub-basin (Carpentaria, Surat, Western Eromanga, and Central Eromanga) and for the entire GAB.

Average Annual Cycle and Deseasonalization of GWS and Rainfall
The average annual cycles of GWS and rainfall for each sub-basin in the GAB were assessed to investigate seasonal variation in GWS and rainfall. GWS variation and rainfall amplitudes were observed to find the peaks for respective months and helpful in finding the phase lag between rainfall and GWS variation. In addition to assessing the average annual cycles, deseasonalization of GWS and rainfall time series was undertaken to assess the human footprints and climatic contributions to the basin's hydrology.
Here, deseasonalizing means eliminating the annual seasonal signals from the GWS and rainfall time series.
To further understand the drivers of GWS variation across the basin, the GWS-rainfall relationship was analyzed using least square method. To this end, the spatial patterns of GWS variation and rainfall trends [66] between different time periods such as 2002-2008, 2009-2012 and then from 2012-2017 were executed using least square analysis.

Relationship between Water Storage Components (TWS, GWS, ET) and Rainfall Using Cross-Correlation Analyses
To further understand how the GWS varied in response to climate and non-climatic factors, cross-correlation analysis was used to evaluate the strength of agreement between rainfall, ET, and the water storage components TWS and GWS [59,61]. By implementing cross-correlation method on various water storage components such as rainfall, ET, GWS and TWS, the GAB sub-basins with maximum correlation can be obtained (i.e., between rainfall and GWS, rainfall and TWS, ET and GWS), to observe the GWS-rainfall relationship and lag months corresponding to the maximum correlation. The time-lagged correlation between the water storage components (e.g., GWS variation and rainfall) helped understand the response time (i.e., time taken by groundwater storage to recharge from rainfall). Spatio-temporal variations in groundwater storage relating to rainfall provided insight into climate variability and human intervention impacts in the GAB, as reported by some past studies [3][4][5]10,35,[37][38][39].

Understanding GWS Changes-Rainfall Relationship Using Multi-Linear Regression Analysis
Multiple linear regression was used to model the trend, annual and semi-annual components of variation in GWS and rainfall. Mathematically, it can be represented [18,59] as, where X(t) represents GWS or precipitation at time t, β0 represents constant, β1 shows linear trend, β2 and β3 shows periodic components (annual signal), while β4 and β5 show periodic components (semi-annual signals). The term β6 is the amplitude of GWS variation linked to climate signals. The variable E represents the normalized sequence of respective climate index (i.e., after eliminating long-term average) and ϕE shows the phase difference between GWS temporal variation and respective climate mode. ε is the error and shows the deviation between model outputs and observations. To validate the strength of MLRA, other parameters such as root mean square error (RMSE) and coefficient of determination (R 2 ) were assessed grid by grid. Using the multi-linear regression analysis (MLRA), trends and harmonic components of GWS changes such as rainfall patterns over the GAB areas were obtained [61,67]. This approach helps in assessing the GAB region in a more accurate way using RMSE and R 2 parameters.

Spatial and Temporal Changes in Rainfall across the GAB
The gridded rainfall data over the GAB for a 15-year window (i.e., 2002-2017) is analyzed using PCA. The mean monthly gridded rainfall pattern observed from rainfall data ranges from 5 to 180 mm with higher-level over the northern (i.e., in Carpentaria sub-basin) part and lesser over the western (Western Eromanga sub-basin) GAB (Figure 1b). The mean distribution of rainfall over GAB shows that much of it falls in the northern region compared with the typically dry southwest region of the basin (Figure 1b).
The total variability in rainfall observed by the first three principal component (PC) modes is 78.67% (Figure 3). In the first principal component, EOF loadings (or spatial patterns) are relatively strong and account for 63% of the cumulative variability in rainfall and are majorly localized over the Northern (Carpentaria sub-basin) part of the GAB (Figure 3a). TWS variations in the first principal mode are also localized in the northern region (Carpentaria), suggesting rainfall as a driver of water storage in the region. The second mode (Figure 3b) represents 9.81% of the cumulative variability. Figure 3b shows semi-annual rainfall distribution concentrated over the northern (Carpentaria sub-basin) part of the GAB. The third PC mode (Figure 3c) accounts for approximately 6% of the cumulative variability. The temporal patterns of rainfall associated with the spatial loadings (spatial patterns) in the northern region show strong annual variability (PC-1/EOF-1, Figure 3d), contrasting the multi-annual and seasonal variations observed in PCs 2 and 3 (Figure 3e,f). Changes in terrestrial water storage is analyzed using the PCA technique. From the scree plot analysis of terrestrial water storage (TWS), the first three PC modes accounted for 92.75% of the total variance and are chosen as the dominant patterns (or leading modes) of variability in GRACE hydrological signals over the GAB. The first PC mode (Figure 4a,d) of TWS shows relatively stronger spatial (EOF loadings) and temporal patterns in the northern region (Carpentaria sub-basin). Other orthogonal modes of variations in TWS (Figure 4b,c) show multi-annual variations (PC-2) and short-term GRACEhydrological signals (PC-3) and accounts for the remaining variability (28.79% and 3.26%, for the second and third modes, respectively). The GRACE-hydrological signals in these modes are associated with parts of Surat and Western Eromanga sub-basins (Figure 4b,c). As with temporal evolutions of TWS in the Carpentaria sub-basin, they both show pronounced peak amplitudes between 2010 and 2012 (Figure 4d-f). When the temporal and spatial patterns are jointly interpreted, changes in TWS show declines in the temporal series associated with parts of the Surat sub-basin since 2012 (Figure 4e). Looking closely at the Carpentaria sub-basin, a similar trend is noticeable but is gradual (Figure 4b).

Changes in Groundwater Storage
The spatio-temporal patterns in the first PC mode for groundwater storage variation reveal annual groundwater storage localized over the northern region of the basin. This is the dominant pattern of GWS accounting for approximately 59% of the total variability in the basin (Figure 5a,d). The second PC mode captures about 21% of the total variability in GWS and depicts relatively higher multiannual variation in the southeast region of GAB (Figure 5b). The corresponding temporal series (Figure 5e) associated with this spatial pattern shows peak amplitudes that coincide with the 'big wet period' and a steady decline after this period (2013)(2014)(2015)(2016)(2017). Figure 5c,f show a 7.11% total variability in the third PC mode and represents the multi-annual variations over the GAB. In this mode, the north and southwest part (i.e., in Carpentaria and Western Eromanga sub-basin) of the GAB depicts higher EOF loadings.   (Figure 6a). TWS and GWS variation obtained for the Carpentaria sub-basin show similarity with each other in terms of magnitude and time. However, there is some delay between GWS variation and rainfall (Figure 6b), which is discussed later in the manuscript. For the Surat sub-basin, water budget indicators (TWS, GWS, rainfall and ET) show temporal variations that suggest they are related to one another such that rainfall acts as a key input to the hydrological system and results in a consistent response of TWS followed by GWS (Figure 6c). Considering the Central Eromanga sub-basin, the variation in water storage also shows higher amplitude between 2010 and 2012 (Figure 6d), as observed for the Surat sub-basin (Figure 6c). GWS variation results for different GAB sub-basins and the whole of the GAB strongly represent the complexity of hydrological processes within the GAB. The fluctuations in rainfall in the Surat sub-basin, for instance, appear to be in the opposite phase with GWS, in contrast with the Carpentaria, where the temporal evolutions (patterns) of both data maintain the same phase, largely due to similar annual cycles. The Eromanga sub-basins (central and west) are typically dry, and the response of groundwater to periods with pronounced strong peak in annual rainfall (e.

Average Annual Cycles and Deseasonalization of GWS and Rainfall
Looking at the individual sub-basin responses, apart from the Carpentaria (Figure 7a), the seasonal annual cycles of GWS are in opposite phase with rainfall (Figure 7c,e,g). Deseasonalized GWS and rainfall time series are evident when annual seasonal cycles were removed from GWS signals. To track human footprints/factors other than rainfall in GWS variation, the long-term trend and seasonal signals are observed and GWS variation time-series with non-climatic factors is obtained. Here, long-term trend provides information on what influences water storage between 2002 and 2017 during water years (e.g., flood or other natural variability occurrence). From this, it can be deduced that an increase in rainfall will likely cause an increase in TWS or GWS and the lack of rainfall can cause a decline in GWS. It is likely that this increase or decrease in GWS is climateinduced (Figure 7b,d,f,h). However, GWS usage during the March-November period, when rainfall is largely limited, may have an impact on GWS variations (seasonal signal or human footprints). Isolation of GWS seasonal cycles from its averaged times-series can help in finding the inclusion of human-induced factors in GWS variation (Figure 7b,d,f,h). There are significant differences between non-deseasonalized and deseasonalized GWS signals for some of the sub-basins (indicated by the red color) where it is possible to tease apart some factors that are driving the GWS variation in the basin (Figure 7b,d,f,h). Observing the deseasonalized GWS time series for the Carpentaria, rainfall is strongly associated with GWS variation (Figure 7a). In other sub-basins, deseasonalized GWS time series shows that GWS variation is associated with factors other than rainfall (Figure 7d,f,h). Western Eromanga, an arid region in the GAB (Figure 1b), shows the existence of other non-climatic factors associated with varying GWS (Figure 7h). However, there is significant difference between non-deseasonalized and deseasonalized GWS time series in the Surat and Central Eromanga sub-basins depicting that GWS variation in these regions are likely driven by the influence of climatic and non-climatic factors (e.g., human extraction, industrial and agricultural use), especially in the Surat sub-basin (Figure 7d,f).
Isolating GWS seasonal cycles from the averaged GWS time series showed that GWS variation in the Carpentaria sub-basin has a strong GWS annual component and is largely driven by annual rainfall (Figure 7b). The Surat, Central Eromanga and Western Eromanga sub-basins show multi-seasonal oscillations that capture extreme events such as the big wet period during 2010-2012 period (Figure 7d,f,h). Regarding Figure 7d,f,h, it seems the declines during the 2012-2016 period appears to be sharper in the deseasonalized GWS. Overall, the steady declining trends resulting from a probable human groundwater use are evident in these regions and appear to be more pronounced during the post big wet period.

Trends in Ground Water Storage Variations
The spatial patterns of trends in GWS and rainfall over GAB reflect a complexity of geology and hydrological processes in the basin. We found that long-term GWS variation over the southeast region is inconsistent with rainfall variation (Figure 8a

Response of Land Water Storage to Climate Variability
Rainfall and evapotranspiration are major factors causing GWS variations [11]. For that reason, the response of TWS and GWS to rainfall and ET is assessed. Figure 9 represents the maximum correlation coefficients (r) value between the two variables (for example, GWS and rainfall) and the lags at which GWS and rainfall show maximum correlation (Figure 9a,b). From the observed r value, it is clear that rainfall drives GWS variation for more than 50% of the GAB. For example, GWS variation shows a relatively high correlation with rainfall in the northern and southeast regions of the GAB (Figure 9). This is also confirmed from the deseasonalized trends in the Carpentaria sub-basin (Figure 7a,b). The north and southeast parts of the GAB show that rainfall precedes GWS variation (lags ranging from approximately 2-12 months) with correlation coefficients ranging from 0.50 to 0.70 (Figure 9a,b). As can be seen from Figure 9b  There are some parts of the GAB, such as the southwest region, that are characterized by small rainfall amounts. The correlation between rainfall and GWS variation in the southwest region is lower than other regions, ranging from 0.20 to 0.30 (Figure 9a). The southwest region's rainfall is lagging behind GWS variation with a phase of approximately more than a year (Figure 9b). This is expected as the southwest region is characterized by a dry climate [5]. Relatively low values of correlation coefficients (r) are observed for the southwest region of between 0.10 and 0.30. Low values of r could imply weaker influence of rainfall over groundwater recharge. Overall, rainfall drives the hydrology (i.e., GWS, TWS) in the northern region (Figure 9a,c). Deseasonalized and average annual cycles of GWS and rainfall also validate these results for each of the sub-basin (Figure 7).
Rainfall has a strong association with variation in TWS over more than half of the GAB (Figure 9c). Rainfall commonly precedes TWS variations by 2-3 months and is consistent in the GAB except for the southwest region where rainfall is considerably limited (Figure 9d). A strong correlation exists between ET and GWS variation for most of the GAB except for some of the southwest regions (Figure 9e). The north and southeast regions of the GAB show a negative phase lag in 2 to 5 months between ET and GWS variation (ET leads GWS variation). The southwest regions show a positive phase lag between GWS variation and ET (GWS variation precedes ET) of 15 months (Figure 9f). This is consistent with the observed relationship between GWS variation and rainfall for the southwest region of the GAB (Figure 9a).

Understanding Drivers of Groundwater Variability
Multiple linear regression analysis (MLRA) is implemented to understand drivers of variations in GWS across the GAB. GWS and rainfall indicates consistency in average annual amplitudes (Figure 10a,f) and also coincides with the leading orthogonal modes of variability in both data (Figures 3 and 5). However, GWS in the north of Eromanga shows an annual signal (Figure 10a,b). GWS trends indicate complete dissimilarity with the spatial patterns of rainfall trends for Central and Western Eromanga sub-basins (Figure 10b,c,g,h). In contrast with rainfall, the only sub-basin where significant trend in GWS exist is the Surat sub-basin (Figure 10g,h). The disparity in GWS and rainfall trends (Figure 10c,h) observed in the Surat sub-basin possibly reflects the complexity of geology and climate across the GAB [3]. Root mean square error (RMSE) values in the southern regions (arid/semi-arid region) of the basin are low for GWS and rainfall and relatively higher in the northern regions (Figure 10d,i). Values of R 2 (i.e., GWS) ranging from 0.80 to 1.00 are observed over the northern region of basin (Figure 10e). Such higher values of R 2 indicate that GWS variation in the northern part of GAB basin is likely affected by rainfall.

Changes in Terrestrial Water Storage
TWS is an indicator of freshwater availability and plays a vital role in quantifying important hydrological changes [68,69]. The temporal variation of TWS and GWS indicate that the GAB and its sub-basins have experienced increases in TWS variation between 2009 and 2011 ( Figure 6). A strong La Niña occurred between 2010 and 2012 [70], causing flooding in the northeast and southeast regions of Queensland. This climatic event was likely responsible for increases in terrestrial water storage between 2010 and 2012 in the northeast and southeast regions of the GAB [70], suggesting rainfall to be a substantial driver of variation in TWS.
TWS variation in the GAB was largely related to climate factors such as rainfall and evapotranspiration. Rainfall and evapotranspiration values were high in the north and southeast regions of the GAB compared with the southwest region where they are relatively lower. Furthermore, cross-correlation between TWS and rainfall were high in the north and southeast regions and imply rainfall to be a significant factor behind variations in TWS for these regions. ET may also cause variation in TWS. As an outgoing flux and water budget indicator, ET can serve as a sink term that induces strong changes in hydrological conditions apart from rainfall. Given the semi-arid nature of the south-east region (Surat sub-basin), and the considerably strong correlations between GWS and rainfall, ET is likely to be crucial to the water balance of this region. All-in-all, the influence of ET on GWS is likely to be much greater in parts of the GAB where rainfall is low and ET is high.

Understanding Changes in Groundwater Storage over the GAB
Our study presents the spatiotemporal patterns of rainfall in the GAB and evaluates its impact on hydrological stores (TWS and groundwater). We observed that rainfall has a leading role in driving GWS variation in the northern region of GAB, and other non-climatic factors (e.g., human water extraction) have also been recognized as driving GWS changes in south regions of the GAB, especially southeast (Surat sub-basin) [35,39]. Factors causing variation in GWS may decrease the number of springs and impact groundwater-dependent ecosystems.
The GAB is a complex hydrological basin. Based on our deseasonalized GWS and rainfall results, GWS variation in the southeast (Surat) region and the northern (Carpentaria) region showed that climate conditions and likely other factors, such as anthropogenic groundwater extraction, are associated with changes in GWS (Figure 7). GWS variation, rainfall trends and the cross-correlation results confirm the observations made from deseasonalized GWS and rainfall patterns (Figures 8 and 9). The inconsistencies that we observed in RMSE values across the southeast region (Figure 10d,j) may be due to interannual changes in rainfall or the complex hydro-geology in the basin [71]. Low values of R 2 (coefficient of determination) over the southern region could suggest non-climatic influence or human intervention for these parts of the basin (Figure 10e,j) [61]. This is supported by the MLRA results for rainfall and GWS variation, which suggest influence of human use on GWS variation. For instance, pastoral industries, coal, and gas power generation plants exist in the Surat sub-basin and depend on groundwater resources [72]. These groundwater extractions consequently impact groundwater storage in the basin. Kent et al. [38] exhaustively analyzed the water extraction in the southeast GAB from 1900 to 2015 and helped finding the temporal patterns of groundwater decline in the region.
GWS variation showed negative trends in the 2005-2006 and 2008-2009 periods before a La Niña [70] event which caused widespread flooding in the northeast and southeast regions of Queensland in 2011 ( Figure 6). The sharp rise in GWS over this period is likely linked to this strong La Niña event (Figure 5a). High rainfall was associated with a sharp rise in GWS after 2009 (Figure 5d). Given this, rainfall is most likely an essential contributor to recharging groundwater storage. The cross-correlation results between GWS variation and rainfall (r = 0.5-0.7) highlight that rainfall is a main factor behind varying GWS for more than half of the GAB where rainfall precedes GWS variation by two to four months (Figure 9a,b).
We observed a relative temporal difference between TWS and GWS variations. There is a noticeable water loss factor (e.g., evapotranspiration) creating a relative difference between GWS and TWS variation. To understand the changes of TWS, ET over the GAB was evaluated between 2002 and 2017. Additionally, the GWS variation and ET were cross-correlated with each other (Figure 9e,f). MLRA results of rainfall and GWS confirm the results obtained by cross-correlation analysis and suggest that GWS variation is driven by climatic variability over tropical parts of the GAB (i.e., Carpentaria) and by the combination of climatic and non-climatic factors over arid and semi-arid GAB (i.e., Western Eromanga, Surat and Central Eromanga).
Changing climate conditions impact the world in terms of increased or decreased level of rainfall, water storage, water inundation and drought [27,[73][74][75][76][77][78]. In past years, research on the GAB has highlighted that varying climatic conditions and unsustainable use of groundwater (e.g., through boreholes) has likely resulted in the decline of groundwater stores and impacted sensitive ecosystems [6,8,10]. The water storage in semi-arid or arid regions largely depends on rainfall [79][80][81]. GWS variations in the basins are likely caused by climate conditions and events such as ENSO (El Nino Southern Oscillation) [82,83]. The GWS variation in the GAB is influenced by changes in annual and seasonal rainfall in the GAB. Recent studies observed a strong relationship between climate conditions and water storage variation [5,10]. Thus, quantifying climate variability factors such as rainfall effects on GWS variation across the GAB could be advantageous for groundwater management strategies by finding the locations with strong effects of rainfall on GWS and the locations deprived of rainfall.
We can better manage groundwater-dependent ecosystems once we have a better understanding of groundwater storage variation in the GAB. Springs are groundwaterdependent ecosystems and are a permanent water source that sustains wildlife and habitat including endemic vegetation and fauna [84]. However, pastoral and mining activities may have adverse impacts on many springs, as observed in southeast regions of the GAB [2]. Identifying groundwater storage variation hotspots help us understand how climate effects and human impact vary spatially and temporally across the GAB and is necessary for the better management of groundwater-dependent ecosystems including springs [85].

Conclusions
The main conclusions of the study are outlined below: a.
GWS varied most in Carpentaria sub-basin and some parts of the Surat sub-basin. The relatively high amount of rainfall within the Carpentaria and south-east region of the Surat sub-basin is the main factor driving the observed variability in GWS in these regions. Multi-annual variations in GWS are prominent in the Surat sub-basin and coincide with variation in rainfall. This amplification of GWS variation is observed as a result of changing climate conditions and potentially human water extraction in the Surat sub-basin. b.
GWS variation over the GAB between January 2009 and March 2012 period shows the Surat sub-basin had positive and strongest rise (40 mm/year) in GWS coinciding with a wet period following the Millennium drought. However, the Surat sub-basin showed negative trends before and after the 2009-2012 period. Overall, rainfall and GWS variation trends are consistent for the Carpentaria sub-basin and inconsistent for south regions of GAB, which indicates that other important drivers of variation in GWS exist apart from rainfall. c.
In general, the rainfall-GWS variation relationship indicates a time lag of two to three months for more than half of the GAB. The GWS-rainfall relationship is strong for more than half of the GAB; however, it is low for some parts of GAB (e.g., Western Eromanga and some parts of Surat sub-basin). Such regions show a phase lag of approximately 12 months and highlight the probable effects of non-climate factors on GWS variation in addition to climatic variation.
d. The ET correlates more with GWS variation than the rainfall in the Surat sub-basin (Figure 9a,e). This observation indicates that ET is an important factor in the recharge processes in this low rainfall region (Figure 1b). If rainfall declines in this region, then it could be a problem for recharge in the Surat sub-basin. e.
GWS variation in the southern regions of GAB (i.e., in Surat, Western Eromanga and Central Eromanga sub-basins) showed a weaker relationship with climate (i.e., rainfall). This relationship could potentially be due to the combined effects of human water extraction and complex hydro-geological processes.