Responses of groundwater to precipitation variability and ENSO in the Vietnamese Mekong Delta

Groundwater is a critical component of water resources and has become the primary water supply for agricultural and domestic uses in the Vietnamese Mekong Delta (VMD). Widespread groundwater level declines have occurred in the VMD over recent decades, reflecting that extraction rates exceed aquifer recharge in the region. However, the impacts of climate variability on groundwater system dynamics in the VMD remain poorly understood. Here, we explore recent changes in groundwater levels in shallow and deep aquifers from observed wells in the VMD and investigate their relations to the annual precipitation variability and El Niño–Southern Oscillation (ENSO). We show that groundwater level responds to changes in annual precipitation at time scales of approximately 1 year. Moreover, shallow (deep) groundwater in the VMD appears to correlate with the ENSO over intra-annual (inter-annual) time scales. Our findings reveal a critical linkage between groundwater level changes and climate variability, suggesting the need to develop an understanding of the impacts of climate variability across time scales on water resources in the VMD.


INTRODUCTION
Groundwater is the largest freshwater resource on Earth and a critical supply for human systems (Anderson et al. 2005). It has become the primary source of water used for domestic, industrial, and agricultural sectors in many countries (Giordano 2009). Globally, groundwater supports more than 40% of irrigation and is the sole water source for crops in most arid regions (Shiklomanov 2000;Siebert et al. 2010). Groundwater also plays a pivotal role in sustaining ecosystem services (Griebler & Avramov 2015) and surface water systems through flows into lakes and baseflow to rivers in the dry seasons (Tharme 2003).
In key food-producing regions, groundwater is economically important for a variety of reasons. First, it offers reliable and flexible access to a cost-effective source of freshwater with minimal infrastructure requirements that irrigation canals can hardly match. Second, groundwater is generally less prone to environmental pollution and degradation than surface water. Finally, its general continuity of water supply during dry seasons provides an essential buffer against episodic droughts. Despite its critical importance, groundwater often attracts insufficient management and much less attention in comparison to more visible surface water supplies in the rivers and reservoirs (Famiglietti 2014). precipitation variability and ENSO across time scales. The rest of this paper is organized as follows. Section 2 describes the study area, data used, and the analysis methodology. Section 3 discusses the main results of this study followed by conclusions in Section 4.

Study area
The VMD is located in the southern part of Vietnam, bounded by the Gulf of Thailand to the southwest, the East Sea to the south and southeast, and Cambodia to the north (see Figure 1(a)). It covers an area of about 40,000 km 2 and comprises a complex, dense network of canals (À50,000 km long in total or 80 m/ha density) and dykes. This network has been expanded continuously over time due to agriculture and human settlement (MRC 2015). Rice is the dominant crop in the VMD cultivated on À1.9 million ha, representing about 50% of the agricultural production land (Sebesvari et al. 2012). The surface of the VMD is relatively flat, with elevation ranging from 0.3 to 2.0 m above the mean sea level (MSL). The stratigraphy of the VMD is highly heterogeneous, including seven distinct aquifers separated by various aquitards composed of predominately marine silt, silty clay, and clays (SVGMD 2004;Wagner et al. 2012; see Table 1 and Figure 1(b)). The average thickness of these aquifers varies from several to almost a hundred meters. The lithology of each aquifer consists of fine to coarse sand, gravel, and pebbles (Bui et al. 2015).
The climate in the VMD is influenced by two monsoon systems: the ISM and WNPM, causing two distinct dry (from December to April) and wet (from May to November) seasons. The warmest temperatures and highest moisture occur during the wet season when the ISM brings humid air from the Indian Ocean. The annual rainfall varies from 1,300 to 2,500 mm, and about 80-90% of the rainfall amount occurs during the wet season. Similarly, river flow in the VMD shows pronounced seasonality. Discharge in Tan Chau station varies by from 5,000 to 17,000 m 3 /s from dry to wet seasons (Figure 2(a)). Approximately 75-85% of the total annual river flow occurs in the wet season, putting a high pressure of water demand on the groundwater system during the dry season. Groundwater in the VMD is exploited mainly for domestic (42%), agricultural (40%), and industrial (18%) uses (Bui et al. 2015).  Nguyen et al. (2004) and Pham et al. (2019). Please refer to the online version of this paper to see this figure in color: https://doi.org/10.2166/nh.2021.024.

