Introduction

Aquifer systems play an important role in the global hydrological and biogeochemical cycles, as well as in the sustainability of ecosystems. (Zektser and Loaiciga 1993; Jackson et al. 2001; Sophocleous 2002; Griebler and Avramov 2015). Groundwater is a vital resource that supplies water to billions of humans, serving domestic, farming, and industrial activities (Giordano 2009; Siebert et al. 2010). With continued growth in population and economic activities, as well as changes in water consumption patterns, water demand is expected to increase significantly in the future. Meanwhile, changes in temperature and precipitation, driven by climate change, will further affect water availability, resulting in water shortages around the world, especially in developing countries. Such demand will likely affect not only surface water use but also groundwater use. If not properly managed, these water resources will be depleted quickly. The evolution of groundwater utilization has induced declined groundwater levels in numerous parts of the world (Konikow and Kendy 2005; Gleeson et al. 2012; Wada et al. 2012). Konikow (2011) estimated that there was a global depletion in groundwater of 41.4 km3/year from 1900 to 2008, posing a potential threat to water security and potentially leading to a global reduction in agricultural productivity and energy production, and also increasing the risk of conflicts around the world (Famiglietti 2014). Approximately 1.7 billion people live in locations where groundwater resources are under stress (Gleeson et al. 2012).

Groundwater resources, despite their significance, are frequently undermonitored, and precise assessment of groundwater-storage change has been challenging due to a lack of in-situ groundwater level data (Famiglietti et al. 2011). Inadequate geographical and temporal coverage by monitoring networks is a general problem, but stringent data-sharing rules (for political reasons) with other regions of the world may further hinder investigations on groundwater-storage change using in-situ well data (Chen et al. 2016). In terms of modelling techniques, the lack of in-situ groundwater level monitoring reduces the accuracy of simulated GWS change, especially over long timescales. Additionally, parameter values are frequently unknown or difficult to quantify across a large area, making reliable modelling of long-term GWS change challenging (Chen et al. 2016).

On March 17, 2002, the American National Aeronautics and Space Administration (NASA) and the German Aerospace Center (Deutsches Zentrum für Luft- und Raumfahrt-DLR) launched the Gravity Recovery and Climate Experiment (GRACE), the first committed satellite time-variable gravity mission that provides an unprecedented technique for tracking large-scale mass variations in the Earth system (Tapley et al. 2004). GRACE provides a monthly spatiotemporal terrestrial water storage anomaly (TWSA) that includes snow water, surface water, soil water, and groundwater storage (Voss et al. 2013). Because GRACE is unable to separate individual components of TWS change, changes in these components must be assessed using other data sources, such as land surface models (LSMs), in order to apply GRACE datasets to estimate changes in GWS (Rodell and Famiglietti 2002; Chen et al. 2016). Using advanced land surface models, the NASA Global Land Data Assimilation System (GLDAS) integrates satellite-based and in-situ observation data to offer optimal fields of land surface states and fluxes such as soil moisture, snow depth and other variables (Rodell et al. 2004). Many studies have successfully used GLDAS in combination with GRACE datasets to estimate GWS change over regions around the world (Rodell et al. 2009; Tiwari et al. 2009; Castellazzi et al. 2018; Singh et al. 2019; Gao et al. 2020; Wang et al. 2020). GRACE datasets have also been used in SE Asia to assess the variation of TWS in the Lancang-Mekong River Basin (Jing et al. 2020), to evaluate TWS and identify flood events in the Tonle Sap Basin (Tangdamrongsub et al. 2016), and to estimate the variation of subsurface water storage (soil moisture + groundwater) in the Lower Mekong River Basin (Pham-Duc et al. 2019).

Groundwater is a significant resource in Cambodia, and it has constantly been used for farming and domestic purposes alongside surface water (Phalla and Paradis 2011). Access to groundwater is crucial to increase the adoption of double cropping (Chan et al. 2004; Johnston et al. 2013). Furthermore, during the dry season, groundwater accounts for more than half of overall drinking water consumption (Sandqvist et al. 2008). Apart from domestic and agricultural uses, it is progressively being used in the industrial sector (ODC 2016b). Due to the extensive use of groundwater for agricultural and industrial purposes, its sustainability may be compromised by excessive withdrawal (Johnston et al. 2013). Between 1996 and 2008, an average groundwater level reduction of 14 cm/year was reported in the country’s southeastern regions (Prey Veng and Svay Rieng provinces), indicating that over-pumping may already be a concern in this region (IDE 2009). Erban and Gorelick (2016) discovered that the groundwater-irrigated area in the Cambodian Delta is expanding at a rate of more than 10% per year, and that if this rate continues, the groundwater level will fall below the lift limit of suction pump wells used for household purposes by more than 1.5 million people throughout much of the area within 15 years. Thus, assessing and understanding the GWS changes, stresses and resilience of the groundwater system is crucial for effective groundwater resource planning and management in the country. However, information about groundwater in the country is currently scarce; two major research projects have been completed in the last 20 years, but there has been no national plan or cooperative research (ODC 2016b). According to the Mekong River Commission (MRC) Strategic Plan 2016–2020 report, MRC had plans to execute a range of activities to establish practical knowledge on surface water and groundwater capacity and to assess the potential for agricultural water use in the Lower Mekong River Basin. Activity 1.7.2: “Conduct a survey of current groundwater use and the potential of new developments” and the Activity 1.7.3: “Conduct a study on groundwater sustainable yield management for crop production” are two major groundwater activities. However, due to the COVID-19 pandemic, the groundwater component of the irrigation study was not finished and, up to now, the study has only shown evidence of raising awareness and has promoted shared understanding among the member countries (MRC 2022). Meanwhile, remote sensing-based GRACE datasets can complement existing monitoring networks and modelling studies, as well as helping to compensate for gaps in spatial and temporal GWS-change coverage. However, these datasets have not yet been used to quantify GWS change, aquifer stress, and groundwater system resiliency in Cambodia. Due to a scarcity of surface-water datasets, the effect of surface-water storage (SWS) change signals from GRACE-based TWS signals is rarely addressed. The present study also investigates the impact of changes in SWS in Tonle Sap Lake, South-East Asia’s largest freshwater lake, on deriving changes in GWS in Cambodia, which will provide a good understanding of lake-water dynamics over GWS change in Cambodia.

