Groundwater Hydrograph Decomposition With the HydroSight Model

Groundwater, the most important water resource and the largest distributed store of fresh water in the world, supports sustainability of groundwater-dependent ecosystems and resilient and sustainable economy of the future. However, groundwater level decline in many parts of world has occurred as a result of a combination of climate change, land cover change and groundwater abstraction from aquifers. This study investigates the determination of the contributions of these factors to the groundwater level changes with the HydroSight model. The unconfined superficial aquifer in the Gnangara region in Western Australia was used as a case study. It was found that rainfall dominates long-term (1992–2014) groundwater level changes and the contribution rate of rainfall reduced because the rainfall decreased over time. The mean rainfall contribution rate is 77% for climate and land cover analysis and 90% for climate and pumping analysis. Secondly, the increasing groundwater pumping activities had a significant influence on groundwater level and its mean contribution rate on groundwater level decline is -23%. The land cover changes had limited influence on long-term groundwater level changes and the contribution rate is stable over time with a mean of 2%. Results also showed spatial heterogeneity: the groundwater level changes were mainly influenced by rainfall and groundwater pumping in the southern study region, and the groundwater level changes were influenced by the combination of rainfall, land cover and groundwater pumping in the northern study region. This research will assist in developing a quantitative understanding of the influences of different factors on groundwater level changes in any aquifer in the world.