Data
Groundwater level and precipitation observations over the VMD are used for analyses in this study. Groundwater level observations in Ho Chi Minh (HCM) city, the largest city of Vietnam, are also included. The locations of the monitoring wells are shown in Figure 2(b). Specifically, monthly groundwater level data observed at different aquifers over the period 2000-2018 for 130 wells throughout the region are obtained from the Vietnam National Center for Water Resources Planning and Investigation (NAWAPI). Groundwater wells are relatively well distributed across the VMD and classified into two depth categories: shallow and deep. Shallow wells are defined as those in the Holocene (qh) and Upper Pleistocene (qp 3 ) aquifers, whereas deep wells are those in the remaining aquifers: Upper-Middle Pleistocene (qp 23 ), Lower Pleistocene (qp 1 ), Middle Pliocene (n 22 ), Lower Pliocene (n 21 ), and Middle Miocene (n 13 ). Monthly observations of precipitation over the same period at 10 meteorological stations in the region were obtained from the Vietnam Meteorological and Hydrological Administration (VMHA). Monthly local precipitation time series of each well are obtained from the nearest meteorological station using a Voronoi diagram based on the network of 10 meteorological stations in the VMD (see Figure 2(b)). Numbers in parentheses show the percentage to the total number of wells in an aquifer. We use the Oceanic Niño Index (ONI), defined as the 3-month running means of sea surface temperature (SST) anomalies in the region 5°N-5°S and 120°W-170°W, as an indicator of the ENSO. It is a measure of the departure from normal SST in the central-to-eastern equatorial Pacific Ocean, smoothed out across multiple months to reduce noise and spikes. The monthly ONI during the study period is obtained from the National Oceanic and Atmospheric Administration (NOAA) at https://psl.noaa.gov/data/climateindices/list/. The El Niño (warm) phase is defined when the ONI is above þ0.5°C, whereas the La Niña (cold) phase is defined when the ONI index is below À0.5°C. A neutral phase is defined when the ONI is between À0.5°C and þ0.5°C (see Figure 3).

Groundwater trend analysis
Trend analysis is applied to investigate recent changes in the groundwater level in the VMD. For all wells, monthly time series of groundwater levels are first used to perform the nonparametric Mann- Kendall test (Mann 1945;Kendall 1955) to determine the existence of monotonic trends in the variation of groundwater level. Due to the brevity of observations, we consider groundwater records with a 10% significance level (p , 0:1) in rejection of the null hypothesis of no trend to be significant. Second, the linear slope of the time trend of groundwater level is calculated using the Theil-Sen method (Sen 1968). This trend estimator is typically used in conjunction with the Mann-Kendall test and is robust to nonlinearity and outliers. Due to the presence of missing or unusable data, trend analysis is applied only for wells that have at least 10 years of observation. Wells meeting this requirement numbered 116, including 50 in the shallow and 66 in the deep aquifers.
To estimate the trend of groundwater levels over the entire VMD, the Barnes analysis (Barnes 1964) is used for the spatial interpolation of groundwater level trends from the discrete wells in both depth categories. First, a weighted spatial average using the sum of Gaussian decay functions around each measurement point is generated. Then, the initial estimate is improved by adding the calculated error surface until a convergence factor is met. Barnes analysis with a region of interest set to À50 km is applied only to groundwater wells showing a statistically significant trend. Barnes analysis offers a weighted spatial analysis and includes error analysis to refine the result, which provides better interpolations than the inverse distance weighting and some other methods.

Precipitation control on groundwater
A vector autoregression (VAR) process for bivariate time-series data at the annual time scale is used to evaluate the causal relationship between local precipitation variability and groundwater level in the VMD. We assume here that annual groundwater level change is a linear function of past lags of itself and local precipitation anomaly. Specifically, a p-lag VAR model for groundwater level change based on the lagged groundwater level change and precipitation anomaly values is given as: where is the vector of present and lagged annual precipitation anomaly, is the vector of present and lagged annual groundwater level change, vector of white noise, and t represents time. The annual groundwater level change is calculated as the difference in groundwater elevation between two subsequent years. Precipitation anomaly is obtained by removing the monthly mean of each month during the study period. To isolate the signal of global warming, all groundwater level change and precipitation anomalies are also detrended prior to analysis by removing a linear trend from data using a least squares best-fit line. Equation (1) is solved using the ordinary least squares method (Ltkepohl 2007). The models for x 2,t are fitted using a lag of up to 4-year (p ¼ 4) using well records of at least 10 consecutive years. The x 1,t is not modeled but provides data so that the impulse response of X 2 given X 1 could be computed. A moving average (MA) representation of the VAR process is then used to calculate the impulse response of annual precipitation anomaly on groundwater level change, measuring the dynamic reaction of the groundwater level when precipitation is shocked by an impulse. The strength of the impulse response is measured as the coefficient in the MA representation of the VAR process (see details in Ltkepohl (2007)). Finally, for results with stable VAR coefficients, a Granger causality test (Granger 1969) at the significance level of 10% is used to determine the causal relationship between annual precipitation and groundwater level. The null hypothesis here is that the coefficients of the past precipitation values b k 21 in the regression Equation (1b) is zero. In other words, we quantify the contribution of the past lags of precipitation anomaly in explaining the variations of present groundwater level change.