Based on the foregoing discussion, the present study aims to assess the GWS variation, stresses and resilience of the groundwater system in Cambodia using GRACE datasets. The specific objectives of this study include (1) assessing GWS change from GRACE and GLDAS datasets, (2) assessing the effect of Tonle Sap Lake water-storage dynamics on the GWS change in Cambodia, (3) evaluating the performance of GRACE-derived estimates with respect to observation well data, and (4) evaluating groundwater/aquifer stress (GAS) and aquifer resilience (AR) in Cambodia.

Study area and data

Site description

Cambodia is located in South-East Asia, and its climatic condition is tropical and dominated by monsoon air masses. The months of May–October are associated with the southwest monsoon, with the highest precipitation occurring between September and October, and the remaining months are associated with the northeast monsoon, which comes with dry and cooler air. The average annual rainfall in the country ranges between 1,000 and 1,500 mm, and the heaviest amounts fall in the southeast (GlobalSecurity 2019). The topographical condition can be divided into three main terrains: mountains and plateaux of the border areas (mountains in the west, north, north-east, and south, as well as the Eastern Plateau); central plains that occupy almost three-quarters of the land area; and the southwestern coastal plain, as shown in Fig. 1a. Cambodia has an area of 181,035 km2 with a population about 16 million (WorldAltas 2020) and it consists of six river basins based on HydroBASINS Level 5 (Lehner and Grill 2013) with the Tonle Sap Lake basin being the largest at 83,657 km2; the areas of other major basins are 26,095 km2 (basin 1), 26,059 km2 (basin 2), 6,698 km2 (basin 3), 15,238 km2 (basin 4), and 19,449 km2 (basin 5) (Fig. 1a). Tonle Sap Lake and Mekong River are major surface water bodies in the country (Fig. 1a). The Tonle Sap Lake is the largest freshwater lake in South-East Asia, and its volume varies significantly annually in response to rainfall and flooding from the Mekong River. The permanent size of the lake is approximately 2,400 km2 (Kummu and Sarkkula 2008). The catchment area of the Tonle Sap Lake, including the permanent lake area, is the area of the lake–floodplain system. The volume of the system varies from 1.8 to 58.3 km3 in the driest month and 31.1 to 73.9 km3 in the peak water level month (Kummu and Sarkkula 2008).

Fig. 1
figure 1

a Digital elevation model (DEM) map showing six river basins and significant water bodies in the Cambodia regions, including 24 observation wells and 1 gauging station in Tonle Sap Lake; b Boundaries of the four subbasins covering the observation wells

Full size image

Data

GRACE data

The unique GRACE mission has two identical satellites orbiting around the Earth at about 500 km altitude with an inclination of 89.5°, apart from each other by approximately 137 miles (Tapley et al. 2004). Three different processing centers continuously release different solutions using various methods and parameters, including the Center for Space Research (CSR) at the University of Texas in Austin (USA), GeoforschungsZentrum Potsdam (GFZ) in Germany, and the Jet Propulsion Laboratory (JPL) in Pasadena, California (USA) (Cooley and Landerer, 2019). GRACE datasets are comprised of three different datasets. The level 1 dataset is calibrated and time-tagged raw data, the level 2 gravity field data has a set of spherical harmonics (SH) coefficients which refer to the exterior potential gravity field of the Earth system, and the level 3 products are changes in Earth surface mass which are presented by an equivalent unit of water thickness in cm. For detailed data processing of each level of data, please refer to Case et al. (2010), Bettadpur (2018), and Cooley and Landerer (2019). The most recent versions of GRACE TWSA products include spherical SH versions produced using a standard SH approach (Landerer and Swenson 2012), and mascon data versions processed using mass concentration blocks (Scanlon et al. 2016). However, according to Cooley and Landerer (2019), the current gridded mascon data are encouraged to be used instead of SH data because they do not require destriping or smoothing and suffer fewer leakage errors; furthermore, using the coastline resolution improvement (CRI) filter and gain factors, this approach allows for better separation of land and ocean signals; additionally, calculating basin averages for hydrology applications shows general agreement between both solutions for large basins. However, mascon data have higher spatial resolution for smaller regions; therefore, the monthly GRACE TWSA level 3 product by JPL for land release 06 version 2 using mascon solutions (Wiese et al. 2018) with a spatial resolution of 0.5° × 0.5° from April 2002 to March 2017 (180 months) were used in this study, and these anomaly values are relative to a 2004–2009 time-mean baseline. The data were downloaded from the NASA GRACE Data Analysis Tool. The gain factor was retrieved from PODAAC (Physical Oceanography Distributed Active Archive Center at JPL) Drive. Due to the aging batteries on the GRACE satellites, active battery management began in 2011. During specific orbit periods over several consecutive weeks, no range data were collected, and hence no gravity fields could be computed. These gaps occur approximately every 5–6 months and persist for about 4–5 weeks, resulting in approximately 12% of the data (22 months) missing over the study period, and these missing data were filled by averaging the two consecutive months preceding and following the missing month in this study.

GLDAS data