INTRODUCTION
Groundwater is an important natural resource that supplies water to humans (Döll, 2009), especially as the primary source of drinking water for over two billion people (Famiglietti, 2014). In arid and semiarid regions, groundwater is also used for agricultural irrigation. At the global level, 50% of the domestic water supply, 40% of the industrial water supply and 20% of the irrigation water supply originate from groundwater (Zektser and Lorne, 2004). However, global groundwater depletion has been increasing since the 1960 (Wada et al., 2010). Many countries and regions now face serious problems of excessive groundwater depletion, such as North Africa, North China, North America, the Middle East, South and Central Asia, and Australia (Konikow and Kendy, 2005). For example, groundwater depletion in the United States during 1900-2008 is estimated totals approximately 1,000 cubic kilometers (km 3 ) (Konikow, 2013). Groundwater depletion rate in North China based on GRACE was 2.2 ± 0.3 cm/yr from 2003 to 2010, which is equivalent to a volume of 8.3 ± 1.1 km 3 /yr (Feng et al., 2013). Groundwater is often poorly monitored and managed (Famiglietti, 2014). This has led to notable social and economic impacts (Gleeson et al., 2012). As a result, more efforts and attention are required to better understand and manage groundwater resources.
Continuous groundwater level decline in unconfined and confined aquifers over a long period is an important depletion indicator (Zektser and Lorne, 2004). Groundwater level decline can not only impose negative influences on local economic and social development, but could also impose significant influences on natural streamflow and groundwater-dependent ecosystems (Wada et al., 2010;Fu et al., 2019). It has been widely recognized that groundwater level fluctuation is influenced not only by natural processes, such as precipitation, evaporation and river water stages (Zhou et al., 2020), but also by anthropogenic activities, such as groundwater abstraction (Shapoori et al., 2015a) and land use and land cover change (Yue et al., 2016;Cheng et al., 2017;Abiye et al., 2018). Therefore, it is necessary to detect these drivers and decompose groundwater hydrographs into individual drivers to support groundwater resource management and regional development.
A well-known method is time series analysis which has been widely adopted in groundwater hydrology to interpolate, simulate and predict groundwater levels and further quantify the effects of various drivers of groundwater level fluctuation (Peterson and Western, 2014;Peterson and Western, 2018;Obergfell et al., 2019). This method is relatively simple and requires few parameters, and the model is easily constructed and produces reliable results. The Hydrograph Analysis: Rainfall And Time Trends (HARTT) model proposed by Ferdowsian et al. (2001), Ferdowsian et al. (2002) and Ferdowsian and Pannell (2009) is a representative time series model for groundwater level modelling. Another time series approach model is the transfer function noise (TFN) model developed by von Asmuth et al. (2002) and von Asmuth et al. (2008). This model is based on predefined impulse response functions and simulates groundwater head observations by weighting input historical forcing data and estimating the noise component. TFN models, in contrast to the HARRT model, do not require groundwater level observation data obtained at regular time steps, and they do not assume a stationary climate (Peterson and Western, 2014). Therefore, the TFN model is readily applied in groundwater hydrograph modelling because groundwater observation data are not always acquired at regular time steps. This method has been applied in many studies, including the estimation of groundwater recharge (Obergfell et al., 2019) and aquifer hydraulic properties (Shapoori et al., 2015b) and even the decomposition of the observed groundwater head into different hydrological stresses Shapoori et al., 2015c). In this study, the TFN model developed by Peterson and Western (2014) was adopted to separate the contributions of land cover change and groundwater pumping from that of rainfall variation to the observed groundwater level fluctuation. The Gnangara groundwater system is one of the most important groundwater sources in Western Australia and is vital to the local drinking water supply, supplying over 40% of Perth's drinking water each year (Merz, 2009), as well as supporting nationally significant groundwater-dependent ecosystems, such as lakes, wetlands, woodlands and cave ecosystems (Western Australian Planning Commission, 2001). The Gnangara groundwater resources are a key factor to achieve sustainable social and economic growth in the region, and they provide approximately 60% of the potable water supply to the city of Perth (with a population 1.5 million people) (Elmahdi and McFarlane, 2009). However, the groundwater levels in this region have been found to decline in the unconfined (superficial aquifer with groundwater moves slowly) and confined aquifers (Leederville with a maximum thickness of more than 600 m and Yarragadee aquifers with a maximum thickness of more than 2000 m), over the last 40 years due to a combination of rainfall decline, groundwater abstraction for public and private water supply purposes (Iftekhar and Fogarty, 2017), evapotranspiration, and interception by pine plantations (Bekesi et al., 2009;Strobach, 2013). It has been widely accepted that the sustainability of the Gnangara groundwater resources and groundwater-dependent ecosystems is greatly threatened by the continued decline in groundwater levels (Xu, 2008). However, the relative contributions of these factors to the groundwater level decline remain uncertain. Therefore, quantitative estimation of the impact of the individual drivers on groundwater level dynamics is necessary for scientific and effective groundwater resource management in the Gnangara region.
The HydroSight program employed in this study is a highly flexible statistical toolbox developed by Peterson and Western (2014). It integrates multiple models, such as soil moisture and time series models, into one model provides an efficient means to build a wide range of groundwater time-series models and separate the impacts of different factors from climate on groundwater level without knowledge of programming knowledge (http://peterson-tim-j.github.io/HydroSight/). Based on this tool, this study focuses on long-term groundwater level change analysis in unconfined superficial aquifer to identify the controlling factors and quantitatively decompose groundwater hydrograph into different drivers, such as rainfall variability, land cover change (represented by the normalized difference vegetation index (NDVI)) and groundwater pumping (for public and private water supply purposes), in the Gnangara region. This study has implications for groundwater resource management and regional development.
area of 2,200 km 2 and is bounded by the Gingin Brook and Moore River to the north, the Swan River to the south, the Darling Scarp, Ellen Brook and Swan Valley to the east, and the Indian Ocean to the west (Davidson, 1995;Western Australian Planning Commission, 2001). This region experiences a Mediterranean climate with hot, dry summers and mild, wet winters. The longterm (1900-2017) mean annual rainfall is 758 mm and the mean annual potential evapotranspiration reaches 1,386 mm (calculated based on the daily rainfall and potential evapotranspiration which is extracted from the SILO Data Drill database, see the detail in 2.2.1). Approximately 80% of the rainfall and 23% of the potential evapotranspiration are concentrated from May to September (Figure 2). A dry climate period following 1968 was evaluated by cumulative deviation from the mean rainfall (CDFM) technique ( Figure 2). Over the last 47 years , rainfall has declined by 13% below the long-term average , which has influenced the local ecosystem (Wilson et al., 2012).
The long-term mean annual rainfall was considered to determine the spatial distribution of rainfall in the Gnangara region via the empirical Bayesian kriging in ArcGIS 10.5 ( Figure 3A). Figure 3A shows that the rainfall is high in the south and central part of the Gnangara region (up to 865 mm) and is low in north part of the Gnangara region (as low as 688 mm).
In the Gnangara groundwater system, plantation forestry is one of the major land use types (land use map in 1992 and 2018 extracted from Bureau of Rural Sciences (2006) and Australian Bureau of Agricultural and Resource Economics and Sciences (2018), Supplementary Figure S1). Approximately 17,000 ha of the area contains currently mature maritime pine (Pinus pinaster) plantations, located in the centre of the system (Forest Products Commission, 2009). Mature pines consume much water because the transpiration rate of pines is 23% higher than that of the native banksia woodland surrounding the pine plantation (Carbon et al., 1982;Tremayne, 2010). Urbanized areas mainly occur in the southwest and southeast of the system. Other land use types, such as pastures, are found along the eastern and northern margins of the system (Australian Bureau of Agricultural and Resource Economics and Sciences, 2018). The long-term mean annual NDVI at various sites was adopted to determine the spatial distribution of the NDVI in the Gnangara region via the empirical Bayesian kriging in ArcGIS 10.5 ( Figure 3B) (with Root Mean Squared Error of 0.003). Figure 3B shows that high NDVI values largely occurred in the north and east of the Gnangara region, which are mainly covered with pine and banksia plantations and pastures. Low NDVI were mainly found in the south and middle west of the Gnangara region where urban areas are located.
The groundwater resources in the Gnangara region are mainly from winter rain. A high proportion of the precipitation infiltrates into the soil and recharges local groundwater due to the sandy soils with a high permeability occurring in the study area (Western Australian Planning Commission, 2001). The groundwater resources are usually abstracted for public water supply purposes, private water use, park and garden watering, industrial and commercial use, horticultural and agricultural irrigation, and domestic use. Approximately 42% of the extracted groundwater is used for the public water supply (Department of Water, 2009a). The Gnangara groundwater system comprises four different hydrogeological aquifers and the 3D Aquifer Visualization of the Gnangara region can be found in http://www.bom.gov.au/water/ groundwater/explorer/3d-aquifer-visual.shtml. The shallowest, unconfined superficial aquifer (its top surface is commonly termed the Gnangara Mound) which stretches across the coastal FIGURE 2 | Long-term averages of the monthly rainfall, annual rainfall and potential evapotranspiration, and dry climatic periods by cumulative deviation from mean rainfall (CDFM) in the Gnangara region.
Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 736400 plain, with an average thickness of 45 m, a maximum thickness of 75 m and a depth to groundwater ranging from 3 to 20 m. From east to west, the sediments of the superficial aquifer generally vary from being predominantly clayey adjacent to the Darling Fault and Gingin Scarp, to a sandy succession in the central coastal plain area, and to sand and limestone within the coastal belt. The hydraulic properties of the superficial aquifer vary significantly depending on geology. The horizontal hydraulic conductivity increased from western and eastern margin to central with values of 0.1, 15, and 50 m/day in Guildford Clay, Bassendean Sand, Tamala limestone. The shallow, semiconfined Mirrabooka aquifer which is mainly occurs in the southern and eastern regions of the Gnangara region and varies in thickness to a maximum thickness of 160 m, and its horizontal hydraulic conductivity of this aquifer ranges from 4 to 10 m/day. The deep, partially confined Leederville aquifer below the superficial aquifer, extends beneath the entire coastal plain except in the north near the Swan Estuary and in the south-east corner, and is typically several hundred meters thick, consisting of discontinuous interbedded sandstones, siltstones and shales with horizontal hydraulic conductivity of the sandstone beds about 10 m/day and that of the siltstone and shale beds about 1 × 10 -6 m/day. The Yarragadee aquifer is the deepest and major confined aquifer underlying the Perth Region and extending to the north and south within the Perth Basin. It is a multilayered aquifer often more than 2000 m thick, consisting of discontinuous interbedded sandstones, siltstones and shales with the average horizontal hydraulic conductivities range between 1 × 10 -6 and 10 m/day., and it offers a vast storage and a robust supply of groundwater (Davidson and Yu, 2008;Department of Water, 2009a;Environmental Protection Authority, 2017). The Gnangara Mound developed because the vertical rainfall infiltration rate (about 1 m/day) exceeds the horizontal groundwater flow rate (ranges from 50 m/year to 1000 m/year) in the aquifer (Davidson, 1995;Davidson and Yu, 2008). Groundwater recharge of the superficial aquifer is highly variable and depends on the local rainfall, land use, and geological conditions (Davidson and Yu, 2008;Department of Water, 2009a). The superficial aquifer is predominantly recharged by rainfall in winter with some upward recharge, occurring from the underlying Leederville and Yarragadee aquifers (Department of Water, 2009a). Groundwater is naturally discharged into wetlands, rivers, springs and ocean, and groundwater undergoes vegetation transpiration and leaks into underlying aquifers during groundwater movement. Additional discharge is associated with groundwater abstraction (Davidson and Yu, 2008;Department of Water, 2009a).

Climate Data
Climate data , including the daily precipitation and FAO56 potential evapotranspiration (FAO Penman-Monteith equation, Allen et al. (1998)), were extracted from the SILO Data Drill database created by the Queensland Government Department of Science, Information Technology and Innovation (DSITI). This dataset provides 0.05°gridded daily data and are constructed by spatially interpolating the observational data with methods of a thin plate smoothing spline and ordinary kriging (Jeffrey et al., 2001) and can be accessed on the Internet at https://www.longpaddock.qld.gov. au/silo/. Minimum and maximum temperatures, solar radiation and vapor pressure prior to 1956 were interpolated by an anomaly interpolation technique (Zajaczkowski and Jeffrey, 2020). At each observation site, the acquired climate data exhibit a daily time-step and start at 20 years prior to the first observation date of the groundwater level.

Land Cover Data
The NDVI (Normalized Difference Vegetation Index) -high resolution gridded (0.05°* 0.05°grid) monthly NDVI dataset (1992-2017) used as the land cover change data in this study was obtained from the Australian Bureau of Meteorology website (http://www.bom.gov.au/metadata/catalogue/view/ ANZCW0503900404.shtml). The satellite data originated from the Advanced Very High-Resolution Radiometer (AVHRR) instruments onboard the National Oceanic and Atmospheric Administration (NOAA) series of satellites operated by the US (http://noaasis.noaa.gov/NOAASIS/ml/avhrr.html). The NOAA-11, -14, -16 and -18 satellites were considered. The NDVI was used to examine the vegetation cover changes because the vegetation and non-vegetation area is easy to identify and historical land use information cannot be obtained. The NDVI value ranges from −1 to +1, where positive values indicate vegetation features and negative values indicate non-vegetation features (Gandhi et al., 2015).

Groundwater Monitoring Data
Groundwater monitoring data were retrieved from 325 observation bores in the unconfined superficial aquifer within the Gnangara region ( Figure 1). The observation records are irregular and the data from 1992 to 2014 are used in this study. This dataset was provided by the Department of Water and Environmental Regulation, Government of Western Australia. The details of cite description can be downloaded for free from https://water.wa.gov.au/maps-and-data/monitoring/waterinformation-reporting.

Groundwater Abstraction
A groundwater abstraction dataset (2,324 pumping sites) collected from the different aquifers (mostly pertaining to from superficial aquifer) in the Gnangara region, from 1992-2014, was used in this study. These groundwater pumping data referred to the extraction for public and private water supply purposes. The groundwater abstraction data were obtained from the Department of Water and Environmental Regulation, Government of Western Australia.

The Transfer Function Noise Model
The original TFN model was developed by von Asmuth et al. (2002) to simulate the groundwater level. This model includes three components: a deterministic component simulating the groundwater level due to the combined effect of all external factors, a residual series of the groundwater level and a constant component of the local drainage level. The Pearson type III distribution function was introduced by von Asmuth et al. (2002) to establish the precipitation and evapotranspiration impulse response functions. A revised version of the Pearson type III distribution function was adopted and five weaknesses in use for climatic stressors had been improved by Peterson and Western (2014). The first modification was to minimize the parameter covariance. The second modification was the calibration reproducibility was improved by undertaking a log 10 transformation of the parameters. The third modification was to reduce the impulse response function value at the start of the climate record. The fourth modification was to address the integration of the function from the first climate observation to negative infinity. The fifth modification was minimizing rounding errors in the numerical estimation of the integrals. Finally, Eq. 1 details the linear TFN model comprising two components of precipitation and evaporation: However, the linear TFN model does not simulate the groundwater head well because of the nonlinearity between precipitation and groundwater head (Peterson and Western, 2014). Therefore, a vertically lumped soil moisture model (Eq. 3) was introduced into Eq. 1. The model is highly flexible and contains one to five parameters (KlemeŠ, 1986). transformed parameter related to the maximum vertical soil saturated conductivity k sat [LT −1 ], E t [LT −1 ] is the potential evapotranspiration rate at time t, α is a dimensionless parameter controlling the fraction of precipitation available for infiltration as the catchment wetness, β is the antilog (log10) transformed dimensionless parameter controlling the free drainage as the catchment wetness, and c is a dimensionless parameter controlling soil evapotranspiration.
In the process of transforming the precipitation and potential evapotranspiration time series data into the required TFN model input, the precipitation P τ in Eq. 1 can be replaced with the free , and the potential evapotranspiration E t can be replaced with the soil functions θ P and θ E for precipitation and evapotranspiration, respectively, were adopted in the time series model. To consider for the groundwater level response to groundwater pumping, Shapoori et al. (2015a) and Shapoori et al. (2015c) added a third integral to the time series model (Eq. 4). Shapoori et al. (2015c) showed that models with and without groundwater evaporation produced similar results. Therefore, the second integral of groundwater evaporation was omitted in the models of this study.
where Q τ [L 3 T −1 ] is the daily pumping rate. The response function θ H F can either be the Hantush's equation for a leaky aquifer (Hantush, 1956;Eq. 5) or the Ferris and Knowles' well equation for a nonleaky aquifer (Ferris and Knowles, 1963; Eq. 6). The Ferris and Knowles' well equation without groundwater potential evapotranspiration was selected in this study due to the groundwater levels in most bores occurring deeper than 1 m according to tests conducted by Shapoori et al. (2015a): where a, b and c are parameters that only have a physical meaning if the basic Hantush assumptions are satisfied. According to Peterson and Western (2014), there are 84 nonlinear TFN models based on 16 soil moisture models available within HydroSight. Eqs 3, 7 (Peterson and Western, 2014), Eq. 8 (Peterson and Western, 2014), and Eq. 9 (Siriwardena et al., 2011) are four model structures used in this study and is named as structure "cccc," "c1c1," "c0cc," "inf101." Structure cccc means all parameters are calibrated. Structure c1c1 means fixing α c 1. Structure c0cc means fixing α 1. Structure inf1c1 means fixing α c 1 and β 0. Free drainage was chosen to transform the precipitation time series data into the required TFN model input. The Ferris and Knowles' well equation was selected as the response function.

Model Calibration and Evaluation
Following Shapoori et al. (2015a) and Shapoori et al. (2015c), the first 70% and last 30% of the observation records were designated as the calibration and evaluation periods, respectively. To identify individual parameter which produces the best possible fit to the hydrograph, the default calibration methods in the toolbox of Shuffled Complex Evolution with Principal Components Analysis -the University of California at Irvine (SP-UCI), developed by Chu et al. (2011), were used. SP-UCI is a global optimization algorithm based on the Shuffled Complex Evolution (SCE-UA) method (Duan et al., 1992) to address highdimensional and complex problems.

Model Performance Assessment
The Akaike information criterion with correction (AICc) and Bayesian information criterion (BIC) (Burnham and Anderson, 2004) were used to account for the number of model parameters.
The lower AICc and BIC indicates the better model. Moreover, the Nash-Sutcliffe efficiency (NSE) has been widely used to assess the goodness fit of a hydrograph (Ritter and Muñoz-Carpena, 2013;van der Spek and Bakker, 2017). Therefore, to assess the TFN model performance, the NSE (Nash and Sutcliffe, 1970) and unbiased NSE (Shapoori et al., 2015a) were adopted in the calibration and evaluation periods, respectively. NSE ranges from −∞ to 1. The model is more accurate when the NSE value is closer to 1. At NSE 1, the modelled data are a perfect match to the observed data (only if the measurements are free of errors). At NSE 0, the modelled data are as accurate as the mean of the observed, and if NSE <0, the modelled data are less accurate than the mean observed data.

Radius of Influence
Based on the current data on the study area, the Dupuit equation (Dupuit, 1863) was selected to calculate the radius of influence of pumping wells. In the confined aquifer, the Dupuit equation (Eq. 10) was used to calculate the radius of influence of pumping wells.
In the unconfined aquifer, the Dupit equation (Eq. 11) (Dupuit, 1863) was used to calculate the radius of influence of pumping wells.
r 10 s 1 log 10 r 2 −s 2 log 10 r 1 r 10 where r [m] is the radius of influence of a pumping bore, s 1 and s 2 [m] are the drawdowns in the observation bores, r 1 and r 2 [m] are the distances between the observation and pumping bores, and H [m] is the thickness of the aquifer. In this study, the thickness of the unconfined aquifer was determined according to Smith and Pollock (2010).

Model Structure Comparison
To determine the most applicable model structure, four soil moisture model structures (structures cccc, c1c1, cocc and inf101) were tested in all bores for climate-only analysis. For AICc and BIC, there are not obvious differences between these four model structures, but structure cccc showed worse than other three model structures ( Figure 4A). The climate-only analysis means the groundwater level is considered only influenced by rainfall. The NSE values of 325 bores using these four soil moisture models are listed in a box plot is shown in Figure 4B. During the calibration period, model structure inf101 performed the worst, with the lowest mean, median values of the NSE. The model structure cccc performed the best with the highest mean value of NSE in calibration but performed the worst with the lowest mean value of NSE. The model structure c0cc performed as good as model structure c1c1. Therefore, either model structure c1c1 with less parameters than c0cc was appropriate in this study. In this study, model structure c1c1 (Eq. 7) was chosen to model all the bores in the Gnangara region in all subsequent groundwater level analyses.

Climate Only
Sixty percent of the bores in climate-only (C) analysis during the calibration period produced an acceptable performance (NSE >0.5), of which 36% of the bores performed very good, 10% of the bores performed good and 14% of the bores performed satisfactorily. In addition, 40% of the bores performed unsatisfactorily. During the evaluation period, 48 Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 736400 and 52% of total bores yielded an acceptable and unsatisfactory groundwater head modelling performance, respectively. Spatially, bores producing an acceptable model performance in climate-only analysis predominantly occurred in the southeastern Gnangara, followed by the northern Gnangara, and an unsatisfactory performance was predominantly produced in central-western and northwestern Gnangara and the coastal area of the Gnangara region during the calibration period (Supplementary Figure S2). The spatial distribution of model performance during the evaluation period is similar with that during the calibration period. However, an unsatisfactory performance increased in south area of the Gnangara region (Supplementary Figure S2).

Climate and Land Cover
The model performance was improved during both the calibration and evaluation periods when the NDVI was added to the TFN model. The acceptable performance increased from 60 to 62%, especially the good and satisfactory performance level with an improvement of 2% during the calibration period. In addition, the unsatisfactory performance level was reduced by 2% during the calibration period. However, the acceptable performance level decreased from 48 to 47% during the evaluation period.
Among the 325 bores in the Gnangara region, an acceptable performance was largely attained in the south, southeast and north of the Gnangara region during the calibration and evaluation periods in the climate and land cover (C + L) analysis (Supplementary Figure S2). An unsatisfactory performance during the calibration and evaluation periods was predominantly distributed is mainly attained in the centralwestern and northwest of the Gnangara region (Supplementary Figure S2). Figure 5 shows the spatial distribution and magnitude of the NSE improvement considering the land cover data were considered. The improvement was calculated by subtracting the NSE value of the climate-only analysis from the NSE value of the climate and land cover analysis. The model performance had decreased when the NSE change value was smaller than −0.1. The model performance remained stable when the NSE change value varied between −0.1 and 0.1. The model performance had improved when the NSE change value was larger than 0.1. During the calibration period, the NSE value at 17 bores had improved with the NSE change values ranging from 0.1 to 1.85 and 18 bores had been reduced by values ranging from −8.62 to −0.1. Most of the bores (290) remained stable, with NSE change values ranging from -0.1 to 0.1. During the evaluation period, the NSE value at 21 bores had improved with NSE change values ranging from 0.1 to 0.85, and at 19 bores, the performance had decreased with NSE change values ranging from -0.52 to −0.1. Most of the bores (285) remained stable, with NSE change values ranging from −0.1 to 0.1. In calibration and evaluation period, the NSE improved mainly at bores in northern Gnangara and the NSE declined mainly at bores in the northern parts of the Gnangara region.

Climate and Groundwater Pumping
According to the empirical Dupuit equation, the radius of influence of the pumping wells in the confined and unconfined aquifers ranged from 708 to 9,303 m, and 1,302 pumping sites were determined to influence 271 observation sites. At these 271 observation sites, the groundwater pumping stressor was included to assess its contribution to the groundwater level changes. The corresponding modelling results of the C analysis and C + L analyses were separated and compared to climate and pumping analysis (C + P) results during the calibration and evaluation periods.
To validate the pumping well influence on the groundwater head, models with multiple bores and only the nearest bore were applied. The NSE of these two models are shown in Figure 6. During both the calibration and evaluation periods, the models containing one bore and multiple bores did not exhibit notable differences, but the model with one bore attained a slightly better performance than the model with multiple bores. Therefore, the model with one bore was applied in the climate and pumping stressor analysis.
During the calibration period, 60% of the bores produced an acceptable modelling performance, and 40% of the bores produced an unsatisfactory modelling performance in the C analysis ( Figure 7). The acceptable model performance level had been improved by 3 and 19% in the C + L and C + P analyses, respectively. Especially for the very good performance rate had improved by 21% in the C + P analyses (Figure 7). Moreover, unsatisfactory model performance rate has been reduced by 2 and 18% for C + L and C + P analysis. During the evaluation period, the acceptable model performance rate increased slightly from 48% in the C analysis to 50% in the C + P analysis. In the C + L analyses, the acceptable model performance rate was even reduced by 3% (Figure 7). In general, the pumping analysis results were slightly better than the land cover analysis results during both the calibration and evaluation periods.
Among the 271 bores in the Gnangara region, an acceptable and unsatisfactory model performance predominantly occurred in the south and north (coastal region), respectively, of the study area in the C analysis (Supplementary Figure S3). In central and north area, the model performance level increased when the pumping stressor were considered (Supplementary Figure  S3). The model performance level was improved at a few bores in north area for C + L analysis. During the evaluation period, the unsatisfactory performance level was much higher that during the calibration period, especially in the bores of central, coastal and north part (Supplementary Figure S3). Figure 8 shows the spatial distribution and magnitude of the NSE improvement when the land cover and pumping data were considered. The improvement was calculated by subtracting the NSE value of the climate-only analysis from the NSE value obtained when the land cover and pumping data were included in the model, respectively. During the calibration period, most of the bores (140) remained stable, with NSE values ranging from −0.1 to 0.1. At 18 bores, the performance improved with NSE change values ranging from 0.1 to 0.85, and at 17 bores, the performance decreased with NSE change values ranging from −0.52 to −0.1 in the C + L analyses. During the evaluation period, 15 bores with an increased NSE value were found, with the NSE improvement ranging from 0.1 to 1.85, while 16 bores with a reduced NSE value were found with NSE reduction ranging from −0.1 to -8.62 during the calibration period. And 240 bores keep a stable NSE values. In the C + P analysis, the number of bores with an increased NSE in both the calibration (98, range: 0.1-0.98) and evaluation periods (59, range: 0.1-27.27) was larger than that in the C + L analysis. Moreover, the number of bores with an NSE reduction increased to 59 during the evaluation periods (range: 9.82∼−0.1) compared to that during the calibration (15 bores) period. The bores with an NSE improvement in the C + L analysis were largely distributed in the upper part of the bores area during the calibration and evaluation periods (Figure 8). The bores with an NSE reduction and increase in the C + L analysis were mainly foundin north area of the Gnangara region. In the C + P analysis, the areas in the southwestern Gnangara with a high concentration of pumping wells, the model performance level did not exhibit a notable improvement. Areas with an NSE

Contribution in Space
To determine the spatial distribution of the contribution of each stressor to the groundwater level changes, the multiyear average rainfall, land cover, and pumping contribution rate of bores with NSE exceeding 0.5 (acceptable model performance) during both the calibration and evaluation periods were used to generate the spatial distribution of the contribution rate in the Gnangara region by the Inverse Distance Weighted method in ArcGIS 10.1 (Figures 9, 10). As shown in Figure 9, the multi-year mean contributions of rainfall and land cover to the groundwater level changes ranged from 44 to 99% and from -53 to 43%, respectively. A large rainfall contribution rate was mainly found in central and southern Gnangara region. In contrast to the rainfall recharge contribution, the land cover contribution can either be positive or negative. A positive land cover contribution (green area) was largely distributed in upper central and southeastern Gnangara region. In certain areas, such as the southwestern and centralwestern areas (light yellow area), the land cover contribution was very small. In the south-western, north margin, central-eastern parts (red areas) of Gnangara region, the contribution was negative. Figure 10, the multi-year mean contributions of rainfall and pumping to the groundwater level changes ranged from 1 to 99% and from −1% to −99%, respectively. Large rainfall recharge contribution to the groundwater level changes mainly occurred in the southwestern and northwestern parts of the Gnangara region. In contrast to the rainfall recharge contribution, the pumping contribution to the groundwater level changes was negative. The bores most affected by pumping occurred in the north-eastern parts of the Gnangara region.

Contribution in Time
Over time, the multi-site mean contributions of climate change and land cover on groundwater level were 90 and 2% with ranges of 89-90% and 2.1-2.3%, respectively. The multi-site mean contributions of climate change and pumping were 77% and -23% with ranges of 68-99% and −1%∼−32%, respectively. In the whole Gnangara region, the trend of the rainfall contribution slightly decreased from 1992 to 2014 ( Figure 11A). The land cover contribution rate was stable from 1992 to 2014. The pumping contribution rate decreased from 1992 to 2014. However, the pumping contribution rate is negative, that means the pumping activities lowers the groundwater decline. Moreover, the trend of the contribution of factors on groundwater level were calculated. Seventy-four percent of the FIGURE 9 | Spatial distribution of the rainfall and land cover (as suggested by the NDVI) contributions to the groundwater level changes in the Gnangara region.
Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 736400 FIGURE 10 | Spatial distribution of the rainfall and pumping contributions to the groundwater level changes in the Gnangara region.
FIGURE 11 | Contribution rates in the climate-only (C), climate and land cover (C + L), and climate and pumping (C + P) analyses to the groundwater level changes (A). Rainfall, NDVI and groundwater pumping in the Gnangara region (B).
Frontiers in Environmental Science | www.frontiersin.org July 2021 | Volume 9 | Article 736400 13 rainfall contribution rate trend showed decline (< −0.01) over time and 26% of the rainfall contribution rate trend showed stable (0) for C + L analysis. Eighteen percent of the land cover contribution rate trend showed decline over time, 70% of the land cover contribution rate trend showed stable, and 12% of the land cover contribution rate trend showed increase (≥0.01) for C + L analysis. Sixty percent of the rainfall contribution rate trend showed decline over time and 40% of the rainfall contribution rate trend showed stable for C + P analysis. Fifty-eight percent of the pumping contribution rate trend showed decline over time, 42% of the pumping contribution rate trend showed stable for C + P analysis.

Model Parameters
The variability of critical fitted parameters of the model was presented in a box plot ( Figure 12) and the spatial distribution of the fitting parameters was presented in Figure 13. Figure 12 showed that the Ksat, β, and Scap between C + L and C + P analysis is slight. Parameter of Ksat (maximum vertical conductivity) ranged from 0.01 to 10 m/day with mean of 0.47 and 0.61 m/day for C + L and C + P analysis, respectively. Parameter of β(power term for drainage rate) ranged from 1 to 10 with mean of 3.59 and 3.56 for C + L and C + P analysis, respectively. Parameter of Scap (soil moisture storage capacity) ranged from 0.01 to 1 m with mean of 0.18 and 0.21 m for C + L and C + P analysis, respectively. Spatially, Figure 13 showed about 60% bores has the Ksat of 0.01-0.2 m/day for both C + L and C + P analysis. And larger Ksat (more than 2 m/day) was found in the northern Gnangara for C + L analysis and in both southern and northern Gnangara for C + P analysis. Parameter of β in spatial showed slight difference between C + L and C + P analysis and the larger β is mainly occurred in the north area of the Gnangara region. 84 and 76% bores has the Scap of 0.01-0.2 m. The bores with Scap larger than 0.4 m is dispersedly distributed in the study area.

Groundwater Hydrograph Decomposition Evaluation
The spatial distribution of the acceptable model performance level (in the southern Gnangara) in the climate-only analysis was consistent with the area with a high rainfall. Additionally, the areas with a large rainfall recharge contribution were largely found in the area with a high rainfall (Figure 3), while a small rainfall recharge contribution occurred in the area with a low rainfall. The rainfall always keeps higher contribution (47-99 and 1%-99%) than land cover (−53-43%) and groundwater pumping (−1% to −99%) contribution after the NDVI and groundwater pumping were included into the model, respectively. Generally, the area rainfall more than 770 mm have the contribution rate more than 75%. These results indicated that the spatial distribution of the rainfall recharge contribution to groundwater was consistent with the rainfall distribution based on a comparison of the rainfall recharge contribution graph (Figures 9, 10) and rainfall graph were compared (Figure 3). Temporally, the rainfall contribution rate continually decreased and were higher than land cover and pumping contribution rate ( Figure 11A), which is closely related to the rainfall reduction over time ( Figure 11B). All of the above results indicate that the climate (precipitation) is the main factor influencing the groundwater level in the Gnangara region (dominating the groundwater level changes). Yesertener (2005) and Gallardo (2013) also proposed that climate change is the primary factor influencing the groundwater level decline in the Gnangara region based on the cumulative deviation from main rainfall (CDFM) analysis method. This finding has been verified, and the contribution of climate change to the groundwater level changes in time and space have been quantified. The model performance was improved at some bores in northern Gnangara when the land cover data were added to the model, indicating that the land cover change also exerts an impact on groundwater level changes. However, the land cover contribution on groundwater level was stable over time and the contribution rate kept at 2%. Although the land cover had changed in some area over time, the NDVI changes for the whole region is little ( Figure 11B). This demonstrated that the land cover changes had influence on local long-term groundwater level changes, but the influence is limited and is stable over time. The locations of the sites with an improved model performance generally agreed with the vegetated area in the northern Gnangara region (suggested by high NDVI) ( Figure 3B). In the central of the northern Gnangara region, the land cover contribution rate was higher and the NDVI was lower than surrounding areas, indicating a negative relationship between land cover changes and groundwater level changes. These indicated that the groundwater level in the northern Gnangara region was closely related to the local vegetation conditions. Although the land cover contribution rate is positive in central of the northern Gnangara region, the positive trend over time was decreasing because of the increased areas of conservation and  Economics and Sciences, 2018). However, in the plantation forests area in land use map 2018, the land cover contribution rate is up to −53%. This is likely closely associated with the high density of pine plantations in these areas (Yesertener, 2007;Gallardo, 2013). In the southwestern area of the Gnangara region, urban use is the main land use, but the contribution is negative, that means other factors such as groundwater abstraction influenced the groundwater level.
The model performance was significantly improved after the pumping data were added to the model. In the south area of the study area, the model performance remained high. In the western part of the study area, which is the coastal region, the model performance remained unsatisfactory, although the land cover or pumping data were included. In the coastal areas, the climate, land cover and pumping exert limited impacts on the groundwater level changes, and other factors such as seawater intrusion, may impose major impact on the groundwater level changes. The limited changes in the groundwater level along the coast also occur due to the extremely high hydraulic conductivity of the Tamala limestone and proximity to the discharge zone (Costall et al., 2020). Other possibility of limited groundwater level changes may be from the boundary conditions (Morgan et al., 2012). As the C, C + L, and C + P analysis results indicated, the model performance was improved when the pumping data were considered, especially during the calibration period. However, the model performance was slightly improved when the land cover data were considered, and the performance at some sites was reduced. Moreover, the large pumping contribution (negative) to the groundwater level changes occurred in the northeast of the Gnangara region and did not always occur near the pumping sites, indicating that the distance between the bores and pumping wells and the density of the pumping wells are not critical elements in the determination of the pumping contribution to the groundwater level changes. In the southwestern Gnangara, the contribution rate mainly ranged between −15% and −40% ( Figure 10). However, the groundwater pumping contribution continually increased over time due to the sustained groundwater abstraction ( Figure 11B). It is obvious that groundwater pumping is highest in the 2011, and the pumping contribution on groundwater level is also highest with contribution value of −32% ( Figure 11B). Results indicated that sustained groundwater pumping over time has a significant influence on groundwater level decline. Depletion of groundwater levels is a global phenomenon and is defined as long term water level declination caused by sustained groundwater pumping over time (Tularam and Krishna, 2009). More attention should be paid to manage the groundwater pumping activities to ensure the sustainable use of groundwater.
The fitted parameter of the maximum vertical conductivity of the model in this study has a value range 0.01-9.8 m/day and the maximum vertical conductivity of most of the bores is concentrated in range of 0.01-2 m/day. Salama et al. (2005) and Department of Water (2009b) calculated the hydraulic conductivity in the Gnangara, ranging from 0.56 to 6.38 m/day and from 0.01 to 5 m/day, respectively. The similar results indicated the results of fitted parameters of the maximum vertical conductivity is accepted and the results of the model are credible. Moreover, the fitted parameter of Ksat and Scap showed slight difference in C + L and C + P analysis, indicating that the maximum vertical conductivity and soil moisture storge capacity are important parameter of the model.

Uncertainties and Limitations
The approach adopted in this study is based on various assumptions and has limitations, which constrain the results of the groundwater hydrograph decomposition. It should be noted that this method assumes that the drivers (climate, land cover, and groundwater pumping) of the groundwater level changes are independent of each other, and the land cover and pumping were added into the model, respectively. However, the considered drivers interact with each other, which may complicate the evaluation results. Climate change (rainfall change) could lead to changes in vegetation, thus affecting groundwater recharge, and groundwater pumping may lead to groundwater level changes, thereby negatively impacting vegetation (Şen, 2015). Many factors influence groundwater level changes, but only three major drivers were considered to establish the model in this study. Other drivers such as bush fires, pine clearing, and groundwater evaporation, were not examined. And the surface water-groundwater interaction is also an important cause of groundwater level changes in many areas (Wang et al., 2014). Moreover, this model assumed that the aquifer is homogeneous, but the thickness, permeability and water-bearing structure of the aquifer affect the accuracy of the modelling results.
In this study, model was constructed for each individual observation bore of the Gnangara region. Then, for the whole Gnangara region, the groundwater level changes and its main factors can be assessed by interpolating the results of each bore. In the model construction, an optimal model structure fitted for the study area was used for each model to reduce the model running time and improve the model operation efficiency. In fact, one model structure cannot ensure the satisfactory model results of each site. Therefore, model performance in spatial is not always satisfactory. However, in the process of interpolation, only the satisfactory model performance was used to reduce the errors in spatial.
Other limitations arise from the dataset used. NDVI was selected as the representation of the land cover impact on the groundwater level changes. The NDVI may not be accurate enough to express all of the land cover situation. However, the vegetation and non-vegetation area is easy to identify and historical land use information cannot be obtained. Furthermore, the groundwater level data at certain sites were fragmented or the number of observations was insufficient, but the established models allowed the simulation of irregular water level observations.

CONCLUSION
Understanding and interpreting changes in groundwater level is essential for long term management of both a groundwater resource and urban development. The present study quantitatively clarifies the impacts of three drivers, namely, climate change, land cover change and groundwater pumping (for public and private use) on the long-term groundwater level changes with HydroSight model. And the unconfined aquifer of Gnangara region in Western Australia was used as a case study.
Based on the three established independent models (climateonly analysis, climate and land cover analysis, and climate and pumping analysis models), climate always plays the most important role and positively contributes to the observed groundwater level changes. And groundwater decline in this region is mainly caused by the reduction in rainfall recharge over time. In the whole region, the contribution of the groundwater pumping to groundwater level decline is larger than that of land cover change. Temporally, from 1992 to 2014, the contribution of rainfall on the groundwater level of the Gnangara region decreased because the rainfall decreased over time. The mean pumping contribution rate is −23%, and the impact of groundwater pumping on the groundwater level decline continually increased from 1992 to 2014 because of the sustained groundwater pumping. The land cover changes had influence on long-term groundwater level changes, but the influence is limited and is stable over time with contribution rate of 2%. Spatially, in the southern Gnangara region, the groundwater level changes were mainly influenced by rainfall and pumping activities. In the northern Gnangara region, the groundwater level changes were influenced by the combination of rainfall, land cover and groundwater pumping.
The results of this study suggest that the improved groundwater hydrograph decomposition method is effective and can be easily applied in other regions due to its highly flexible. And this method improved the efficiency of data utilization, especially for the region which the groundwater head record is irregular. And the best-fit model for a certain study area can be obtained by trying different model structures. The findings of this study have important implications for research on the influence of various drivers on groundwater level changes and also provide notable guidance for local governments to rationally allocate and utilize groundwater resources. In areas where the groundwater level is mostly affected by groundwater pumping, other water resources should be utilized, such as rainfall runoff collected during the wet season for park irrigation, seawater desalinized, and surface water quality improved, while the groundwater abstraction reduction could ease the stress resulting from groundwater level decline.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.