Principal component analysis
Since it is difficult to directly analyze the high-dimensional groundwater system, a principal component analysis (PCA; see Abdi & Williams (2010) also known as empirical orthogonal function (EOF) analysis (Lorenz et al. 1956)) is used to reduce the groundwater datasets to a few modes that capture the dominant trends and patterns of the groundwater level dynamics. PCA includes three main components: eigenvectors (EOFs; representing the spatial patterns), principal components (PCs; representing the dominant variation of the time series in all wells), and eigenvalues (explaining the variance of the PCs). In this study, we extract the first two leading PCs to assess the relationship of groundwater level dynamics in each depth category and the ENSO.

Wavelet analysis
The wavelet coherence analysis (WCA) is used to explore the cross-correlation of each PC and the ONI time series across frequencies. The WCA combines continuous wavelet transform (CWT) with cross-spectrum analysis (Addison 2002). Our study employs the Morlet analytic wavelet that allows robust estimation of the power spectra and cross-spectra in the frequency domain even from relatively short time series while allowing for reasonably good frequency localization. Consider a time series x [ (x 1 , x 2 , ::, x n ) [ R n , the Morlet CWT of x is defined as: where C Ã (t) is the complex conjugate of the Morlet wavelet C(t) ¼ p À1=4 (e i2pf 0 t À e À(2pf 0 ) 2 =2 )e Àt 2 =2 , n is the scale parameter, and f 0 is the central frequency of the Morlet wavelet.
, often used in practice, such that the magnitude of the second highest peak of the wavelet is half the magnitude of the central peak (Addison 2002). In short, the CWT projects the time series x into the time-frequency domain by convolving the time series with a dilated and translated mother wavelet. For time series x and y with different wavelet transform w x (n, t) and w y (n, t), the wavelet coherence of these time series is given by: where Á h i represents a localized smoothing operation in both time and wavelet scale performed on the constituent transform components (Torrence & Compo 1998;Grinsted et al. 2004). The relative phase relationship between x and y is shown as Hydrology Research Vol 52 No 6, 1285 small arrows, with in-phase (positive) pointing right (0°), antiphase (negative) pointing left (180°), and x leads (lags) y by 90°p ointing downward (upward). The significance level for each scale is computed by the Monte Carlo method using only values outside of the cone-of-influence (Grinsted et al. 2004). Wavelet transform analysis is blind to linear trends and a more effective method than the Fourier transform to analyze multiple time scales of hydro-meteorological time series, which can truly reflect the periodicity of the given series and the distribution characteristics of the main cycles (Addison 2002).

Trends of groundwater level
The distribution of the slope of the time trend of groundwater level for all wells in the VMD during 2000-2018 is presented in Figure 4 and Table 1. The results show that most of the groundwater wells (101/116, equivalent 83.6%) exhibit negative slopes of change, indicating considerable decreases of groundwater levels over the entire delta. This dominant declining trend could primarily be attributed to the rapid increase of groundwater extraction in the region. Meanwhile, only 10/116 wells (8.6%), all in the shallow aquifers (7/26 wells in the qh and 3/24 wells in the qp 3 ), show a modest rising trend during the study period. In addition, 9/116 wells (7.8%) show no statistically significant time trend of groundwater level, mostly in the shallow aquifers. These results are consistent with previous studies (Wagner et al. 2012;Duy et al. 2021) and the fact that most of the major groundwater extraction wells in the VMD are located in deep aquifers where groundwater storage is more stable and less prone to pollution than groundwater in shallow aquifers. Consequently, given the slow groundwater recharge from the upper aquifers, groundwater levels tend to decline more rapidly in deep aquifers in the VMD. Figure 5 further shows the time trend of groundwater level interpolated from wells that exhibit statistically significant rates of change. For shallow aquifers, the magnitudes of the slopes are largest and relatively uniform in the central part of the delta with an average value of À0.1-0.2 m/yr ( Figure 5(a)), also overlapping with regions where rice crops are intensively cultivated (Nguyen et al. 2015). The decrease of groundwater levels in shallow aquifers may be correlated with the increasing groundwater use for irrigation during the dry seasons in the VMD, and this could be worth further exploration. We also observed no statistically significant time trend of groundwater level in the Southwest (Ca Mau province) and Northwest (along the border between Vietnam and Cambodia) in the shallow aquifers (see Figure 1(a) for these locations). For deep aquifers, groundwater level decline is much larger than that observed in the shallow aquifers over the delta. The magnitudes of the slopes appear to increase from the central of the delta toward the Southwest (Ca Mau province) and Northeast (Long An province and Ho Chi Minh city) regions with the largest rate of change is approximately 1 m/yr ( Figure 5(b)). In these regions, deep groundwater is used extensively for aquaculture and domestic uses. The smallest rates of groundwater decline in deep aquifers are observed in the Northwest (An Giang province) at about 0.2 m/yr. Overall, statistically significant declining trends of groundwater level are observed widely over the delta, especially in deep aquifers.
The temporal variation of total annual precipitation averaged over the VMD and the monthly groundwater level averaged over all groundwater wells for shallow and deep aquifers are presented in Figure 5(c) and 5(d). We note here that no significant trend of annual precipitation is detected in the VMD during the study period. Since aquifers can act as filters, it is obvious that groundwater levels in the shallow aquifers exhibit more pronounced seasonality signals than those in the deep aquifers. This filtering can also lead to long recharge travel time, depending on local geology, land use and land cover, and other factors affecting infiltration and recharge rates. Therefore, considerable decreases of the groundwater level suggest that the primary drivers of long-term groundwater level change can be attributed to the increase of deep groundwater abstraction in the VMD. A summary of groundwater level change is presented in Table 1.

Impacts of local precipitation anomaly on groundwater level change
Figure 6(a) shows the impulse response of local precipitation variability on annual groundwater level change as a function of time. Our results suggest that the groundwater response is largest in the first year of the precipitation change for both depth categories. The groundwater responses to precipitation is damped out to almost the fifth year. Indeed, the time required for groundwater recharge to travel from the land surface to the groundwater depends on soil moisture conditions and the properties of each aquifer, ranging from months to years (Russo & Lall 2017). Our results, therefore, may reflect the rapid transmission of groundwater recharge in the VMD region. The spatially interpolated map of the first VAR coefficients b 1 21 with stable VAR results, representing the response strength of groundwater level change to precipitation anomaly is further shown in Figure 6(b). We observe positive values of b 1 21 in most wells in the shallow aquifers, indicating a positive correlation between local groundwater responses to changes in precipitation anomaly at annual time scales. However, it is not the case for deep aquifers where negative values of b 1 21 are found in the south and some locations in the east, implying a negative correlation between groundwater level and precipitation anomaly in these areas. The Grainger causality tests indicate that only 21% (10%) of the individual wells in the shallow (deep) aquifer found annual precipitation is Granger-causal for predicting groundwater level at the significance level of 10%. This result suggests that annual precipitation causes changes in groundwater level in only a small area of the VMD. The remaining contribution is expected to come from precipitation associated with large-scale climate variability at longer time scales discussed below. Figure 7 shows the wavelet analysis of the first and second PC time series obtained from 50 wells in the shallow aquifers. The normalized PC1 time series exhibits an obvious declining trend of these wells over time (Figure 7(a), red solid line). In shallow aquifers, PC1 alone explains 62.6% of the groundwater level's variability. The corresponding squared magnitude of wavelet coefficient jWj 2 indicates a dominant periodicity of 1-year period of the PC1 (Figure 7(c)). The PC2 explains only 15.3% of the system's variability, but showing a higher correlation with the ONI or ENSO dynamics (Figure 7(b)). PC2 time series also exhibits dominant variations over longer (2-4 years) time scales in comparison to PC1 (see Figure 7(d)).

Groundwaterclimate connection
The wavelet coherences between ONI and the leading two PC time series are shown in Figure 7(e) and 7(f), respectively. The significant coherences between ONI and PC1 are found in seasonal and annual scales (2004, 2006-2007, 2010-2011, 2017) with phase differences mainly from 120°to 160° (Figure 7(f)), indicating that ONI leads PC1 from 3 to 6 months. We further observe that ONI also exhibits a high correlation with PC1 at time scales of 2-4 years from 2006 to 2016, with a phase difference ranging from 60°to 90°(indicating ONI leads PC1 between 9 and 12 months). This fact likely reflects the influences of ENSO on the increasing groundwater demand and pumping in the VMD during this period. For PC2, the significant coherences are found mostly at seasonal and annual scales (2002-2005, 2007-2008, 2011-2012, and 2014-2017) with phase differences from 0 (in-phase) to 90°(ONI leads 4-6 months).
Wavelet analysis of the first two PCs obtained from 66 wells in deep aquifers through PCA is presented in Figure 8. Similar to shallow aquifers, PC1 also exhibits a declining trend of groundwater levels in the deep aquifers of the VMD (Figure 8(a)). PC1 and PC2 explain 60.7% and 10.5% of the groundwater system variability, respectively. The squared magnitude of the wavelet coefficient jWj 2 indicates stronger periodicity of groundwater level in deep aquifers at annual and inter-annual scales compared with that in the shallow ones (Figure 8(c) and 8(d)). Both PCs show coherent patterns of jWj 2 around The wavelet coherence between ONI and the leading two PC time series of groundwater levels at deep aquifers is shown in Figure 8(e) and 8(f), respectively. At the annual time scale, the significant coherences between ONI and PC1 series are found only in 2008-2010 with a phase difference of about 80-90°, indicating that PC1 lags ONI by approximately 3-5 months. For PC2, significant and larger coherences are also found during 2008-2010 with a phase difference of about 200-220°. This phase difference suggests that the ONI lags the PC2 by 7-8 months. At the inter-annual time scale, significant coherences between ONI and PC2 are found during 2006-2014 with ONI leads PC2 by 12-18 months. Overall, groundwater levels in deep aquifers are highly correlated with ENSO as both PCs show large coherences at the inter-annual time scales, suggesting strong groundwater-climate connections in the deep aquifers over the VMD. Considering the expected attenuation of climate and recharge signals filtered through the upper aquifers, our findings suggest an indirect but strong influence of ENSO on groundwater demand via pumping that often has immediate local impacts on groundwater mass change.

CONCLUSION
In this study, we have investigated the changes in the groundwater levels in shallow and deep aquifers of the VMD and their relations with local precipitation variability and the ENSO. Our primary conclusions are as follows.
1. Statistically significant and notable declining trends of groundwater levels are found over the entire VMD in both shallow and deep aquifers. These trends are primarily attributed to the increase of groundwater extraction rates in the VMD, especially in the deep aquifers. 2. In addition to groundwater pumping, the Granger causality test indicated the influence of annual precipitation variability on groundwater levels, with the largest strength response of groundwater to local annual precipitation variability of approximately 1 year. However, the results showed the influence only in a small area fraction (10-20%) of the delta, which may also be obscured by groundwater pumping. 3. In shallow (deep) aquifers, the leading two PCs account for 78% (70%) of the total variance and significantly associated with the ENSO on the intra-annual (inter-annual) time scales. On average, the ENSO signal leads groundwater dynamics in shallow (deep) aquifers by 6-12 months (8-18 months).
The analysis herein suggested that climate variability, including local precipitation and the ENSO, plays an important role in the dynamics of groundwater in the VMD. Given projected changes in future temperature and precipitation patterns in which the VMD will be warmer with longer and drier summers, it is very likely that the groundwater system in the VMD in the future will face grand challenges to meet the increasing water demand. This also has implications for potential impacts of the changes in the ENSO-related precipitation variability under global warming on the groundwater system. Our results highlight a critical need to develop adaptive strategies for groundwater abstraction management in the VMD. Due to limited measurements of groundwater in the delta, it is not feasible to investigate the connection of groundwater change with major climate modes at longer time scales, i.e., the Pacific Decadal Oscillation (PDO). However, this work can be extended to study the connection of groundwater in the VMD with major climate modes of variability at shorter time scales such as the Madden Julian Oscillation (MJO). Our findings also encourage potential future research on understanding the contributions of monsoon, i.e., ISM and South Asian Summer Monsoon (SASM) that play influential roles on wet-season precipitation in the VMD (Irannezhad et al. 2020), on the groundwater dynamics.