Global Land Data Assimilation System datasets are produced by NASA, in collaboration with National Centers for Environmental Prediction (NCEP), and National Oceanic and Atmospheric Administration (NOAA). GLDAS used an advanced land surface model and data assimilation technique to generate optimal fields of land surface states and fluxes by integrating satellite-based data and observed data (Rodell et al. 2004). GLDAS drives multiple offline (not coupled to the atmosphere) land surface models, integrates a massive amount of observation-based data, and executes globally at coarser resolutions (2.5° or 250 to 1 km) enabled by the Land Information System (LIS), producing near real-time results (typically within 48 h of the present; Rodell et al. 2004; Kumar et al. 2006). Previous research has shown that the NOAH model correlates better with the GRACE TWSA than Community Land Model (CLM), Mosaic, and Variable Infiltration Capacity (VIC) model (Leblanc et al. 2009; Grippa et al. 2011; Yang et al. 2013), and it is also recommended by Jet Propulsion Laboratory as well (JPL 2018). Therefore, this study used monthly GLDAS products processed using the last version 2.1 of the NOAH model (Rodell et al. 2004) with a spatial resolution of 0.25° × 0.25° for soil moisture storage (SMS) in four soil layers (0–10, 10–40, 40–100, and 100–200 cm) with a total depth of 200 cm of soil. These datasets were obtained from the NASA Giovanni online achieve. Climatic information from the same GLDAS products, including precipitation (P) and evapotranspiration (ET), were also used to evaluate GAS and AR in the study area from 2003 to 2016. These datasets were obtained from the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC).

Observed groundwater level data and aquifer formations

Due to the limitations of observed groundwater level data throughout the study area, only observed groundwater level elevation (h) in 24 monitoring wells located in some parts of major basins 4 and 5 (as shown in Fig. 1a), at a monthly frequency from 2006 to 2008, were used to validate the remote-sensing-based estimates from the GRACE and GLDAS datasets. These data were originally collected for a groundwater study in Cambodia by the Seila Task Force Secretariat (STFS), the National Committee for Sub-National Democratic Development, in collaboration with the Ministry of Water Resources and Meteorology (MoWRAM; IDE 2009). The observed groundwater levels ranged from –2.2 to 12 m relative to mean sea level (msl), with a mean of 3.5 m msl, as presented in Fig. 2. However, there were some missing values in the groundwater level data from observation wells for 3 months in 2006, during April, October, and November, and these missing values were filled using interpolation methods as GRACE data. Information on the aquifer material where the monitoring wells are located was also taken from the previous groundwater study report. Nine of the total wells are located in the old alluvium aquifer and the remaining wells are in young alluvium. The old alluvium is composed primarily of coarse-grained sand and gravel and has a thickness of up to 200 m, whereas the young alluvium is composed primarily of fine-grained clay and silt and has a thickness of 10–40 m (IDE 2009).

Fig. 2
figure 2

Time series monthly in-situ groundwater level between 2006 and 2008

Full size image

Observed surface water level in the Tonle Sap Lake

A recorded daily time series of water levels from April 2002 to March 2017 at Kampong Luong station in Tonle Sap Lake was collected from Mekong River Commission (MRC) online portal, for estimating the surface-water storage in Tonle Sap Lake; the location of the gauging station is shown in Fig. 1a. The daily time series lake water level for the entire study period is depicted in Fig. 3. The average lake water level was 4.5 m msl, with the maximum and minimum values of 12.1 and 0.8 m msl, respectively. Similar to missing data in groundwater level datasets, lake water level data were absent for some periods and were filled using the same interpolation procedure as GRACE data. Following gap filling, these water levels were substituted into an area-volume function to determine the extent and volume of the lake water. Table 1 summarizes all the required datasets.

Fig. 3
figure 3

Daily time series Tonle Sap Lake water level at Kampong Luong station between April 2002 and March 2017

Full size image
Table 1 Data for analysis, the related specifications, spatial and temporal resolution, time period and source
Full size table

Methodology

Deriving GWSA from GRACE TWSA datasets

Because the GRACE TWSA comprises surface water, soil moisture, snow, and groundwater storage components, obtaining the groundwater storage anomaly GWSA is as simple as subtracting other components from the TWSA, as indicated in Eq. (1) (Frappart and Ramillien 2018):

$${\mathrm{GWSA}}={\mathrm{TWSA}}-\left({\mathrm{SWSA}}+{\mathrm{SMSA}}+{\mathrm{SNSA}}\right)$$
(1)

where SMS is total soil moisture in the storage depth of the soil, SWS represents all surface-water storage, including lakes, reservoirs, weirs, and channel storage, SNC is total snow water storage, and letter A represents their anomaly values.

In this study, all the collected datasets were preprocessed. Preprocessing generally involves filling gaps in missing data, transforming data to a standard/needed datum projection, and resampling to a required resolution. GRACE TWSA in each grid was then multiplied by gain factors to reduce leakage error and resampled from the original resolution of 0.5° × 0.5° to 0.25° × 0.25° in order to obtain a consistent resolution with GLDAS data; additionally, the nearest neighbour interpolation method (Johnsy 2021) was used for resampling TWSA data in order to preserve their original values.

GLDAS-based SMS values for each soil layer were summed together to get the total soil moisture storage. To be consistent with TWSA, SMSA values were also estimated using the same baseline period of 2004–2009. Monthly SMSA values were estimated after subtracting baseline average values from monthly time series SMS values.

The relationship between Tonle Sap Lake water level (SWL), inundated area (A, in km2), and water volume (V, in km3) was used to convert the water level into the equivalent depth of water over the lake in centimeters. The open water area and volume of the lake water vary as a direct function of SWL and are expressed as follows (Kummu et al. 2014):

$$A=5.5701\times {{\mathrm{SWL}}}_{{\mathrm{KL}}}^3+1.374\times {{\mathrm{SWL}}}_{{\mathrm{KL}}}^2+470.29\times {{\mathrm{SWL}}}_{{\mathrm{KL}}}+1680.2$$
(2)
$$V=0.7307\times {{\mathrm{SWL}}}_{{\mathrm{KL}}}^2-0.3544\times {{\mathrm{SWL}}}_{{\mathrm{KL}}}+0.9127\kern0.5em$$
(3)

where SWLKL is water level depth [m] at Kampong Luong station.

A series of steps were carried out to estimate SWSA from lake water levels data as follows:

  1. 1.

    Observed daily water level from April 2002 to March 2017 was averaged to get monthly water level.

  2. 2.

    The estimated monthly water levels from step 1 were substituted into Eqs. (2) and (3) to obtain monthly A and V, respectively.

  3. 3.

    An equivalent depth of water was calculated from V by dividing by the average monthly area (Aavg). Aavg was calculated from estimated monthly A values using a simple arithmetic mean function.

  4. 4.

    Computed monthly time series equivalent depths of water from step 3 were then converted into SWSA by removing baseline time mean values.

Because the spatial extent of surface water was unknown, the Aavg area was divided with a 0.25° × 0.25° grid (i.e., ca. 27.75 km × 27.75 km) to estimate the number of grids covered over the surface water. The locations of surface water were assigned to the grids that are nearest to the permanent area of Tonle Sap Lake. Then the estimated monthly SWSA value for individual months was assigned to all the grids covering the Tonle Sap Lake area.

Since there is an absence of snow in the study area, the snow cover storage component was not considered in the present study. Canopy water storage fluctuations have been shown to be well below the detection limit of GRACE (Rodell et al. 2005; Liesch and Ohmer 2016). Therefore, both snow cover storage and canopy cover storage components were not considered in the present study; thus, GWSA values were simply obtained by removing SWSA and SMSA values from TWSA values as shown in Eq. (4):

$${\mathrm{GWSA}}={\mathrm{TWSA}}-\left({\mathrm{SWSA}}+{\mathrm{SMSA}}\right)$$
(4)

Trend estimation from GWSA signals

The nonparametric Mann-Kendall and Theil-Sen’s slope tests were adopted to estimate the trends and slope of the trend of the time series signals, respectively. Trend analysis of the time-series GWSA values was separately performed at pixel, basin, and country scales. Each of the time series GWSA values was used to estimate the trend and slope of the trend for each pixel in the study area. Consequently, spatially averaged GWSA values resulting from the basin and country boundaries were used to estimate the trend and slope of the trend at basin and country scales, respectively.

Mann-Kendall (MK) test

The MK test was applied to detect significant changes in GWSA time series datasets. This test is based on the test statistic Stest, which is defined as follows (Chatterjee et al. 2016):

$${S}_{{\mathrm{test}}}=\sum_{i=1}^{n-1}\sum_{j=i+1}^n\operatorname{sign}\left({x}_j-{x}_i\right)$$
(5)
$$\operatorname{sign}\left({x}_j-{x}_i\right)=\left\{\begin{array}{c}1\ {\mathrm{if}}\ {x}_j-{x}_i>0\\ {}0\ {\mathrm{if}}\ {x}_j-{x}_i=0\\ {}-1\ {\mathrm{if}}\ {x}_j-{x}_i<0\end{array}\right.$$
(6)

where x1, x2,…, xn are n data points and xj is the data point at time j. Where the Stest value is positive, an increasing trend is indicated, and where the Stest value is negative, a decreasing trend is reflected. When n is equal to or higher than 10, the statistic Stest has an approximately normally distributed mean (E) and variance (var), as demonstrated in Eqs. (7) and (8), respectively:

$$E\left({S}_{{\mathrm{test}}}\right)=0$$
(7)
$$\operatorname{var}\left({S}_{{\mathrm{test}}}\right)=\frac{n\left(n-1\right)\left(2n+5\right)-\sum_{i=1}^m{t}_i\left({t}_i-1\right)\left(2{t}_i+5\right)}{18}$$
(8)

where m represents the number of tied ground (a set of sample data having the same value) values and ti represents the number of data points in the ith group. The standardized test statistic Z is calculated as follows:

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{\operatorname{var}(S)}}\ {\mathrm{if}}\ S>0\\ {}0\ {\mathrm{if}}\ S=0\\ {}\frac{S+1}{\sqrt{\operatorname{var}(S)}}\ {\mathrm{if}}\ S<0\end{array}\kern3em \right.$$
(9)

The null hypothesis (H0), which means that an insignificant trend exists, is accepted if Z is not statistically significant, i.e. −Zα/2 < Z < Zα/2 where α is the selected significant level; a p-value of 0.05 was chosen as a threshold of significance in this study.

Theil-Sen’s estimator

The slope of change in groundwater storage was calculated using Theil-Sen’s estimator. According to Lavagnini et al. (2011), this approach is much less sensitive to outliers than simple linear regression approaches since it calculates the median slopes of lines fit through pairs of points in the dataset. This method has been applied to quantify magnitudes of the trends in hydrology as well as climate data records, and it is often used in conjunction with the MK test (Li et al. 2014). The Theil–Sen’s estimator by Sen (1968) and Theil (1992), which is used to estimate the slope of n pairs of data points, is defined by the following relation (Eq. 10):

$${T}_i=\frac{x_j-{x}_i}{j-i}$$
(10)

x j and xi are data values at time j and i (j > i) respectively. Sen’s estimate of slope is the median of the N values of Ti, which is derived as follows:

$${Q}_i=\left\{\begin{array}{c}{T}_{\frac{N+1}{2}}\ {\mathrm{if}}\ N={\mathrm{odd}}\\ {}\frac{1}{2}\left({T}_{\frac{N}{2}}+{T}_{\frac{N+2}{2}}\right){\mathrm{if}}\ N={\mathrm{even}}\end{array}\right.$$
(11)

The sign of Qi denotes trend direction, while its value represents the trend’s steepness.

Validating GRACE-based GWSA with respect to observation well data

The comparison between remote-sensing-estimated GWSA and in-situ GWSA was performed at subbasin scale. There are 10 subbasins that cover the monitoring well locations, and they were dissolved into only four subbasins for validation purposes, as shown in Fig. 1b. Another set of GRACE TWSA and GLDAS SMS data covering validation subbasins from 2006–2008 were prepared for validation purpose. TWSA and SMS values were converted to anomaly values with respect to this 3-year mean time value so that resulting GRACE-derived GWSA values were able to consistently be compared with in-situ well data.

The monthly in-situ GWSA was estimated using observed groundwater level depth, h, in the monitoring well and storativity (S) of the aquifer (Todd and Mays 2005; Bhanja et al. 2017, 2018) as follows:

$${{\mathrm{GWSA}}}_{{\mathrm{obs}}}={h}_{{m}}\times S-{h}_i\times S$$
(12)

where GWSAobs is in-situ-based GWSA, hm and hi are the mean groundwater level and groundwater level that are recorded from an observation well at any time i, respectively, and S is storage coefficient/storativity, the ability of an aquifer to release a volume of water from the storage per unit of decline of hydraulic head per unit of the area of the aquifer. As the aquifer formation where the wells are located was known, the S value used in this study was approximated based on the Johnson (1967) information as presented in Table 2.

Table 2 Storage coefficient (S) for different materials (Johnson 1967)
Full size table

Because the young alluvium was constituted primarily of fine-grained clay and silt, its S value was determined by averaging the S values of clay and silt. Similarly, the S value of the old alluvium was determined by averaging the S values of coarse sand and gravel, as it was formed primarily from coarse-grained sand and gravel (IDE 2009). The procedure to calculate in-situ GWSA values is as follows:

  1. 1.

    At first, monthly groundwater level was multiplied by the storage coefficient for each well to get monthly GWSobs.

  2. 2.

    Monthly GWSobs values for each well were converted to monthly GWSA using 2006–2008 baseline mean values. Another set of GRACE TWSA and GLDAS SMS data was prepared, covering validation subbasins from 2006 to 2008, and then these were converted to anomaly values with respect to this 3-year mean time value so that the resulting GRACE-derived GWSA values for validation with in-situ GWSA are consistently compared with each other.

  3. 3.

    The obtained in-situ GWSA values from the previous phase were interpolated over subbasin boundaries using the inverse distance weighting method to construct spatially continuous raster grids. This approach estimates the unknown point based on the concept that the farthest known point with a record gets less weight than the closer one, and its formulation is given by Eq. (13) (GISGeography 2020):

$${{\mathrm{GWSA}}}_{{\mathrm{int}}}=\frac{\sum_{i=1}^n\frac{{{\mathrm{GWSA}}}_i}{d_i^2}}{\sum_{i=1}^n\frac{1}{d_i^2}}$$
(13)

where GWSAint is an unknown point value, GWSAi is a known point value, d represents the distance between the known and the unknown point, and n is the total number of known points.

The performance of remote-sensing estimates with respect to in-situ groundwater level measurements was assessed using statistical measures such as root mean square error (RMSE), and correlation coefficient (r). The RMSE and r values between monthly GRACE-based and in-situ GWSA for each subbasin were computed using the following formulations (AgriMetSoft 2019):

$${\mathrm{RMSE}}=\sqrt{\frac{\sum_{i=1}^n\left({x}_i-{y}_i\right)}{n}}$$
(14)
$$r=\frac{n\sum_{i=1}^n\left({x}_i{y}_i\right)-\sum_{i=1}^n\left({x}_i\right)\sum_{i=1}^n\left({y}_i\right)}{\sqrt{n\Big[\sum_{i=1}^n{x}_i^2-{\left(\sum_{i=1}^n{x}_i\right)}^2}\left]\times \sqrt{n\Big[\sum_{i=1}^n{y}_i^2-{\left(\sum_{i=1}^n{y}_i\right)}^2}\right]}$$
(15)

where xi represents the observed data, yi is the estimated data and n is the total number of observed data.

Estimating monthly GWS change

Monthly times series groundwater storage change (GWSC) values from 2003 to 2016 were computed from estimated monthly time series GWSA values using second-order central difference, as given by Eq. (16), since it is more numerically stable than the forward/backward difference (Biancamaria et al. 2019).

$${\mathrm{GWSC}}(m)=\frac{{\mathrm{GWSA}}\left(m+1\right)-{\mathrm{GWSA}}(m)}{2}$$
(16)

where GWSC(m) is the value of GWSC at month m; GWSA(m + 1) and GWSA(m – 1) are successive and preceding ‘m’ months values of GWSA.

Evaluating GAS and AR

In this study, the GAS values result from the subtraction of groundwater released from the renewable groundwater resource as defined by the following expression:

$${\mathrm{GAS}}={\mathrm{Renewable}}\ {\mathrm{groundwater}}-{\mathrm{Groundwater}}\ {\mathrm{release}}$$
(17)

Here, renewable groundwater is represented by potential groundwater recharge, which was estimated by subtracting monthly time series P from ET values considering only positive values, i.e., P – ET > 0.

Only estimated GWSC values with a negative sign were used to represent groundwater release (Eq. 17), which, in this study, refers to the amount of groundwater released from the aquifer by all mechanisms (e.g., extraction or natural groundwater flow). Therefore, GAS ≥ 0 means that no stress occurs because renewable groundwater is greater or equal to groundwater release, while GAS < 0 means that stress happens because release exceeds the renewable groundwater of the aquifer.

Resilience (AR) measures how quickly a system returns to a satisfactory state (i.e., recovered from shocks) following an unsatisfactory state (Hashimoto et al. 1982). In this study, groundwater resilience was determined according to Loucks (1997) and Mays (2013) as follows:

$${\mathrm{AR}}=\frac{{\mathrm{No}}.\kern0.5em {\mathrm{of}}\ {\mathrm{a}}\ {\mathrm{s}}{\mathrm{a}}{\mathrm{tisfactory}}\ {\mathrm{condition}}{\mathrm{s}}\ {\mathrm{following}}\ {\mathrm{an}}\ {\mathrm{unsatisfactory}}\ {\mathrm{condition}}\ \left({\mathrm{GAS}}\right)\ }{{\mathrm{Total}}\ {\mathrm{no}}.\kern0.5em {\mathrm{of}}\ {\mathrm{unsatisfactory}}\ {\mathrm{conditions}}\ \left({\mathrm{GAS}}\right)}$$
(18)

If the GAS value shifts from negative to positive, this is considered a recovery event. This study classified AR into five levels—0–0.25: very low resilience; 0.25–0.4: low resilience; 0.4–0.6: moderate resilience; 0.6–0.75: high resilience; and 0.75–1: very high resilience.

Results and discussion

Trends in GWSA values deriving from GRACE TWSA

The spatially averaged values and linear fitted trends in each of the component anomaly values for Cambodia for each month from April 2002 to March 2017 are depicted in Fig. 4. Monthly TWSA values were between –300 and +650 mm, and SMSA values ranged between ±200 mm. For the surface water component, the Aavg value of the lake water was 4,978 km2, which is equal to about the area of seven raster grid cells. The monthly SWSA values were then assigned to only these seven grids closest to the permanent Tonle Sap Lake location and, hence, these assigned grids are only located in basin 1. This assumption of surface-water locations, therefore, generates spatial uncertainty for grid-scale results; nevertheless, this assumption may not affect basin and national-scale results because estimated monthly surface water area fluctuated ≤17,388 km2, which is around ≤20% of the catchment area for basin 1. As a result, the surface-water component was considered for the basin 1 and country scales but not for the grid-scale analysis. The SWSA values at the national scale were calculated by averaging the SWSA values in each grid throughout the entire country’s boundary. The spatial average SWSA values ranged from –100 to +400 mm.

Fig. 4
figure 4

Spatial average values of anomalies of the terrestrial water storage (TWSA), total soil moisture in the storage depth of the soil (SMSA), surface water storage (SWSA) and groundwater storage (GWSA) for the whole study area, with their fitted linear trend, the slope of change in GWSA, and the p-value from the MK-test and Sen’s slope estimation

Full size image

The spatially averaged GWSA values varied between –250 and +300 mm, and it was noticed that SWSA and SMSA fluctuated in a seasonal cycle, whereas GWSA fluctuated in a complicated pattern, indicating that it is not just climatic factors that have contributed to this variation in groundwater storage. Furthermore, a decreasing linear trend in TWSA values was observed with a change of slope of 0.03 mm/month, with changes in groundwater storage accounting for approximately 70% (linear change slope = –0.0234 mm/month) of the variation in TWSA values, while changes in soil moisture and surface water storage accounted for only 10% (slope = –0.0037 mm/month) and 20% (slope = –0.0068 mm/month), respectively (Fig. 4). The MK-test and Sen’s slope estimation results for the trend and slope of change in GWSA are also shown in Fig. 4. A significant decreasing trend (p-value 0.05) in GWSA values was observed, with a slope of change of 0.63 mm/month, indicating that GWS declined significantly at a rate of approximately –113 million m3/month over the last recent 15 years, implying that groundwater storage has been overextracted in the study area.

The spatial average TWSA, SWSA, SMSA, and GWSA values, as well as the slope of change and linear trend in GWSA within each river basin, are depicted in Fig. S1 in the electronic supplementary material (ESM). Except in basin 2, GWS was found to be significantly decreased in basins 1, 3, 4, 5, and 6. The calculated slopes ranged from –0.14 to –1.51 mm/month. The most significant decline was found in basins 3 and 4, with declining rates up to –1.51 mm/month, almost 5 times higher than the slope values in basins 1 and 6. Similarly, GWS in basin 5 has reduced at a rate of roughly –1.05 mm/month. In basin 2, on the other hand, the rate of decreasing trend was assessed to be only 0.1 mm/month. These results indicate that groundwater extraction has occurred in some basins, such as basins 3 and 4, where the highest negative slopes were estimated, while a balance between recharge and discharge has been maintained in basin 2, where the slope of change is insignificant. It was noticed that GWSA values within most of the basins fluctuated according to the seasonal climate pattern (increasing in the wet season and decreasing in the dry season), but only within basin 1 did GWSA values reach their positive peak in a dry month (e.g., February 2005 and March 2006) and negative peak in a wet month (e.g., October 2005 and 2006), which means that this basin may be recharged by means other than direct precipitation infiltration such as losing streams (e.g., losing surface water from Tole Sap Lake), preferential flows, or through fractures.

The grid level analysis for slope of change and the significance test over Cambodia is presented in Fig. S2 of the ESM. According to this figure, the majority of the western regions (latitudes 10–13.5°N and longitudes 102–104.5°E) had p-values greater than 0.05, indicating that there were no significant trends in GWSA signals in these locations, whereas the remaining regions had significant trends, indicating that changes in GWS in these regions are critical. Significant decreasing trends in GWSA were found in the eastern areas and north-west borders, with rates ranging from –0.5 to –2 mm/month, whereas significant increasing trends were observed in the north-east borders, which is a mountainous area, with rates of +1.5 mm/month as presented in Fig. S2a of the ESM. Because the grids covered by basin 2 showed both significant increasing and decreasing trends, the resulting slope of change when all pixels were averaged over the basin’s boundary was not significant in the basin level analysis.

Influences of the surface water in Tonle Sap Lake on deriving GWS signals

The spatially averaged monthly time series GWSA values with their slope of change resulting from taking and not taking surface water in the lake, taking into consideration the basin 1 and country scales, are illustrated in Fig. 5a and b, respectively. When surface water is not included in deriving GWSA from TWSA, the slope of change in GWSA values decreases to about –0.73 mm/month for the country scale, which means that when surface water is not included in deriving GWSA from TWSA, there is a 15% overestimation of change in GWSA values for the country scale. However, when the basin scale was assessed, there are very substantial changes in both patterns and slopes of change between these two cases. The slope of change of GWSA values without surface water decreases to –0.69 mm/month, while the slope of change with surface water is only –0.43 mm/month. Therefore, if surface water is not taken into account, the decreasing rate of GWS change in basin 1 will be overestimated by 60%. In general, it has been shown that for the SWS component, even in the largest river basins in the world (such as the Mississippi River basin, USA), SWS changes have been shown to be at least an order of magnitude less than those of groundwater and soil moisture changes (Rodell and Famiglietti 2001; Rodell et al. 2007). However, this assumption is not valid in extreme flooding years or in basins that are extremely wet (i.e., rainforest regions). Consequently, SWS changes are often not considered in GRACE studies (Yeh et al. 2006; Rodell et al. 2007; Tiwari et al. 2009). In Cambodia, due to Tonle Sap Lake being the largest freshwater lake in Southeast Asia, the SWS component has been found to significantly impact the derivation of the GWS component from GRACE data. Therefore, the effect of large water bodies such as Tonle Sap Lake, must be taken into account in order to interpret the TWSA signals appropriately, as this could significantly change the statistical results and lead to misinterpretation of the effects on the groundwater storage signals.

Fig. 5
figure 5

Spatial average GWSA values with slopes of change estimated with and without considering surface water storage in the Tonle Sap Lake for: a basin 1 scale, b country scale. Red and blue lines represent linear fits of GWSA values with and without surface water, respectively

Full size image

Comparison and validation of GRACE-based and in-situ-based GWSA

The comparison of GRACE-based GWSA and in-situ-based GWSA values estimated from in-situ groundwater level data are illustrated in Fig. 6. This figure shows scatter plots of monthly GWSA derived from GRACE against in-situ field-based observation well GWSA values within 4 subbasins from 2006 to 2008. The scatter points for all subbasins were found to be reasonably scattered around the 1:1 line. The fitted linear models between these scatter points were also very well aligned with the 1:1 line, indicating that the remote-sensing-based GWSA and in-situ GWSA values reasonably agreed with each other. A good match of r = 0.58 was observed between the satellite-derived and ground-based measurements in the High Plains Aquifer in the USA, using 2,700 monitoring wells from 2003–2005 (Strassberg et al. 2007). Using 236 wells between 2003 and 2007, Shamsudduha et al. (2012) found reasonable matches between the two estimations in Bangladesh. Strong correlations (r > 0.70) between the two GWS estimations (GRACE and observed data from >15,000 groundwater observation wells between 2005 and 2013) have also been seen in 8 out of 12 Indian basins (Bhanja et al. 2016). The GRACE-derived groundwater storage change values agree reasonably with in-situ well observations (1,180 wells from May 2005 to November 2016) with correlations of up to 0.89 for four basins in India (Sarkar et al. 2020). Similar to the results here, estimated r values of 0.88, 0.82, 0.85, and 0.86 were found for the four subbasins of this study, confirming that remote-sensing-based GWSA estimates are strongly related to in-situ-based GWSA estimates. However, when it comes to the RMSE indicator, the estimated RMSE values range from 70 to 86 mm. This error could be caused by uncertainty in the approximation of the storage coefficient values used to convert water level to storage.

Fig. 6
figure 6

Scatter plots of GRACE monthly GWSA estimations vs observation well datasets for four subbasins in the study area. The blue line is the linear fit of the scatter point, whereas the grey line represents the 1:1 line

Full size image

Mean monthly GRACE-derived GWSC

The mean monthly spatial GWSC values estimated from GRACE-derived GWSA between 2003 and 2016 are presented in Fig. 7. The GWSC values were found to be in the range of –150 to 200 mm/month. The largest decreases in GWS were found in November and December, at the start of the dry season, with a maximum decreasing rate of –150 mm/month in the south-west mountainous region and the region surrounding the Tonle Sap Lake, which might have caused reductions in water storage/stream flow in the lake. Similarly, the highest decreasing rate was found in the remaining dry months until April at a mountainous area at the north-east boundary. Regarding the rainy season, the months of May and June, which are the beginning of the wet season, groundwater levels declined but not significantly. The greatest increase in groundwater storage was seen in July, August, and September, with a maximum increasing rate of >200 mm/month, and groundwater began to decrease at the end of the wet season, i.e., October.

Fig. 7
figure 7

Mean monthly GWSC values between 2003 and 2016, estimated from GRACE-derived GWSA values

Full size image

Despite the fact that overall decreasing trends in GWSA values in the south-west mountainous region and the Tonle Sap Lake region were not significant, it was discovered that extreme events, i.e., the greatest decrease and increase in mean monthly groundwater (Fig. 7), occurred at these locations, indicating that although these regions had the most abundant groundwater in the rainy season, they also have the highest risk of water scarcity or groundwater depletion. In addition, the effects on water storage/flow in the Tonle Sap Lake, causing threats to the lake’s ecosystem, should be taken into account. Meanwhile, severe droughts have hit Cambodia for many years, with the frequency of drought varying from province to province, with the most affected provinces being Kampong Speu, Takeo, and Battambang (ODC 2016a), all of which are located in these regions.

Potential groundwater recharge and stress on the aquifer

The spatial mean monthly PR values computed from the GLDAS datasets from 2003–2016 are shown in Fig. 8. Estimated PR values ranged roughly from 0 to 350 mm/month, with only up to 50 mm/month during dry months (November–April) but ranging from 50–350 mm/month during rainy periods (May–October), with a peak in September and significant spatial heterogeneity. Overall, the highest PR values (up to 350 mm/month in September) were detected in the southern and mountainous areas, indicating that these are the parts of the aquifer that are the most viable renewable source in the research area.

Fig. 8
figure 8

Spatial mean monthly potential recharge (PR) values between 2003 and 2016, estimated from GLDAS products

Full size image

After subtracting estimated groundwater release values (i.e., negative GWSC values) from estimated renewable groundwater values (i.e., potential groundwater recharge), the resulting spatial mean monthly GAS values for each month between 2003 and 2016 are shown in Fig. 9. It was discovered that the aquifer of the research area experienced groundwater stress during the dry period, with GAS values of up to –200 mm/month, with the most severe stress occurring at the beginning of the dry season through November to December. The strongest groundwater stress was encountered in the southwestern mountainous region, the area surrounding Tonle Sap Lake, and the mountainous area at the northeastern boundary, which had the highest mean monthly decreasing rate in GWS (Fig. 7). During the rainy season, positive GAS values of up to 300 mm/month were found, suggesting that the available renewable groundwater was more than the extracted amount from the aquifer, and hence no groundwater stress emerged. The highest positive GAS values (GAS ≥ 300 mm/month) were observed in September in the regions receiving the highest potential recharge (Fig. 8). According to Oeurng (2020), during the wet season, potential groundwater recharge is more than sufficient in most of the country. Likewise, the results of this study reveal that the research area has had an abundance of renewable groundwater during the rainy season, while the aquifer experienced groundwater stress during the dry season, and this stress may be even stronger in reality because the actual amount of renewable groundwater could be less than the potential renewable groundwater values used in this study.

Fig. 9
figure 9

Spatial mean monthly groundwater aquifer stress (GAS) values between 2003 and 2016

Full size image

Figure 10 illustrates the monthly time series GAS values for each river basin from 2003 to 2016. During the dry season (November–April), all basins faced groundwater stress, with basin 6 experiencing the most stress (GAS ≥ –200 mm/month) and basin 2 experiencing the least stress (GAS ≥ –100 mm/month). The highest negative GAS values were detected in December 2011 in all basins (excluding basin 2), indicating that the most serious groundwater stress occurred during this period, while Cambodia also experienced its worst seasonal flooding in a decade during the same year (Tritt 2011). Therefore, 2011 was the worst year in which the country was hit by two of the worst water-related hazards consecutively; however, lower groundwater stress was found in the last 4 years (from 2013 to 2016). During the rainy season (May–October), GAS values were positive in the majority of the basins (basins 2–6), indicating that no groundwater stress occurred in the aquifer; however, even during the wet season, negative GAS levels were detected in basin 1 (e.g., August 2009 and 2011). Nevertheless, because this study used potential recharge from climate data, and this basin might be replenished by additional means other than direction infiltration of rainfall (as discussed in section ‘Trend estimation from GWSA signals’), and thus the temporal pattern of groundwater release was different from the renewable source, the resulting aquifer stress in this basin should be interpreted carefully.

Fig. 10
figure 10

Spatially averaged GAS values within each river basin between 2003 and 2006

Full size image

Resiliency of the aquifer

The results of the quantification of the aquifer resilience are illustrated in Fig. 11. Approximately 17% and 4% of the area possessed low resilience (0.25 < AR ≤ 0.4) and moderate resilience (0.4 < AR ≤ 0.6), respectively, which were mostly seen along the country’s border, while the remaining roughly 78% of the whole region had very low resilience (AR ≤ 0.25). When basin-scale analysis was performed, the AR values ranged only between 0.21 and 0.29. Three out of all the river basins (basins 1–3) possessed a very low recovery ability, and the remaining basins exhibited low resiliency. Even though having the highest groundwater stress, basin 6 was the most resilient, while basin 2 had the weakest recovery ability, despite getting the lowest stress (Fig. 10). Overall, the Cambodian aquifer exhibited a very low ability to swiftly recover from stress during the 14 years of study. This clearly indicates the impact of recent developments on the groundwater system in the region. Therefore, continuous monitoring, planning and regulation of the groundwater system against present and planned developments in the region is necessary for the sustainable management of groundwater systems in Cambodia (Table 3).

Fig. 11
figure 11

Spatial groundwater aquifer resilience (AR) values between 2003 and 2016

Full size image
Table 3 Resilience of the aquifer for each river basin
Full size table

Conclusions

For better monitoring and management of the groundwater dynamics in Cambodia, remote-sensing-based datasets were used to evaluate GWS change, aquifer stress and resilience in the country. The effects of surface water in Tole Sap Lake on deriving GWS change from GRACE TWSA datasets likewise were evaluated with the comprehensive validation of GRACE-based GWSA with well observation data as were the aquifer stress and its resiliency to the stress. Important conclusions from this study are drawn as follows:

  1. 1.

    The correlation coefficient values between the remote-sensing and in-situ GWSA values are up to 0.88, implying that the remote-sensing estimates and ground-based observations have good consistency. This correlation analysis shows the potential of GRACE and GLDAS datasets for monitoring and assessing groundwater storage dynamics throughout Cambodia.

  2. 2.

    At the national scale, groundwater storage has been decreasing at the rate of –0.63 mm/month, i.e., ca. –113 million m3/month, over the period April 2002–March 2017, with basins 3 and 4 experiencing the highest declining rate of more than 1.4 mm/month, and thus effective measures must be implemented to avoid depletion of the groundwater resources in them.

  3. 3.

    Understanding the impact of Tonle Sap Lake water storage on groundwater storage change is critical to understanding the surface-water and groundwater interactions and dynamics over the years. The statistical results significantly change when the surface-water component from Tonle Sap Lake is not considered in the GWSA analysis; thus, neglecting the storage change of large water bodies in the GWSA analysis will propagate errors in the analysis, particularly at the basin scale.

  4. 4.

    During the dry season, the north-east border and south-west mountainous regions, as well as the region surrounding Tonle Sap Lake, had the greatest decrease in groundwater storage. Therefore, if no proper management is implemented, the greatest risks of water shortage will be during the dry season or groundwater storage depletion in these locations, as well as a reduction in water storage or flow in the Tole Sap Lake, which could lead to the destruction of the lake’s ecosystem.

  5. 5.

    During the rainy season, Cambodia has had a surplus of potential groundwater recharge; however, the aquifer has faced groundwater stress during the dry season, with basin 6 suffering the most stress and basin 2 experiencing the least aquifer stress over the 14 years studied.

  6. 6.

    The Cambodian aquifer had a very low ability to quickly recover from stress, with basin 6 being the most resilient and basin 2 being the weakest, emphasizing that appropriate management is urgently needed to ensure the sustainability of the country’s groundwater resources.

Despite the fact that this study validated these remote-sensing-based estimates using 24 observation well records, denser and more extensive networks of observation well information with a relatively longer period of records and storage coefficient information are expected to reduce the uncertainty in converting groundwater level to storage volume and to validate GRACE and GLDAS-based estimates further. In addition, more GRACE datasets from different versions and processing centers should be analysed in order to select the one with the best performance for the study area. Furthermore, information regarding the spatial extent of the lake’s surface water should be identified (through, for example, satellite image data) and added in future studies to show more clearly the spatial effects of the surface-water component on deriving GWSA signals from GRACE TWSA data. Finally, given the study’s findings of decreased groundwater storage, the effects of these change drivers, such as climate change, land use/land cover change, and groundwater pumping, should be investigated for better understanding and for more effective management of groundwater systems in terms of climatic, land use/land cover, and human intervention factors.