# The recharge illusion: How seasonal surface loading can hide continuing groundwater resource declines in confined deltaic aquifers

Felix Dörr, Jonas Bauer, Gabriel C Rau, Elias Lewi, Viet Tran Hoan, Le Thi

Mai Van, Remi Valois, Anke Steinel, Franz Nestmann, Stefan Norra

► To cite this version:

Felix Dörr, Jonas Bauer, Gabriel C Rau, Elias Lewi, Viet Tran Hoan, et al.. The recharge illusion: How seasonal surface loading can hide continuing groundwater resource declines in confined deltaic aquifers. Journal of Hydrology, 2026, 665, pp.134674. 10.1016/j.jhydrol.2025.134674. hal-05444583

HAL Id: hal-05444583

https://hal.inrae.fr/hal-05444583v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# The recharge illusion: How seasonal surface loading can hide continuing groundwater resource declines in confined deltaic aquifers

Felix Dörr<sup>a,*</sup>, Jonas Bauer<sup>a,†</sup>, Gabriel C. Rau<sup>b,‡</sup>, Elias Lewi<sup>c,§</sup>, Viet Tran Hoan<sup>a,d</sup>, Le Thi Mai Van<sup>d</sup>, Remi Valois<sup>e,¶</sup>, Anke Steinel<sup>f,§</sup>, Franz Nestmann<sup>g</sup>, Stefan Norra<sup>h</sup>

## ARTICLE INFO

This manuscript was handled by Renato Morbidelli, Editor-in-Chief, with the assistance of Eddie Banks, Associate Editor

### Keywords:

- Loading efficiency

- Mekong delta

- Barometric efficiency

- Passive subsurface characterization

- Groundwater head dynamics

## ABSTRACT

Reliable estimates of groundwater recharge are critical for managing water resources and preventing groundwater depletion. However, groundwater recharge estimates can be distorted when loading signals are mistakenly interpreted as signs of groundwater recharge. This study introduces a novel extension to existing frequency-domain methods to quantify loading imprints in groundwater measurements by disentangling superposing signals from Earth, atmospheric, ocean, and river tides, using the $K_1$ tidal constituent to remove the river tide influence. This approach enables the application of frequency-domain methods for calculating the loading efficiency in deltaic groundwater systems, where superposing multiple forcings govern tidal constituents in groundwater heads. Using multi-factor regression deconvolution analysis, surface loading signals are removed from groundwater head time series. Applying this method to groundwater head data from the southern Vietnamese Mekong Delta strongly suggests that the observed seasonal rise in groundwater heads does not reflect recharge. Instead, the apparent upward trend in uncorrected data is attributed to mechanical surface loading from seasonal surface water accumulation linked to expanding aquaculture. The loading-corrected time series show a consistent decline during both dry and rainy seasons, aligning with observed annual depletion rates. This indicates a deficiency of net effective groundwater recharge throughout the year, reflecting unsustainable groundwater exploitation practices. The methods and findings of this study offer a scientific framework for the analysis of multiple superimposed surface loading processes in deltaic aquifer systems, thereby enhancing groundwater recharge assessments in such hydrogeological environments.

## 1. Introduction

Global groundwater resources are under increasing stress, and the widespread acceleration of depletion trends underscores the urgent need for effective measures to address ongoing declines (Jasechko et al., 2024). Coastal areas are particularly vulnerable to the impact of groundwater overexploitation, as groundwater depletion promotes seawater intrusion and thus further accelerates the loss of available underground freshwater resources (Konikow and Kendy, 2005).

Moreover, in many deltaic lowlands, the depletion of groundwater heads causes notable land subsidence due to the compaction of subsurface sediments (Herrera-García et al., 2021; Davydzenka et al., 2024). The Asian mega-deltas like the Ganges-Brahmaputra Delta, the Yangtze Delta and the Mekong Delta are prominent examples for such subsiding deltas (Svitski et al., 2009). The southern Vietnamese Mekong Delta (VMD) is subject to land subsidence dynamics with rates locally exceeding $0.03 \text{ m a}^{-1}$ (Erban et al., 2014; Minderhoud et al., 2017; Dörr et al., 2024). Given its average elevation of only about 80 cm above mean sea level in 2018 (Minderhoud et al., 2019), such

<footnote><sup>a</sup> Karlsruhe Institute of Technology (KIT), Institute of Applied Geosciences, Kaiserstraße 12, 76131 Karlsruhe, Germany</footnote>

<footnote><sup>b</sup> School of Environmental and Life Sciences, The University of Newcastle, University Drive, NSW 2308 Callaghan, Australia</footnote>

<footnote><sup>c</sup> Institute of Geophysics, Space Science and Astronomy, Addis Ababa University, Arat Kilo Campus, Addis Ababa, Ethiopia</footnote>

<footnote><sup>d</sup> National Center for Water Resources Planning and Investigation (NAWAPI), No. 93, lane 95, Vu Xuan Thieu Street, Sai Dong Ward, Long Bien District, Hanoi 100000, Viet Nam</footnote>

<footnote><sup>e</sup> Avignon University, INRAE, UMR 1114 EMMAH, F-84000 Avignon, France</footnote>

<footnote><sup>f</sup> Federal Institute for Geosciences and Natural Resources (BGR), Stilleweg 2, D-30655 Hannover, Germany</footnote>

<footnote><sup>g</sup> Karlsruhe Institute of Technology (KIT), Institute of Water and River Basin Management, Kaiserstraße 12, 76131 Karlsruhe, Germany</footnote>

<footnote><sup>h</sup> Potsdam University, Institute of Environmental Sciences and Geography, Soil Sciences and Geocology, Campus Golm, Building 12, 14476 Potsdam Golm, Germany</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: felix.doerr@kit.edu (F. Dörr).</footnote>
## Nomenclature

### Acronyms and Abbreviations

BE Barometric Efficiency

BRF Barometric Response Function

ERF Earth Tide Response Function

EORF Earth- and Ocean Tide Response Function

FFT Fast Fourier Transform

GRACE Gravity Recovery and Climate Experiment

HALS Harmonic Least Squares

IRC Instantaneous Response Coefficient

LULC Land-Use Land-Cover

MUFACO Multi-Factor Correction of Groundwater Levels

OTGW Ocean Tide Loading on the Aquifer's Water

OTSE Ocean Tide Loading on the Solid Earth

PLRF Pond Loading Response Function

PT Pressure Transducer

RRF River Response Function

VMD Vietnamese Mekong Delta

pronounced subsidence represents an existential threat to the VMD.

While advancing salinization and increasing pollution often restrain the use of surface water resources in the southern VMD (Bauer et al. 2025; Binh et al. 2025), groundwater remains to serve as a crucial water resource for the region's water dependent socio-economic activities (Hoan et al., 2022, Pham et al., 2023). The overexploitation of groundwater causes an ongoing depletion of the hydraulic heads in the multilayered aquifer-aquitard system with rates of 0.01 to 0.55 m a⁻¹ across its seven aquifers (Duy et al., 2021) and sustainable water resource management remains an unsolved challenge.

To prevent further groundwater depletion and associated land subsidence, the determination of sustainable groundwater extraction rates is essential and therefore accurate and reliable groundwater recharge assessments become ever more relevant for the southern VMD and other subsiding deltas.

Despite long-lasting and multidisciplinary research efforts (Anderson, 1978; Pham et al., 2019; Duy et al., 2021; Bauer et al., 2022), groundwater recharge processes in the VMD remain uncertain at the present time. In some coastal provinces of the VMD, the topmost Holocene aquifer (qh) is unconfined and hence generally capable of receiving groundwater recharge (Wagner et al., 2012). However, in most parts of the southern VMD it is overlain by the Holocene aquitard (Q2). With thicknesses of up to > 20 m (Anderson, 1978; Pham et al., 2019, Fig. 1b) this aquitard constrains significant groundwater recharge by percolation of rain- or surface water resources where present. For the deeper Pliocene, Pleistocene and Miocene aquifers, the intermediate aquitards restrict substantial vertical flow from the qh aquifer to deeper aquifers for which no recharge areas have been identified in the VMD (Pechstein et al., 2018). A potential hydraulic interconnection between the deep aquifers in the southern VMD and outcropping Cambodian alluvial strata remains uncertain and a pending transboundary challenge of essential importance to the VMD's water security (Lee et al., 2018). This is underpinned by recent hydrogeochemical and isotopic analyses across the multilayered aquifer system, suggesting that most groundwater was recharged into the Upper Pleistocene aquifer $qp_{2.3}$ during the last glacial episode when seawater levels were low (Pham et al. 2019, Bauer et al. 2022, Osenbrück et al., 2025). However, these studies on the other hand indicate that local vertical flow or leakage processes through the aquitards may be present (Bauer et al. 2022, Osenbrück et al., 2025).

Hydrodynamic studies (e.g. Duy et al., 2021) investigated seasonal fluctuations in groundwater levels and their correlation with river water levels, suggesting that the deeper aquifers might receive recharge from upstream areas of the VMD. More recently, Van et al. (2023) applied the

water-table fluctuation method (Healy and Cook, 2002) in the southern VMD, a method originally derived for unconfined conditions, to estimate groundwater recharge based on observed seasonally rising groundwater heads and precipitation data.

However, in confined aquifers, groundwater head fluctuations are not solely caused by water mass fluxes, i.e. extraction and recharge, but also by mechanical load acting on the aquifer's confining unit, like for example variations in soil moisture or water levels of surface water bodies (van der Kamp and Maathuis, 1991). Such loading signals bear the potential of being misinterpreted as groundwater recharge, given that they cause an increase of water levels in observation wells (van der Kamp and Schmidt, 2017).

Over the past decade, hydrogeological studies in other subsiding delta systems, such as the Bengal Basin (e.g. Burgess et al., 2017), have identified and investigated loading effects in groundwater observation data. These studies highlight both the importance of accounting for such effects and the limitations of traditional analyses that ignore them, particularly in Southeast Asian mega-deltas. Dörr et al. (2025) uncovered the presence of loading signals in groundwater head time series in the southern VMD, paving the way to assess the contribution of loading signals to apparent seasonally rising groundwater levels in the southern VMD.

The objective of the present study is to identify and quantify the impact of mechanical loading processes due to the accumulation of surface water on seasonally rising groundwater heads in the southern VMD. To achieve this, this study introduces a novel extension to previous frequency-domain methods (Acworth et al., 2016; Rau et al., 2020), which utilize the $S_2$ frequency (2.00 cycles per day) for the calculation of the loading efficiency $\gamma$ (–), using $K_1$ signals for the removal of river tide influence on groundwater head amplitudes. The novel approach expands the applicability of such methods to deltaic groundwater systems, where superposing multiple tidal forcings govern tidal constituents in groundwater heads. While previous methods were only applicable in the absence of river and ocean tide loading (Acworth et al., 2016; Rau et al., 2020), the present method extends the recent inclusion of ocean tide loading (Dörr et al., 2025) by integrating river tide effects and thus, for the first time, enables a quantitative disentanglement of Earth-, atmospheric-, ocean-, and river tide components in groundwater heads.

Evaluating the influence of loading signals on groundwater heads by using the novel extension to existing methods developed in this study enhances the understanding of groundwater response to surface water loading in the study area and supports more accurate groundwater recharge assessments. The insights gained from the case of the VMD are transferable to other delta systems in Southeast Asia and beyond.

## 2. Materials and methods

### 2.1. Hydrometeorological data

To assess the potential contribution of loading effects on rising groundwater heads in the southern VMD, this study utilizes hydrometeorological data from the monitoring site GW1 (Fig. 1) in Ca Mau Province in the southern VMD, located in an area dominated by inland aquaculture (Fig. 1a; Vu et al., 2022; Van Binh et al., 2025) and an extensive multi-channel system (Le et al., 2023; Fig. 1c). The monitoring site GW1 was established by the ViWaT Engineering research project by the Karlsruhe Institute of Technology (KIT), the Vietnamese National Center for Water Resources Planning and Investigation (NAWAPI) and the Department of Natural Resources and Environment Ca Mau (DONRE Ca Mau).

At GW1, the top Holocene aquitard $Q_2$ has a thickness of approx. 32.5 m (Fig. 1b). Like at other study sites in the southern VMD (Pechstein et al., 2018), the development of the aquifers qh and $qp_3$ is not clearly distinguishable at GW1, indicated by the dashed blue sections in the aquifer characterization in Fig. 1b. The introduced color code for the aquifers will be used throughout this study.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Study area and location of the CM3 river station and the GW1 hydrometeorological monitoring site in the southern VMD as well as land use map from 2020 (Vu et al., 2022) (a), aquifer-aquitard categorization and screen depths of observation wells at GW1 (b), rivers and channels in the southern VMD (c) and hydro-geological cross section through the multi-layered aquifer-aquitard system of the southern VMD (modified after Hoan et al., 2022) (d). Note that the depth and height of soil moisture and weather sensors in (b) are not scaled.</figcaption></figure>

The following data from GW1 is utilized for the assessment of loading signals in groundwater heads: (i) hourly to 3-hourly groundwater heads obtained using vented pressure transducers (PTs) in three observation wells (GW1a, GW1b, GW1c) screened in the upper-middle Pleistocene aquifer (qp2-3), lower Pleistocene aquifer (qp1) and middle Pliocene aquifer (n2), respectively (Fig. 2); (ii) 15-minutely observations of the soil's volumetric water content (–) in nine depths between 0.05 and 1.00 m (Fig. 3); (iii) minutely observations from a meteorological station comprising data on air temperature, precipitation, relative humidity, solar radiation, barometric pressure and wind speed (Fig. 4).

The observed hydraulic heads in the confined aquifers qp2-3, qp1 and

n2 (Fig. 2a) show exhibit dynamics, with rising heads between end of May and end of November 2023 (rainy season, see Fig. 4d) and declining heads between end of November 2023 and early May 2024 (dry season, see Fig. 4d).

Fig. 2bI-III illustrates the different resolutions of the PTs of approx. ± 10–15 mm for GW1a and GW1b and approx. ± 2–5 mm for GW1c. The temporal resolution varies between 1 h (GW1b and GW1c, Fig. 2b-II and b-III) and 3 h (GW1a, Fig. 2b-I). The precipitation sensor has a measurement range of 0–400 mm/h and an accuracy of ± 5 % of measurement in the range between 0 to 50 mm/h. Fig. 3 illustrates the soil water content observations at GW1, showing notable variation between the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Hydraulic head observation data of three observation wells at GW1 screened in the upper-middle Pleistocene aquifer (qp<sub>2-3</sub>), lower Pleistocene aquifer (qp<sub>1</sub>) and middle Pliocene aquifer (n<sub>2</sub><sup>2</sup>) from January 2022 to July 2024 (a) and a 20-day zoom-in visualization between October and November 2023 (b).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Observation of soil volumetric water content at GW1 in different depths between 5 cm and 100 cm.</figcaption></figure>

dry and rainy season in depths up to 30 cm (Fig. 3a–c), while below 60 cm a nearly constant volumetric water content is observed (Fig. 3g–i). The sensor was calibrated by a soil specific calibration using laboratory water content analyses of soil samples from each observation depth (Ledieu et al., 1986).

In 5 km distance to GW1 (Fig. 1a), surface water levels are measured at the CM3-Rach Cui station at hourly resolution (Fig. 5). The surface water levels are characterized by seasonal variations with rising water

levels between July and November as well as by variations on tidal scale. These seasonal and sub-diel water level dynamics are induced by ocean tides as well as seasonal variations in the upstream feed of freshwater and comply well with the overall reported river water levels dynamics in the southern VMD (Bauer et al., 2025).
<figure><img src="image_6.png" /><figcaption>Fig. 4. Measurement data of the meteorological monitoring station at GW1 (colored lines) and daily averages (black lines).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Hourly data on surface water levels at the gauging station CM3 Rach Cui.</figcaption></figure>

## 2.2. Water balance and land use land cover classification

For the analysis of seasonal surface loading $\Delta\sigma_L$ (mm), the water balance in the proximity of the study site GW1 is analyzed. In the presented study, the general surface water balance of accumulated precipitation $P$ (mm), accumulated evapotranspiration $E$ (mm) and accumulated runoff $R$ (mm) (Eq. (1)) is considered:

$$ \Delta\sigma_L = P - E - R \quad (1) $$

Given that the site exhibits 32.5 m of aquitard at the surface, no actual unconfined aquifer exists at the site (Fig. 1b). Hence, there is no additional loading from unconfined shallow groundwater level changes on the assessed deep aquifers.

The responding area of an observation well to surface loading can be approximated as a circular area with a radius of at least ten times the depth of the observation well (Smith et al., 2017). For the calculations in this study, radii of ten times the depth of GW1a-c are utilized for the estimation of the observation wells' responding area, complemented by an impact assessment of potential uncertainties on the calculation results. Given that actual evapotranspiration and runoff depend on land use and land cover (LULC), in this study four types of LULC-classes are differentiated within the estimated responding surface area: (i) ponds,

(ii) vegetated land (iii), sealed surfaces as well as (iv) rivers and channels. These four types of LULC-classes differ in their accumulation of precipitation, their evapotranspiration rates as well as their drainage and runoff behavior.

In this study, based on field observations, LULC-areas identified as ponds are considered as undrained water bodies, supposing that they are not continuously connected to the river and channel network and therefore accumulate precipitation and do not contribute to surface water runoff. In the traditional extensive shrimp farming system, the ponds' water management was on tidal basis, linking the ponds' water levels to water levels in rivers and channels by fortnightly opening and closing of the ponds' sluice gates (Johnston et al., 2000). On the other hand, water level assessments of intensive shrimp farming ponds in the study area showed increasing pond water levels within a production cycle from initially 0.8–1.0 m to 1.5–2.0 m due to precipitation as well as the addition of supplemental water (Anh et al., 2010). Sources of used supplemental water were not specified by Anh et al. (2010) and are likely to vary according to the availability of different water resources, e. g. harvested rainwater or surface water in sufficient quality. Despite legal restrictions, a survey from Tra Vinh province suggests that also groundwater is used as supplemental water in shrimp farming (Van Tuan et al., 2024). Given that intensive shrimp farming systems account
for the largest share of recent shrimp production in the VMD (Ho et al., 2025), in this study aquaculture ponds near GW1 are considered as intensive shrimp farming ponds, a presumption that matches field observation and is validated by an evaluation of pond sizes within the analysis.

For ponds $\Delta\sigma_L$ is equivalent to variations in the pond water level and is calculated by the difference between accumulated precipitation and evaporation. The evaporation rate for ponds $E_w$ (mm d⁻¹) is calculated with Eq. (2), a simplified form of the Penman (1948) equation derived by Valiantzas (2006) for the application in settings with measurement data from routine meteorological stations:

$$E_w = 0.051(1 - \alpha)R_s\sqrt{T + 9.5} - 2.4\left(\frac{R_s}{R_A}\right)^2 + 0.052(T + 20)\left(1 - \frac{RH}{100}\right)(a_U - 0.38 + 0.54u) \quad (2)$$

with $\alpha$ (–) as the albedo (for open water surfaces $\alpha = 0.08$ (Allen et al., 1998)), $R_s$ as the measured solar radiation (MJ m⁻² d⁻¹), $T$ as the measured air temperature (°C), $R_a$ as the calculated extraterrestrial radiation (MJ m⁻² d⁻¹) based on the FAO56 method (Allen et al., 1998), $RH$ (%) as the measured relative humidity, $a_U = 1$ according to the original Penman (1948) wind function, and $u$ as the measured wind speed (m s⁻¹).

Areas identified as sealed surfaces like urban settlements do not accumulate water and temporary precipitation loading is completely converted into surface water runoff or evaporation. Therefore, areas classified as sealed surfaces are considered to be fully drained in this study.

The seasonal accumulation of water resources on vegetated land is challenging to determine precisely in the study area, given that several factors such as infiltration, plant water accumulation, inundation, irrigation practices and drainage runoff contribute to this process. On the one hand, the high clay content of the topsoil most likely prevents significant infiltration. In the present study this is assessed by a quantitative evaluation of the soil moisture loading, i.e., the variation of soil water volume, derived from the depth differentiated soil moisture observation at GW1 (Fig. 3). On the other hand, plant- and canopy water accumulation as well as inundation of farmland such as rice paddies add to the overall surface load on vegetated land. Given that these contributions cannot be quantified with the data available, this study does not utilize a single value for $\Delta\sigma_L$ for vegetated land but considers a range of possible values. In this approach a share between 0 % and 50 % of the areas classified as vegetated land are considered to accumulate surface water resources, characterized by the overall water balance dynamics of precipitation and evaporation, equivalent to the loading dynamics of ponds.

For LULC-areas classified as rivers or channels, the loading is directly derived from river water level observations in the proximity of GW1 (Fig. 5). For the LULC-classification, a detailed visual assessment of Google Satellite Imagery, dated to February 2024, was performed, given that remote sensing approaches, for example as used by Vu et al. (2022), have insufficient accuracy for the required site-specific analysis.

## 2.3. Loading signals in the time domain

Water levels in deep groundwater observation wells in confined aquifers, such as GW1a, GW1b and GW1c respond to variations in the water balance at the surface (van der Kamp and Maathuis, 1991). When hydrometeorological processes like soil moisture variations, ponding of surface water, barometric pressure variations or snow accumulation cause a spatially extensive and uniform surface load, the coupled stress-flow equations can be simplified to a 1D problem and the general relationship between aquifer pore-pressure response and surface loading can be written for saturated, linearly elastic, porous media using van der Kamp and Gale (1983):

$$\frac{\partial u}{\partial t} = \gamma \frac{\partial \sigma_L}{\partial t} + D \frac{\partial^2 u}{\partial z^2} \quad (3)$$

where $u$ is the pore pressure, $t$ is time, $\gamma$ is the loading efficiency, $\sigma_L$ is the total stress (mechanical load) acting at the top, $D$ is the hydraulic diffusivity, and $z$ is elevation.

The loading efficiency $\gamma$ (–) describes the partitioning of external stress that is imposed on an aquifer between the sediment matrix and the pore water (Domenico and Schwartz, 1998). The two terms on the right head side of Eq. (3) represent the two governing key processes, namely the undrained pore-pressure response to changes in mechanical load ($\gamma \frac{\partial\sigma_L}{\partial t}$), and the vertical transient flow of groundwater due to the pore-pressure variations ($D \frac{\partial^2 u}{\partial z^2}$) (van der Kamp and Schmidt, 2017).

While variations of the total stress are instantaneously transmitted into the subsurface to depths of kilometers (Schulze et al., 2000), the transient flow is inherently subject to delay, which is characterized by the time factor $T$ (s) as in Eq. (4) (van der Kamp and Schmidt, 2017)

$$T = \frac{b^2}{D} \quad (4)$$

where $b$ (m) is the thickness of the confining aquitard and $D$ (m² s⁻¹) its hydraulic diffusivity. The time factor $T$ scales the time frame for the aquifer-aquitard interaction (van der Kamp and Schmidt, 2017). In an ideal confined aquifer, $D$ converges to zero and the transient flow term is negligible, whereas in leaky aquifer systems a more careful assessment of the contribution of transient flow is necessary. If the duration of a loading signal is much smaller than $T$, the groundwater response will not be significantly affected by flow through the overlying aquitard (van der Kamp and Schmidt, 2017).

In the present study, the loading efficiency $\gamma$ was utilized to remove the loading signals due to mechanical loading ($\Delta\sigma_L$) from hydraulic head measurements in three observation wells at the monitoring site GW1. For that, $\gamma$ was evaluated by three independent methods: (i) the $S_2$ response of groundwater levels to atmospheric tides $BE_{S_2}^{AT}$; (ii) the groundwater level response to selected distinct precipitation events; (iii) the barometric response function (BRF).

The BRF is calculated using a regression deconvolution approach, where changes in groundwater pressure at a given time $t$ are expressed as the convolution of changes in barometric pressure at times before $t$, in the form of instantaneous response coefficients (IRC). The cumulative sum of the IRCs is referred to as BRF (Rasmussen and Crawford, 1997). Toll and Rasmussen (2007) enhanced this method by additionally considering Earth tide effects in the regression as an independent driver for variations of well water levels, represented by the Earth tide response function (ERF). Methods like the multi-factor correction of groundwater levels (MUFACO) (Haehnel et al., 2024) generalize the forcings to consider more processes, e.g., ocean tide and precipitation loading, therefore allowing an individual and more flexible consideration of multiple influences.

In aquifers with near-ideally confined conditions, the BRF typically increases towards a plateau. In that case, the BRF's maximum represents the barometric efficiency $BE$ (–) (Toll and Rasmussen, 2007). This is inherently coupled to the loading efficiency $\gamma$ by Eq. (5) (van der Kamp and Gale, 1983):

$$\gamma = 1 - BE \quad (5)$$

In semi-confined aquifers, the BRF decreases with time, indicating significant contribution of vertical transient flow. Therefore, the shape of the BRF can be used to characterize the aquifer confinement status (Rasmussen and Crawford, 1997; Butler et al., 2011; Spane, 2002). Moreover, once the response function to a certain forcing or multiple forcings are derived, they can be used to remove the imprints of the respective effects in the observed groundwater heads (Rasmussen and Mote, 2007; Toll and Rasmussen, 2007).

In addition to the evaluation of the BRF, the groundwater head re-
sponges to selected distinct precipitation events are used to assess the loading efficiency in this study. Given that short and distinct precipitation events are abrupt increments in the surface loading (Bardsley and Campbell, 2000), slow processes like evaporation as well as surface and subsurface drainage processes can be neglected, Eq. (3) can be simplified and the loading efficiency $\gamma$ can be calculated with the ratio of accumulated precipitation and increase in observed groundwater heads. To minimize the influence of low sensor accuracy at GW1a and GW1b (Fig. 2 b-I and b-II) and to reduce the impact of tidal forcings (Earth-, ocean-, and atmospheric tides), linear trends were fitted to the groundwater head data before and after each precipitation event. The durations of these trends were selected based on the longest available precipitation-free periods. The resulting start and end values provide a robust basis for the quantitative assessment of $\gamma$.

## 2.4. Multi-factor regression deconvolution

In the present study, a multi-factor regression deconvolution is applied to correct groundwater heads and to derive the BRF. This approach accounts for tidal and seasonal drivers influencing groundwater-head variations at the study site. The considered drivers include the ubiquitous drivers of (i) Earth tides and (ii) atmospheric pressure variations, as well as site-specific drivers including (iii) river water levels, (iv) ocean tide loading on the solid Earth (OT$_{SE}$), and (v) calculated variations in pond water levels (see section 2.2). Considering the aquaculture-dominated land use and the dense network of rivers and channels in the study area (Figs. 1a and c), pond water loading and river water levels control the seasonal surface loading signals in groundwater heads. The impact of observed river tides at CM3 on groundwater heads at GW1 is assessed using a magnitude-squared coherence analysis based on measured time series of river water levels and groundwater heads at the respective sites.

Volumetric Earth tide strains are computed for the monitoring site GW1 with PyGTide (Rau et al., 2022), a Python wrapper for ETERNA PREDICT 3.4 (Wenzel 1996), under consideration of longitude, latitude and elevation of GW1. The volumetric strains due to ocean tide loading on the solid Earth (OT$_{SE}$) are calculated with the Fortran-based software package SPOTL (Agnew, 2012) at hourly resolution for GW1. OT$_{SE}$ causes strain waves, propagating from the ocean coast towards the inland by the geologic media. In another independent process (denoted as OT$_{GW}$), ocean tides exert a loading on the aquifer's water (Jacob, 1950; van der Kamp, 1972), causing a pressure wave that propagates from the ocean coast towards the inland by the liquid phase (pore water) only. Dörr et al. (2025) showed that in the southern VMD, groundwater tides due to OT$_{GW}$ reach insignificance (< 0.01 mm) in 6.4 km distance in the middle Pliocene aquifer $N_2^2$. Given that the study site GW1 is in approx. 28 km distance to the coast, OT$_{GW}$ is considered irrelevant for the tidal composition of groundwater at GW1.

Building on previous studies in the VMD, the hydraulic diffusivity $D$ of the overlying aquitards $N_2^2$, $Q_1^1$ and $Q_1^{2-3}$ (see Fig. 1b) can be estimated from mean values of the vertical hydraulic conductivity and the specific storage coefficient based on Minderhoud et al. (2017) to $D(Q_1^{2-3}) = 2.3 \cdot 10^{-3} \text{ m}^2 \text{ s}^{-1}$, $D(Q_1^1) = 5.1 \cdot 10^{-4} \text{ m}^2 \text{ s}^{-1}$, and $D(N_2^2) = 5.5 \cdot 10^{-4} \text{ m}^2 \text{ s}^{-1}$. With aquitard thicknesses of 16.7 m, 9.2 m and 17.8 m at GW1 (Fig. 1b), the time factor $T$ can be estimated with Eq. (4) to be approx. 1.4, 1.9 and 6.7 days, respectively. Based on these time factors, a maximum time lag of 7 days (168 h) was selected in the multi-factor regression deconvolution, to account for potential contributions of vertical flow to the loading-groundwater response.

Within the multi-factor regression deconvolution, uncertainties are assessed in form of the response functions' standard deviations, which are derived from the covariance matrix of the IRC during the regression deconvolution (Haehnel et al., 2024).

## 2.5. Loading signals in the frequency domain

Complementary to the evaluation of the BRF, this study also assesses BE in the frequency domain by an evaluation of groundwater level response to atmospheric tides at $S_2$ frequency (2.0 cycles per day). For that, the underlying frequency patterns in the time series of Earth-, ocean-, river- and atmospheric tides, as well as groundwater heads are analyzed by two independent methods: (i) the Fast Fourier Transform (FFT); (ii) the harmonic least-squares (HALS) estimation. In this study, HALS analyses are conducted on the ten most studied frequencies in groundwater analysis, as outlined in Table 1. Prior to any performed frequency analysis, the time series are detrended with a 3-day moving window to remove lower frequencies and improve the extraction of the relevant tidal constituents (Rau et al., 2020). Given that HALS outperforms FFT in terms of accuracy in amplitude and phase estimation (Schweizer et al., 2021), amplitudes and phases obtained from HALS are used for the assessment of BE. The additional consideration of FFT offers redundancy to the HALS-based amplitude and phase characterization and has the potential to detect amplitudes beyond those listed in Table 1.

Within the HALS estimation of amplitudes and phases, a covariance matrix for the fitted coefficients can be calculated and used to derive the standard deviation (Rau et al., 2020).

## 2.6. Disentangling Earth and river tides for an improved estimation of barometric efficiency

The barometric efficiency (BE) based on the $S_2$ response of groundwater levels to atmospheric tides (denoted as $BE_{S_2}^{AT}$) is evaluated to derive the loading efficiency with Eq. (5). Given that $S_2$ signals are present in Earth-, atmospheric-, ocean-, and river tides (Table 1), the $S_2$ signal in groundwater levels $\hat{z}_{S_2}^{GW}$ needs to be disentangled into its individual drivers before $BE_{S_2}^{AT}$ can be calculated. The mathematical formulation of frequency disentanglement uses complex numbers (denotated by a hat) to account for phase and amplitude information of each signal. For settings without river or ocean tide loading, Rau et al. (2020) introduced a complete tidal disentanglement for $S_2$ signals in groundwater heads, based on $M_2$ signals in groundwater and Earth tide strains using:

$$ \hat{z}_{S_2}^{GW\_ET} = \frac{\hat{z}_{M_2}^{GW}}{\hat{z}_{M_2}^{ET}} \hat{z}_{S_2}^{ET} \qquad (6) $$

Recent extensions of this method included OT$_{SE}$ strains in the analysis by a consideration of the total volumetric strain caused by OT$_{SE}$ and Earth tides (Dörr et al., 2025). To achieve this, the two strain time series can be added up to a single time series, further denoted with the superscript {ET + OT$_{SE}$}, prior to the conducted frequency analysis (Dörr et al., 2025).

To enable the calculation of $BE_{S_2}^{AT}$ in deltaic groundwater systems, the present study introduces a novel extension of previous methods (Rau

<caption>Table 1</caption>

Summary of the ten most relevant frequencies in groundwater analysis after McMillan et al., (2019).

<table><thead><tr><th>Darwin name</th><th>Frequency (cpd)</th><th>Earth Tides</th><th>Atmospheric Tides</th><th>Ocean/River Tides</th></tr></thead><tbody><tr><td>Q<sub>1</sub></td><td>0.893244</td><td>Yes</td><td>—</td><td>Yes</td></tr><tr><td>O<sub>1</sub></td><td>0.929536</td><td>Yes</td><td>—</td><td>Yes</td></tr><tr><td>M<sub>1</sub></td><td>0.966446</td><td>Yes</td><td>—</td><td>Yes</td></tr><tr><td>P<sub>1</sub></td><td>0.997262</td><td>Yes</td><td>—</td><td>Yes</td></tr><tr><td>S<sub>1</sub></td><td>1.000000</td><td>—</td><td>Yes</td><td>—</td></tr><tr><td>K<sub>1</sub></td><td>1.002738</td><td>Yes</td><td>—</td><td>Yes</td></tr><tr><td>N<sub>2</sub></td><td>1.895982</td><td>Yes</td><td>—</td><td>Yes</td></tr><tr><td>M<sub>2</sub></td><td>1.932274</td><td>Yes</td><td>—</td><td>Yes</td></tr><tr><td>S<sub>2</sub></td><td>2.000000</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>K<sub>2</sub></td><td>2.005476</td><td>Yes</td><td>—</td><td>Yes</td></tr></tbody></table>
et al., 2020; Dörr et al., 2025), which for the first time allows a quantitative disentanglement of Earth-, atmospheric-, ocean-, and river tide components in groundwater heads. To remove the river tide contribution to the M₂ and S₂ response in groundwater heads, the K₁ signal in river tides and groundwater heads is utilized. From the ratio of K₁ signals in river tides and groundwater, the phase shift and amplitude damping of the K₁ signal are calculated, which are then assigned to the M₂ and S₂ signal of the river water levels to estimate the river tide influence on M₂ and S₂ signals in groundwater observation wells. It is noted that this is a simplified consideration assuming insignificant K₁ contribution of Earth tides to groundwater amplitudes ($\hat{z}_{K_1}^{GW ET} = 0$) and neglecting the frequency dependency of amplitude damping and phase shifting (Valois et al., 2022; 2024). The applicability of this simplified assumption needs to be validated in any individual application.

Similar to the M₂-based disentanglement of Earth tides (Rau et al., 2020, Eq. (6)), the introduced K₁-based disentanglement of river tides can be formulated with Eqs. (7) and (8) for M₂ and S₂ frequencies.

$$ \hat{z}_{M_2}^{GW,RT} = \frac{\hat{z}_{K_1}^{GW}}{\hat{z}_{K_1}^{RT}} \hat{z}_{M_2}^{RT} \qquad (7) $$

$$ \hat{z}_{S_2}^{GW,RT} = \frac{\hat{z}_{K_1}^{GW}}{\hat{z}_{K_1}^{RT}} \hat{z}_{S_2}^{RT} \qquad (8) $$

The M₂ response of groundwater to Earth tides and OT<sub>SE</sub> can then be calculated by subtracting the river tide contribution (Eq. (8)) from the M₂ signal in groundwater:

$$ \hat{z}_{M_2}^{GW \{ET+OT_{SE}\}} = \hat{z}_{M_2}^{GW} - \hat{z}_{M_2}^{GW,RT} \qquad (9) $$

Then, the disentanglement of Earth-, atmospheric-, ocean-, and river tide components in S₂ groundwater signals can be formulated with Eq. (10).

$$ \begin{aligned} \hat{z}_{S_2}^{GW,AT} &= \hat{z}_{S_2}^{GW} - \hat{z}_{S_2}^{GW \{ET+OT_{SE}\}} - \hat{z}_{S_2}^{GW,RT} \\ &= \hat{z}_{S_2}^{GW} - \frac{\hat{z}_{M_2}^{GW \{ET+OT_{SE}\}}}{\hat{z}_{M_2}^{GW \{ET+OT_{SE}\}}} \hat{z}_{S_2}^{\{ET+OT_{SE}\}} - \frac{\hat{z}_{K_1}^{GW}}{\hat{z}_{K_1}^{RT}} \hat{z}_{S_2}^{RT} \end{aligned} \qquad (10) $$

After the disentanglement of Earth-, ocean-, river- and atmospheric tide signals, BE can be calculated based on the S₂ response of groundwater heads to atmospheric tides $BE_{S_2}^{AT}$ as:

$$ BE_{S_2}^{AT} = \frac{1}{A_{S_2}^r} \text{abs} \left[ \frac{\hat{z}_{S_2}^{GW,AT}}{\hat{z}_{S_2}^{AT}} \right] \qquad (11) $$

with $A_{S_2}^r$ as the amplitude response between water levels in the observation well and the actual aquifer pore pressure. $A_{S_2}^r$ is a function of well and aquifer parameters and for wells screened in confined aquifers with $K > 1 \cdot 10^{-5} \text{ m s}^{-1}$, $A_{S_2}^r$ ranges usually between 0.99 and 1.00 and is therefore negligible in Eq. (11) (Rau et al., 2020).

## 2.7. Summary of applied methods and workflow

The deep analysis and quantification of loading imprints in groundwater heads utilizes time series of observed groundwater heads, computed Earth tide- and OT<sub>SE</sub>-strains, measured barometric pressure, river water levels, and precipitation data. A novel method extension is introduced to derive the loading efficiency $\gamma$ from the S₂ response of groundwater heads to atmospheric tides, using K₁ signals for the removal of river tide signals in groundwater head amplitude spectra. The results are verified using two independent established time domain methods: (i) the barometric response function (BRF) within a multi-factor regression deconvolution, and (ii) the groundwater level response to distinct precipitation events. The surface loading dynamics are

assessed using an analysis of precipitation-evaporation-driven accumulation of surface water resources in aquaculture ponds, river water level dynamics and a land use land cover classification. Using multi-factor regression deconvolution as well as complementary the quantification of the loading efficiency $\gamma$ and the derived seasonal loading dynamics, observed groundwater heads are corrected from all loading influences. Fig. 6 presents a flow chart visualizing the study's methodological approach and a concise schematic illustrating the novel tidal disentanglement workflow.

# 3. Results

## 3.1. Seasonal surface water loading

A detailed site-specific analysis of satellite imagery agrees with the regional land use classification of Vu et al. (2022) on a local scale for GW1, identifying inland aquaculture as a dominant LULC-class in the proximity of GW1. Fig. 7a shows that the majority of the land area is covered by aquaculture ponds and vegetated land, while rivers, channels and sealed surfaces have a minor contribution to the overall land cover. Most ponds are characterized by sizes of approx. 1 ha, confirming a dominant share of intensive shrimp farming ponds (Ho et al., 2025). Fig. 7b visualizes the change in percentage for each LULC-class relative to the distance from GW1, showing that aquaculture ponds and vegetated land each comprise approximately 40–50 %, sealed surfaces account for 6–15 % and rivers and channels contribute the remaining 0–5 %, depending on the distance to GW1. The responding area of each groundwater observation well, estimated using ten times the depth of their screen casing to approx. 0.9 km 1.5 km and 2.4 km for GW1a, GW1b and GW1c respectively, is illustrated through dashed vertical lines in Fig. 7b. This shows that for GW1c, the share of ponds to the total responding area is greatest with 50 %, whereas for GW1a and GW1b it is closer to 45 %.

The loading due to aquaculture ponds is calculated by the water balance of observed accumulated precipitation (Fig. 4d and Fig. 8a) and the calculated evaporation for open water bodies (Eq. (2) and Fig. 8a) under consideration of the observed meteorological parameters (Fig. 4). Fig. 8b illustrates the seasonal variation of the calculated pond loading, characterized by a decrease (unloading) between December and May, when evaporation exceeds precipitation and an increase (loading) in the remaining months, which are governed by monsoon precipitation.

The calculated soil moisture loading (Fig. 8b) is characterized by seasonal variations of maximum 0.06 m and is therefore negligible compared to the calculated pond water loading. Thus, even though vegetated land covers 40–45 % of the area (Fig. 7b), the effect of soil moisture variations on the groundwater heads is negligible compared to pond- and river water loading effects. The comparison of the calculated pond water loading and the observed groundwater heads at GW1 (Fig. 8c) reveals a pronounced similarity between the calculated loading and the groundwater head dynamics during the rainy season. Fig. 8c also illustrates a correlation between river water level dynamics and groundwater head dynamics, which is particularly well pronounced during extreme events (peaks in timeseries) as well as in February 2023 (highlighted by a gray box in Fig. 8c) and less in the overall trend.

## 3.2. Groundwater response to distinct precipitation events

The hydrometeorological monitoring station at GW1 recorded distinct precipitation events of 99.0 mm, 289.5 mm and 63.5 mm on 18.06.2023, between 14.07. and 20.07.2023 and on 05.05.2024, respectively. Fig. 9 illustrates the immediate response of groundwater heads in the three observation wells GW1a, GW1b and GW1c to each of these events. Table 2 summarizes the estimated loading efficiency based on the groundwater response to the three precipitation events, calculated as the simple ratio of recorded precipitation and observed groundwater level rise.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Overview of the workflow and methods applied in the present study.</figcaption></figure>

The assessment of groundwater head responses to three distinct precipitation events yields loading efficiency ranges of $\gamma_{GW1a} = 0.87 – 0.90$, $\gamma_{GW1b} = 0.68 – 0.79$ and $\gamma_{GW1c} = 0.64 – 0.72$ with corresponding averages of 0.89, 0.72, and 0.69, respectively.

### 3.3. Loading efficiency from frequency analysis

The loading efficiency was derived in the frequency domain from $BE_{S_2}^{AT}$, i.e., the barometric efficiency based on the $S_2$ response of groundwater heads to atmospheric tides. $BE_{S_2}^{AT}$ was calculated with the novel method extension introduced in this study (Eq.(10)) utilizing the data presented in Fig. 10 comprising volumetric Earth tide strains (Fig. 10a-I), ocean tide loading on the solid Earth (Fig. 10b-I), river

water levels at CM3 Rach Cui (Fig. 10c-I), barometric pressure (Fig. 10d-I) and groundwater heads in three observation wells (Fig. 10e-I to g-I). The FFT amplitude spectra annotated with HALS amplitudes (Fig. 10e-III to g-III) reveal distinct amplitudes in the groundwater heads at $O_1$, $K_1$, $M_2$ and $S_2$ frequency, resulting from a superposition of Earth- and ocean-tide amplitudes (Fig. 10a-III and Fig. 10b-III), river tide amplitudes (Fig. 10c-III) and barometric amplitudes (Fig. 10d-III). The sum of Earth tide- and $OT_{SE}$-strains are dominated by semidiurnal frequencies with a $K_1:M_2$ amplitude ratio of 0.24 while the river tides are dominated by diurnal frequencies with a $K_1:M_2$ amplitude ratio of 2.57. In the three groundwater wells the $K_1:M_2$ amplitude ratio varies between 0.95 (GW1c, 241 m depth), 1.37 (GW1b, 148 m depth) and 1.47 (GW1a, 90 m depth).
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** LULC-classification in the proximity of GW1, Imagery data: Google, © 2024 Maxar Technologies (a) and the change of each LULC-class in radial distance to GW1(b).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Accumulated precipitation and accumulated calculated evaporation at GW1 (a), derived surface loading due to variations in pond water levels and soil moisture at GW1 (b) and dynamics of pond loading, river loading and groundwater heads at GW1 (c).</figcaption></figure>

The amplitude spectrum of GW1a (Fig. 10e-III) shows an additional amplitude peak at 3.0 cpd, a frequency which is not present in the drivers amplitude spectra (Fig. 10a-III to d-III). Therefore, GW1 is not used for the quantification of $BE_{S_2}^{AT}$, given that such an amplitude response at 3.0 cpd is indicative for nearby groundwater extraction afflicting the observed $S_2$ groundwater amplitude (Dörr et al., 2025). For the calculation of the $BE_{S_2}^{AT}$, atmospheric tides were disentangled from Earth-, ocean- and river tides with Eq. (10) for GW1b and GW1c. $A_{S_2}^T$ was calculated as 0.99 for GW1b and GW1c, respectively. Table 3 summarizes the HALS amplitudes, phases, and corresponding standard deviations at the $O_1, K_1, M_2$ and $S_2$ frequencies for multiple drivers, including Earth tides, $OT_{SE}$, river water levels, and barometric pressure, as well as for groundwater heads at GW1a–c. The uncertainties in amplitudes and phases of the tidal constituents in the drivers and groundwater heads are generally small to moderate, with standard deviations mostly below 1 % in amplitude and 1° in phase for Earth tides, $OT_{SE}$, and barometric pressure. Tidal constituents of the river water levels exhibit slightly higher variability, with standard deviations of <4 % and <2.1° for $O_1, K_1$, and $M_2$, and up to 15 % and 8.5° for $S_2$ amplitudes and phases,

respectively. For the groundwater data, GW1a and GW1b show larger relative standard deviations compared to GW1c. At GW1c, standard deviations range between 2.5–5.0 % in amplitude and 1.7–2.9° in phase, whereas GW1a and GW1b exhibit higher relative deviations of 5–19 % in amplitude and 3–11° in phase.

The magnitude-squared coherence between river tides at CM3 and groundwater heads in GW1c yielded high values of 0.91, 0.93, and 0.97 for the major tidal constituents $O_1, K_1$, and $M_2$, respectively, indicating a significant influence of river tides on groundwater heads.

Fig. 11 visualizes the signal disentanglement of the method extension. Note that the radial axes in the polar plots are in logarithmic scale to provide a comprehensive visualization of the wide range of signal amplitudes. The amplitude damping and phase shift between river tides and groundwater signals at $K_1$ (Fig. 11a) are projected to $M_2$ and $S_2$ frequency, to estimate the river tide contribution in the respective groundwater signals (squares in Fig. 11b and c). The $M_2$ groundwater signals are displayed in Fig. 11b after removal of the river contribution. Fig. 11c visualizes the $S_2$ groundwater signals after removal of the Earth tides, $OT_{SE}$ and river tides, i.e. the groundwater response to atmospheric
<figure><img src="image_11.png" /><figcaption>Fig. 9. Hydraulic head response to three precipitation events (a to c) at GW1 site in the observation wells GW1a, $qp_{2-3}$ (I), in GW1b, $qp_1$ (II) and in GW1c, $n_2^2$ (III).</figcaption></figure>

<caption>**Table 2**</caption>

Loading efficiency derived from hydraulic head responses to three precipitation events at GW1 in the southern VMD.

<table><thead><tr><th></th><th>Date &amp; Duration</th><th>GW1a</th><th>GW1b</th><th>GW1c</th></tr></thead><tbody><tr><td>Precipitation (mm)</td><td>18.06.2023</td><td>99.0</td><td>99.0</td><td>99.0</td></tr><tr><td>Groundwater response (mm)</td><td>[6 h]</td><td>88.0</td><td>66.9</td><td>63.8</td></tr><tr><td>Loading Efficiency (–)</td><td></td><td>0.889</td><td>0.676</td><td>0.644</td></tr><tr><td>Precipitation (mm)</td><td>14.-20.07.2023</td><td>289.5</td><td>289.5</td><td>289.5</td></tr><tr><td>Groundwater response (mm)</td><td>[150 h]</td><td>251.6</td><td>228.1</td><td>202.0</td></tr><tr><td>Loading Efficiency (–)</td><td></td><td>0.869</td><td>0.788</td><td>0.698</td></tr><tr><td>Precipitation (mm)</td><td>05.05.2024</td><td>63.5</td><td>63.5</td><td>63.5</td></tr><tr><td>Groundwater response (mm)</td><td>[2 h]</td><td>56.4</td><td>44.3</td><td>44.9</td></tr><tr><td>Loading Efficiency (–)</td><td></td><td>0.903</td><td>0.708</td><td>0.718</td></tr></tbody></table>

tides. These values yield in a calculation of $BE_{S_2}^{AT} = 0.31$ for GW1b and $BE_{S_2}^{AT} = 0.30$ for GW1c with Eq. (11). Consequently, the loading efficiency calculates to $\gamma_{GW1b} = 0.69$ and $\gamma_{GW1c} = 0.70$ with Eq. (5), respectively. For comparison, $BE_{S_2}^{AT}$ was additionally calculated after Rau et al. (2020) without a consideration of OT$_{SE}$ and river tides. This calculation results in $BE_{S_2}^{AT} = 0.29$ for GW1b and $BE_{S_2}^{AT} = 0.30$ for GW1c.

### 3.4. Multi-factor regression deconvolution

Fig. 12 summarizes the response functions to atmospheric tides (Fig. 12a, utilizing data from Fig. 10d-I), (ii) volumetric strains due to Earth tides and OT$_{SE}$ (Fig. 12b, utilizing the sum of timeseries in Fig. 10 a-I and a-II), (iii) pond water loading (Fig. 12c, utilizing data from Fig. 8b) and (iv) river water levels (Fig. 12d, utilizing data from Fig. 10

c-I), obtained from the multi-factor regression deconvolution for GW1a (Fig. 12a-I to d-I), for GW1b (Fig. 12a-II to d-II) and for GW1c (Fig. 12a-III to d-III). This includes an illustration of the calculated uncertainty of each response function. Given that the response functions are the cumulative sum of the respective IRCs, uncertainties of the IRCs accumulate in the response function with increasing time lag consideration.

The response functions of GW1c are characterized by small standard deviations, whereas the multi-factor regression deconvolution yields large standard deviations for GW1a and GW1b. The regression deconvolution yields a good $R^2 = 0.89$ for GW1c. For GW1a and GW1b moderate to poor $R^2$ of 0.55 and 0.16 are calculated within the multi-factor regression deconvolution.

For GW1c, the BRF approaches a constant value of approx. 0.3 (Fig. 12a-III), from which $\gamma_{GW1c} \approx 0.7$ can be derived with Eq. (5). For GW1a, the BRF varies around an average value of approx. 0.1, suggesting $\gamma_{GW1a} \approx 0.9$ while for GW1b, the BRF approaches a value of approx. 0.3 in the first 48 h of time lag consideration and later approx. 0.4, which would suggest $\gamma_{GW1b} \approx 0.6-0.7$. The response function to Earth tide- and OT$_{SE}$-strains (denoted as Earth- and ocean-tide response function EORF) of all three wells (Fig. 12b) approach constant values <0.01.

The response function from the calculated pond water level, denoted as pond loading response function PLRF, are characterized by an initial rise followed by decreasing trends converging towards a constant value for all three observation wells. For GW1a (Fig. 12c-I), the PLRF peaks at approx. 0.9 and decreases to approx. 0.3, whereas for GW1b (Fig. 12c-II), the PLRF decreases from 0.8 to 0.3. The PLRF of GW1c (Fig. 12c-III) is characterized by less noise and a decrease from 0.75 to approx. 0.4. The river water response functions (Fig. 12d) are characterized by
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Employed measurement data (I, left side) and detrended data (II, middle) used for the frequency analysis at GW1 and the resulting FFT amplitude spectrum (lines) and HALS amplitudes (calculated for the 10 frequencies listed in Table 1, displayed for $O_1$, $K_1$, $M_2$ and $S_2$ frequency as dots) (III, right side), for computed volumetric Earth tide strains (a), computed ocean tide loading on the solid Earth (b), river water levels at CM3 Rach Cui (c), barometric pressure (d), groundwater heads in the aquifers $qp_{2-3}$ (e), $qp_1$ (f) and $n_2^2$ (g).</figcaption></figure>

increasing trends, resulting in plateaus of approx. 0.22 for GW1a, 0.30 for GW1b and 0.17 for GW1c, respectively.

For comparison, the BRF was additionally derived for GW1c with MUFACO considering barometric pressure as the only driver, neglecting Earth-, ocean- and river tides as well as pond loading. This comparative analysis resulted in a distinctly different BRF, characterized by an oscillation and increase up to 0.75 and a poor $R^2$ of 0.37. The poor $R^2$ and the resulting high standard deviation in the response functions at GW1a and GW1b are attributable to accumulating uncertainties of the IRC and are a consequence of the low sensor precision of GW1a and GW1b (Fig. 2b).

### 3.5. Estimated correction of groundwater head data

The observed groundwater head data was corrected from loading signals using the calculated response functions for barometric pressure, volumetric strains due to Earth tides and $OT_{SE}$, precipitation-evaporation based loading of ponds and river water levels (Fig. 12). Fig. 13 summarizes the obtained estimation of the corrected

groundwater heads for the three observation wells at GW1, indicating depleting groundwater levels for all wells, when loading signals are removed from the observation data. For GW1a, GW1b and GW1c, the observed groundwater heads are characterized by increasing trends of 1.29 $m a^{-1}$, 1.10 $m a^{-1}$ and 1.24 $m a^{-1}$ respectively during the rainy season, whilst the corrected heads indicate depleting trends of -0.59 $m a^{-1}$, -0.65 $m a^{-1}$ and -0.44 $m a^{-1}$ respectively.

The standard deviations of the response functions (Fig. 12) are projected into the standard deviation of the corrected groundwater levels (Fig. 13). For GW1c (Fig. 13c) the minor standard deviations in the response functions (Fig. 12 a-d III) lead to a moderate standard deviation in the corrected head, encompassing corrected depletion rates ranging from -0.23 $m a^{-1}$ to -0.66 $m a^{-1}$. For GW1a and GW1b, the high standard deviations of the response functions are transferred into high standard deviations for the corrected heads.

For comparison, the groundwater heads at GW1a, GW1b and GW1c were corrected based on the derived loading efficiencies $\gamma$ and the area-averaged surface loading $\Delta\sigma_L$ for each of the three wells. Based on the three applied approaches to calculate the loading efficiency, $\gamma_{GW1a} =$
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Results of the signal disentanglement using Eqs. (7)–(10). The signals comprise of amplitude and phase of $K_1$ river tides in groundwater signals (a), the removal of river tides in $M_2$ groundwater signals (b) and the removal of Earth tides, ocean tide loading on the solid Earth, and river tides (c). Note the logarithmic scale of radial axes in the polar plots.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Response functions obtained from MUFACO for barometric loading (a) Earth tide- and OT<sub>SE</sub>-strains (b) pond water loading (c) and river water loading (d) for the observation wells GW1a (I), GW1b (II) and GW1c (III) under consideration of maximum time lag of 7 days (168 h).</figcaption></figure>

0.9, $\gamma_{GW1b} = 0.7$ and $\gamma_{GW1c} = 0.7$ were utilized for the correction. Given that the accumulation of water resources on vegetated land cannot be quantified precisely with the data available, the surface loading $\Delta\sigma_L$ and therefore also the corrected groundwater heads are characterized by a value range, depending on the considered share of water accumulation on vegetated land between 0 % and 50 %. Fig. 14 shows that for GW1a, GW1b and GW1c the corrected heads based on $\gamma$ and $\Delta\sigma_L$ comply best

with the deconvolution-based correction for a respective share of water accumulation on vegetated land of 25 %, 46 % and 36 %. These shares imply that depending on the distance to GW1, 54 % – 75 % of the vegetated land discharge surface water out of the wells’ responding area, while the remaining fraction underlies the balance of accumulated precipitation and evaporation.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Corrected groundwater heads for GW1a (a) GW1b (b) and GW1c (c) based on a multi factor regression deconvolution under consideration of barometric pressure, Earth tide- and OT<sub>SE</sub>-strains, pond water loading and river water levels under a maximum time lag consideration of 7 days (168 h).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Corrected groundwater heads for GW1a (a) GW1b (b) and GW1c (c) based on derived loading efficiencies and seasonal surface loading dynamics. For comparison, the corrected groundwater heads based on MUFACO (Fig. 13) are included in the graphs.</figcaption></figure>

# 4. Discussion

## 4.1. Quantification of loading signals in deltaic groundwater systems

The loading efficiency $\gamma$ is derived in the frequency domain using a novel extension of previously introduced methods to assess the $S_2$ response of groundwater heads to atmospheric tides. While previous frequency-domain methods were only applicable in the absence of river and ocean tide loading (Acworth et al., 2016; Rau et al., 2020), the method extension presented in this work covers a disentanglement of multiple influences, such as Earth-, OT<sub>SE</sub>-, atmospheric-, and river tides in groundwater heads. The presented $K_1$-based removal of river tide influences on groundwater head amplitude spectra enables the application of $S_2$-based methods for the calculation of $\gamma$ in deltaic groundwater systems, where superposing multiple tidal forcings govern tidal constituents in groundwater heads.

The results of this novel extension are verified using two additional independent methods: (i) the barometric response function (BRF), and (ii) the groundwater head response to distinct precipitation events. All

three methods yield consistent results of $\gamma_{GW1a} \approx 0.9$, $\gamma_{GW1b} \approx 0.7$ and $\gamma_{GW1c} \approx 0.7$, confirming that the new frequency-based approach is valid.

The calculated $K_1:M_2$ amplitude ratio in groundwater heads lie within 0.95–1.47 and therefore between the $K_1:M_2$ amplitude ratios in Earth tide- and OT<sub>SE</sub>-strains (0.24) and the $K_1:M_2$ ratio in river water levels (2.57). This, together with the conducted magnitude-squared coherence analysis, indicates the presence of river tide imprints in groundwater head amplitude spectra and underlines the necessity to disentangle these influences with the introduced novel approach (Eqs. (7)–(10)).

The method's presumption of insignificant $K_1$ contribution from Earth tide- and OT<sub>SE</sub>-strains to groundwater signals is evaluated by a multiplication of the $K_1:M_2$ amplitude ratio in Earth tide- and OT<sub>SE</sub>-strains (0.24) with the $M_2$ amplitude in groundwater due to Earth tide- and OT<sub>SE</sub>-strains (1.7 mm, 1.3 mm, 3.3 mm, see Fig. 11b). This simplified consideration results in calculated amplitudes of 0.4 mm, 0.3 mm, 0.8 mm, which equals a share of 7 %, 8 % and 17 % of Earth tide- and OT<sub>SE</sub>-strain contribution to the $K_1$ signal in groundwater heads. While these contributions are minor compared to the river tide signals, they
might not per se be considered negligible. However, given that the complementary calculation of BE without consideration of river tides and OT<sub>SE</sub> based on the method of Rau et al. (2020) shows minimal variations of BE (0.00–0.02), inaccuracies due to simplifications within this method are negligible for the application with the given tidal composition at the study site. The uncertainties in the HALS-based estimation of tidal constituents (Table 3) are generally low, and given the consistency with loading efficiencies derived from other methods, uncertainties in the tidal constituents are considered to have a negligible effect on the calculated BE<sub>S2</sub><sup>AT</sup> and γ.

The investigation of rainfall loading responses in deep observation wells is an established method for the assessment of loading signals in groundwater head data (Sophocleous et al., 2006; Bardsley and Campbell, 2000; Jan et al., 2007) and is suitable to verify the results of the novel approach that disentangles multiple influences. In assessments of rainfall loading responses, uncertainties may originate from limitation of the groundwater level and the accuracy of the precipitation monitoring instrument. While more intense precipitation events result in a clearer groundwater level response, there is a concomitant decrease in the accuracy of the recorded precipitation. Moreover, intense precipitation events, which typically occur in the tropical monsoon rainfall environment of Southeast Asia, commonly underly high spatiotemporal variation (Mandapaka et al., 2017), highlighting the need for on-site precipitation records when evaluating groundwater head responses to precipitation loading. The minor to moderate variations in the calculated loading efficiency of ± 0.02 to ± 0.06 across three independent precipitation events (Table 2) align well with estimates of γ based on BRF and BE<sub>S2</sub><sup>AT</sup>. In contrast to the first and third precipitation events, which lasted approximately 2 and 6 h respectively, the second event extended over a substantially longer period of about 150 h. Given the derived time factors and the general assessments of Smith et al. (2017), such long-duration precipitation events may allow for pore pressure dissipation through groundwater flow, whereas the shorter events can be considered instantaneous in terms of loading. Accordingly, the two short events are more appropriate for estimating the loading efficiency γ. The long-duration event yields a notably higher γ at GW1b, whereas γ values at GW1a and GW1c remain in the range of those derived from the shorter events. Overall, the results support the method's applicability in the study area, particularly for short-term precipitation events, and underline its potential for application in other deltaic aquifer systems in Southeast Asia with similar hydrometeorological and hydrogeological settings.

Calculating the BRF is a well-established method for the estimation

of BE and γ (Furbish 1991; Rasmussen and Crawford, 1997; Turnadge et al., 2019). The present study derives the BRF within a multi-factor regression deconvolution, and not only from barometric and groundwater observations. The consideration of all relevant loading signals within MUFACO proved essential, given that a comparative analysis considering barometric pressure as the only driver suggested distinctly different values for BE and γ compared to the complete consideration with additional Earth-, ocean- and river tides as well as pond loading. To the authors' best knowledge, this study is first to consider multiple forcings, such as Earth-, ocean- and river tides as well as seasonal surface water accumulation, when calculating BRF, BE and γ. This allows application in more settings such as deltaic aquifers.

Oscillating BRF-patterns as in BRF<sub>GW1a</sub> (Fig. 12a-I) are indicative for non-constant groundwater extraction influencing the regression deconvolution (Haehnel et al., 2024). Dörr et al. (2025) showed that periodic groundwater extraction can leak into amplitude spectra of observation wells at harmonics, i.e., multiples of S<sub>1</sub>, including S<sub>2</sub> and S<sub>3</sub>. Therefore, the presence of the S<sub>3</sub> amplitude peak (Fig. 10d-III) and the oscillating BRF-pattern at GW1a are indicators for nearby groundwater extraction, afflicting the observed S<sub>2</sub> groundwater amplitude.

Overall, the results demonstrate the presence of loading signals in observed groundwater heads across the three aquifers in the southern VMD and establish a solid quantification of γ. This is supported by independent analyses in both the frequency domain (Fig. 10) and the time domain, covering periods from a few hours (Fig. 9a) to days (Fig. 9b) and months (Fig. 8c). Furthermore, the findings align with recent research (Hoang and Steinel, 2021; Dörr et al., 2025) in identifying loading signals in the southern VMD. For the first time, this study provides a quantification of the loading efficiency of γ for three vertically stacked and highly exploited aquifers in the deltaic aquifer of the southern VMD. Further applications in the study area would benefit from more multi-parameter monitoring stations comprising observations of groundwater dynamics, soil moisture and meteorological parameters, such as the here utilized study site GW1.

## 4.2. Seasonal loading characteristics

The dynamics of observed groundwater heads have an evident correlation with the calculated loading due to pond water level dynamics in the rainy season (May to November, Fig. 8c). Based on the elaborated framework of groundwater-loading response for confined groundwater systems (Eq. (3)), the local aquifer-aquitard stratification (Fig. 1b), as well as the time scale of the groundwater-loading responses (Fig. 9,

<caption>**Table 3**</caption>

Calculated HALS amplitudes and phases for multiple drivers and groundwater heads at GW1 for O<sub>1</sub>, K<sub>1</sub>, M<sub>2</sub> and S<sub>2</sub> frequency including the calculated standard deviation (±σ).

<table><thead><tr><th rowspan="3">Frequency</th><th colspan="6">Drivers</th></tr><tr><th colspan="2">Earth tides</th><th colspan="2">OT<sub>SE</sub></th><th colspan="2">River water levels</th></tr><tr><th>Amplitude ± σ (nstr)</th><th>Phase ± σ (°)</th><th>Amplitude ± σ (nstr)</th><th>Phase ± σ (°)</th><th>Amplitude ± σ (mm)</th><th>Phase ± σ (°)</th></tr></thead><tbody><tr><td>O<sub>1</sub></td><td>2.34 ± 0.02</td><td>10.3 ± 0.6</td><td>1.17 ± 0.01</td><td>152.7 ± 0.2</td><td>63.45 ± 1.52</td><td>310.8 ± 1.4</td><td>—</td><td>—</td></tr><tr><td>K<sub>1</sub></td><td>3.42 ± 0.04</td><td>103.7 ± 0.7</td><td>1.95 ± 0.01</td><td>222.6 ± 0.3</td><td>109.2 ± 2.73</td><td>38.4 ± 1.4</td><td>—</td><td>—</td></tr><tr><td>M<sub>2</sub></td><td>17.49 ± 0.02</td><td>22.4 ± 0.1</td><td>1.94 ± 0.01</td><td>183.7 ± 0.1</td><td>42.42 ± 1.52</td><td>201.6 ± 2.1</td><td>—</td><td>—</td></tr><tr><td>S<sub>2</sub></td><td>7.97 ± 0.02</td><td>182.2 ± 0.2</td><td>0.61 ± 0.01</td><td>323.3 ± 0.5</td><td>10.28 ± 1.52</td><td>244.9 ± 8.5</td><td>12.67 ± 0.09</td><td>256.3 ± 0.4</td></tr></tbody></table><table><thead><tr><th rowspan="3">Frequency</th><th colspan="6">Groundwater response</th></tr><tr><th colspan="2">GW1a</th><th colspan="2">GW1b</th><th colspan="2">GW1c</th></tr><tr><th>Amplitude ± σ (mm)</th><th>Phase ± σ (°)</th><th>Amplitude ± σ (mm)</th><th>Phase ± σ (°)</th><th>Amplitude ± σ (mm)</th><th>Phase ± σ (°)</th></tr></thead><tbody><tr><td>O<sub>1</sub></td><td>2.31 ± 0.44</td><td>276.2 ± 10.9</td><td>1.8 ± 0.28</td><td>273.1 ± 8.9</td><td>2.45 ± 0.12</td><td>261 ± 2.9</td><td>—</td><td>—</td></tr><tr><td>K<sub>1</sub></td><td>5.88 ± 0.78</td><td>359.2 ± 7.6</td><td>3.96 ± 0.5</td><td>355.6 ± 7.2</td><td>4.73 ± 0.22</td><td>356.2 ± 2.7</td><td>—</td><td>—</td></tr><tr><td>M<sub>2</sub></td><td>3.99 ± 0.44</td><td>169.5 ± 6.3</td><td>2.89 ± 0.28</td><td>160.1 ± 5.5</td><td>4.97 ± 0.12</td><td>180.2 ± 1.4</td><td>—</td><td>—</td></tr><tr><td>S<sub>2</sub></td><td>8.39 ± 0.44</td><td>92.9 ± 3</td><td>3.61 ± 0.28</td><td>40.2 ± 4.4</td><td>4.12 ± 0.12</td><td>45.1 ± 1.7</td><td>—</td><td>—</td></tr></tbody></table>
Fig. 12) this correlation strongly suggests, that loading effects cause the apparent seasonal rise in groundwater heads.

Directly after precipitation events, the load of the fallen rainwater exerts a spatially extensive and uniform load (Fig. 9 and peaks in Fig. 8c), as well as equally rising river water levels (Fig. 8c). Given that the river water levels regulate the responding area's drainage dynamics, the correlation of increasing and decreasing groundwater head dynamics after intensive precipitation events (peaks in Fig. 8c) with the river water level dynamics is associated to drainage and temporary loading-unloading dynamics. Therefore, the correlation of river water levels and groundwater heads is particularly well pronounced when heads peak and decrease due to distinct intensive precipitation events and the subsequent drainage (Fig. 8c).

The rise in pond water levels, calculated by means of accumulated precipitation and evaporation, complies well with reported seasonally rising water levels in intensive shrimp farming ponds in the study area, as observed by Anh et al. (2010). However, in the absence of direct monitoring data on pond water levels and records of supplemental water use, a more precise characterization of water level dynamics in the aquaculture ponds is not feasible for the investigated study site. Considering legal restrictions on groundwater abstraction for aquaculture (Van Binh et al., 2025) and the declining quality of surface water, the availability of suitable water sources for supplemental use appears limited. The overall good agreement between the precipitation-evaporation-based estimates and both reported seasonal pond water level trends and the loading efficiencies derived from the PLRF and other approaches in this study supports the validity of using precipitation and evaporation as proxies for pond water levels in the study area. Future investigations would benefit from continuous monitoring of pond water levels and quantification of supplemental water use in aquaculture systems. Given the rapid expansion of coastal aquaculture in Southeast Asia (Luo et al., 2022), precipitation accumulation in such ponds likely represents a significant surface loading process influencing groundwater head dynamics in confined deltaic aquifer systems across the region.

Whilst the accumulation of surface water induces an increase in observed groundwater levels during the rainy season (May to November in Fig. 8c), declining surface water levels in the dry season, e.g. due to pond water evaporation and significantly lower river water levels due to reduced upstream discharge of the Mekong River, equally cause a decline of observed groundwater heads. However, due to the lack of longer (i.e., complete annual) datasets, the present analysis is largely restricted to evaluating loading dynamics during the rainy season, which curtails the full potential of the presented approach. Multi-year time series encompassing several seasonal loading and unloading cycles would provide a valuable basis for future investigations aiming to correct groundwater head variations for surface loading effects in the VMD. In the absence of continuous dry-season observations, unloading dynamics can only be addressed qualitatively for short subperiods in this study. For instance, the short temporary rise in observed groundwater levels at GW1 during the dry season in February 2023 shows an evident correlation with the river water level dynamics (highlighted by a square box Fig. 8c). This is interpreted to be due to unloading caused by a temporary hydraulic connection of aquaculture ponds and the multi-channel system, as a result from pond flushing activities, a common practice for sludge removal between production cycles during this time (Anh et al., 2010).

A comparison of the response functions to river water levels (Fig. 12d) and pond water loading (Fig. 12c) shows a stronger response to the calculated loading of ponds than to the river water levels. Together with the evaluation of pond sizes, field observations, and the elaborated interpretation of the PLRF, this supports the assumption that the majority of the aquaculture ponds within the study area are operated as intensive shrimp farming ponds and accumulate water from precipitation.

The identification of surface loading as a major contributor to rising groundwater heads in the southern VMD enables a more refined

evaluation of groundwater responses to periods of extreme hydrological conditions, such as dry-season droughts or exceptionally wet rainy seasons. Within this framework, severe groundwater depletion during droughts can be attributed not only to increased water demand but also to exceptionally low surface water levels and the associated unloading effects. Similarly, rainy seasons with exceptionally high precipitation and the resulting elevated surface water levels induce a more pronounced seasonal rise in groundwater heads compared to rainy seasons with average rainfall.

### 4.3. Aquifer confinement status: Surface water drainage vs. Aquitard leakage

The pond loading response functions obtained from the multi-factor regression deconvolution are characterized by an initial rise and a subsequent decrease with increasing time lag consideration (Fig. 12c). The initial peaks represent the loading efficiency as they feature the groundwater responses to precipitation events, which exert a uniform loading within the entire responding area of the observation well at first. These peak values comply well with the loading efficiency derived from BRF, $BE_{S2}^{AT}$ and the response to the selected distinct precipitation events. Whilst aquaculture ponds are considered to accumulate precipitation and hence the load, urban areas and parts of vegetated land drain the rainwater into the multi-channel systems, from where it is discharged from the responding area of the observation well. This causes a decrease of the overall load, which manifests as a decrease in the PLRF curve. Understanding the decreasing trends of the PLRF curves as unloading signals due to surface water drainage, the values towards which the PLRF converges are associated with the product of the loading efficiency and the percentage of area, which accumulates water resources over the rainy season. For all PLRF curves (Fig. 12c), this complies well with the derived LULC-classification considering ponds and parts of vegetated land (Fig. 7b) and loading efficiencies, in particular for GW1c (Fig. 12c-III), where due to a good sensor resolution, the uncertainty of the response function is minimal. This suggests that the various surface loading of individual LULC can be simplified to a uniform loading based on their share of land cover. However, unless the accumulation of water resources can be formulated precisely for all LULC-classes, the removal of loading signals based on $\gamma$ and $\Delta\sigma_L$ can only provide a range of possible corrected groundwater heads (Fig. 14) depending on the accuracy of $\Delta\sigma_L$ estimates for each LULC-class and the representativeness of the responding area. With regard to the latter, Fig. 7b shows minor variations in the spatial distribution of LULC-class, thereby indicating that imprecise estimates on the size of the responding area would have little impact on the corrected groundwater heads.

While decreasing trends in response functions are commonly indicative for semi-confined conditions (Butler et al., 2011; Hussein et al., 2013), the obtained decreasing trends of the pond loading-response functions are attributed to surface drainage processes in the study area. A misinterpretation of the decreasing trends in response functions as an indicator for aquitard leakage or aquitard drainage would be precarious in the context of land subsidence due to aquitard compaction. Further investigations on the correlation of groundwater depletion and sediment compaction due to aquitard drainage are required, such as depth-differentiated multi-extensometer compaction monitoring (Dörr et al., 2023).

With the available data, an independent quantitative consideration of these drainage-induced temporary loading-unloading processes is not possible beyond the presented descriptive assessment. Further studies, aiming to enhance the assessment of seasonal loading signals by the calculation of a uniform loading time series for the loading response areas of the groundwater wells, need to consider and verify the water balance of all LULC-classes. This would include a quantification of plant and canopy water accumulation as well as multi-year time series encompassing measurements of pond water levels, temporary
inundation levels of vegetated land as well as investigations on the source and use of supplemental water in aquacultural practices.

For barometric as well as Earth- and ocean tide loading signals, which are predominantly characterized by (semi-)diurnal loading cycles, the vertical transient flow of groundwater due to the pore pressure in Eq. (3) may be neglected, as the estimated time factors of 1.4 to 6.7 days are greater than the (semi-)diurnal loading duration (van der Kamp and Schmidt, 2017). This is substantiated by the shape of the BRF and EORF curves (Fig. 12a and b) which are characterized by plateaus and therefore indicate confined conditions, i.e. no vertical flow from an overlying aquitard due to the processes under consideration.

Conversely, for longer loading cycles like the seasonal pond water loading, the vertical transient flow of groundwater may generate a certain flow. The time factor *T* describes the duration in which 93 % of the consolidation and associated drainage processes are completed (Shukla et al., 2009; Terzaghi, 1943). Therefore, using a time lag of 7 days (168 h) is sufficient to capture potential contributions of vertical flow to the groundwater-loading response within the multi-factor regression deconvolution for the study area, where *T* was estimated between 1.4 and 6.7 days. However, as the shape of the PLRF curve is governed by drainage processes of surface water resources, which equally cause a declining trend in the response function, the present study cannot provide evidence for or against aquitard leakage processes based on the shape of response functions.

### 4.4. Corrected groundwater heads and implications for groundwater recharge

Regression deconvolution has been applied in multiple previous studies to remove diverse influences on groundwater heads, e.g., barometric pressure (Furbish 1991; Rasmussen and Crawford, 1997), Earth tides (Toll and Rasmussen, 2007), soil moisture and variations in surface water storage (van der Kamp and Maathuis, 1991), river water levels (Spane and Mackley, 2011), sea levels (Haehnel et al., 2024) or to evaluate aquifer hydraulic conductivity with atmospheric slug test (Valois et al., 2022; 2023).

The multi-factor regression deconvolution applied in this study allows an individual and more flexible consideration of multiple influences (Haehnel et al., 2024) and is therefore suitable for the assessment and removal of multiple forcings on observed groundwater levels as given in deltaic groundwater systems. To the best of the authors' knowledge, this study is the first to consider a variety of natural forcings in the correction of observed groundwater heads, including Earth, ocean and river tides, barometric pressure variations as well as surface water accumulation in anthropogenically managed environments.

In the study area, the removal of loading signals during the rainy season reveals corrected groundwater heads with downward trends, whilst the observed (uncorrected) groundwater heads show rising trends. The estimated corrected depletion rates during the rainy season of -0.59 m a⁻¹, -0.65 m a⁻¹ and -0.44 m a⁻¹ are below the annual average depletion rates of -0.70 m a⁻¹, -0.97 m a⁻¹ and -0.73 m a⁻¹ for the three assessed observation wells calculated between rainy season onsets 2023 and 2024 (Fig. 2a).

The corrected groundwater heads derived from the multi-factor regression deconvolution are subject to uncertainties. The illustrated standard deviations in the results demonstrate that the measurement precision and resolution of groundwater level sensors exert a pronounced influence on the magnitude of these uncertainties. At observation wells GW1a and GW1b, groundwater heads were recorded with low precision (Fig. 2b), resulting in equally imprecise response coefficients. The uncertainties of the response coefficients accumulate in the response function and ultimately manifest in the uncertainties of the corrected groundwater heads. Fig. 12 and Fig. 13 illustrate that both response functions and corrected heads exhibit substantial uncertainties at GW1a and GW1b, where sensor resolution is low, whereas

uncertainties are minimal to moderate at GW1c, where measurements were obtained with higher precision. This pattern is also reflected in the R² of the regression deconvolution, characterized by good R² of 0.89 for GW1c and moderate to poor R² of 0.55 and 0.16 for GW1a and GW1b respectively. Despite their low R² values and high standard deviations, the response functions and corrected heads at GW1a and GW1b are considered reliable, supported by the consistency of groundwater dynamics (Fig. 2a) and response functions (Fig. 12) with those of GW1c, as well as their overall agreement with other methods applied in this study. In addition to uncertainties arising from sensor resolution, uncertainties in the estimation of pond water loading, derived here from a precipitation-evaporation-based approach, contribute to the overall uncertainty of the corrected heads. To reduce these effects in future studies, continuous monitoring of pond water levels is recommended.

Within a cross-validation framework, the performance of the multi-factor regression deconvolution was evaluated exemplarily for GW1c. For this purpose, a set of response functions was derived from 60 sub-datasets, each covering 120 days of the complete observation period, and subsequently applied to correct the full groundwater head time series. Fig. 15 illustrates that all 60 corrected time series consistently confirm the overall depletion trend, albeit with a less pronounced decline compared to the correction based on the full dataset. This result supports the interpretation of persistent groundwater depletion and, at the same time, emphasizes the importance of long-term records comprehensive model validation for a more accurate quantification of the depletion rate.

In light of the outlined uncertainties, independent validation is necessary to strengthen the presented conclusions regarding persistent groundwater depletion. Such validation could involve analyses of mass variations derived from the Gravity Recovery and Climate Experiment (GRACE) satellite mission (Tapley et al., 2004), combined with observations of pond and river water levels and detailed records of groundwater extraction rates at monthly or higher temporal resolution. Considering legal restrictions on groundwater use for large-scale rice farming and aquaculture (Van Binh et al., 2025), it is reasonable to assume that a substantial share of groundwater extraction meets seasonally independent demands, such as domestic water use. Considering the long-term annual-average groundwater depletion rates of recent years and decades together with the conceptual framework presented here,

<figure><img src="image_17.png" /><figcaption>Fig. 15. Measured and corrected groundwater heads for GW1c, including 60 corrected time series derived within the cross-validation framework.</figcaption></figure>
which attributes seasonally rising groundwater heads to surface loading, a persistent decline in loading-corrected groundwater heads throughout the seasons appears realistic. Thus, the smaller yet pronounced depletion rates obtained from the corrected groundwater heads during the rainy season appear reasonable compared to the annual-average groundwater depletion rates. Depending on the underlying uncertainties, the magnitude of these depletion rates may vary, yet the overall pattern is consistently reproduced.

Overall, the evaluation of uncertainties in the corrected groundwater heads indicates that the derived depletion rates may not accurately represent the actual loading-corrected groundwater dynamics. Nevertheless, they support the overall conclusion that loading-corrected groundwater heads decline during the rainy season, whereas the measured data show a rising trend.

In addition to the loading correction by multi-factor regression deconvolution (Fig. 13), groundwater heads were corrected based on the derived loading efficiency $\gamma$ and the change of the area-averaged surface load $\Delta\sigma_L$ (Fig. 14). Due to the uncertainty about surface water accumulation on vegetated land, groundwater heads corrected by this method are characterized by a high uncertainty. A comparison with the groundwater head correction by MUFACO suggests that 25 % (GW1a), 46 % (GW1b) and 36 % (GW1c) of the vegetated land accumulate water resources based on the balance of precipitation and evaporation, while the remaining share of vegetated land is drained. Given that unintended flooding and inundation of agricultural farmland can be counteracted relatively effectively in the VMD by its complex multi-channel system of irrigation and drainage channels (Le et al., 2023), this is in a reasonable order of magnitude and substantiates the correction based on MUFACO.

The comparison of the two approaches for the correction of loading signals underlines that, unless a precise water balance within the loading response area of an observation well can be formulated, a correction of groundwater heads by multi-factor regression deconvolution is preferable to a correction based on calculated values for $\gamma$ and derived surface loading timeseries.

Overall, the evaluation of loading-corrected groundwater heads indicates that groundwater recharge, if present, is insufficient to compensate for extractions in the assessed aquifers. This aligns with the local and regional thicknesses of the Holocene $Q_2$ aquitard of approx. 20–30 m (Fig. 1b; Anderson, 1978; Pham et al., 2019; Bauer et al., 2022), through which percolation of surface water is considered unlikely at the timescale considered in this study. By the multi-factor regression deconvolution, any significant rise in observed heads was attributed to local loading processes on the ground surface. The loading-corrected time series therefore contain no indication for possible recharge at potential aquifer outcrops, which have been shown to propagate as pressure signals through confined aquifers (Armstrong & Narayan, 1998).

In previous studies (Duy et al., 2021), a relationship between groundwater and river water levels in the VMD was identified using time series analysis. However, as the authors acknowledge themselves, the applied black-box correlation is incapable of exploring the causes for the observed relationship. Therefore, considering the present results, the postulation of possible groundwater recharge to shallow and deep aquifers during the flood season should be reconsidered. The present analysis shows that the observed groundwater dynamics can be consistently explained by mechanical loading responses rather than by actual infiltration and recharge. Data-driven approaches can provide valuable insights into process dynamics, but their interpretation must remain embedded within a consistent physical framework. In the context of this study, that framework is the groundwater response to mechanical surface loading. When recharge is postulated, as in Duy et al. (2021), its plausibility should be evaluated in light of the actual hydraulic gradients and transmissivities controlling recharge flow processes.

Furthermore, the applicability of the water-table-fluctuation method must be carefully reconsidered wherever the investigated aquifer is confined or not developed, such as at the monitoring sites Q177 and

Q199 in the southern VMD analyzed by Van et al. (2023). The water-table-fluctuation method attributes rising groundwater heads to groundwater recharge. While this approach is feasible for unconfined aquifers, its application in confined settings is inappropriate and can lead to a misinterpretation of loading-induced head increases as recharge. Under such conditions, the observed rise in groundwater heads in the southern VMD most likely reflects hydro-geomechanical responses to surface loading rather than actual infiltration or recharge.

Overall, this study challenges previous interpretations by Duy et al. (2021) and Van et al. (2023), which attributed groundwater responses to precipitation and river water levels as indicators of recharge without appropriately accounting for surface loading effects.

This study highlights the importance of identifying and removing loading signals from groundwater head data for strongly confined aquifers in aquaculture-dominated deltaic systems to achieve an accurate understanding of water resource dynamics. It emphasizes that reliable groundwater recharge assessments in such environments must account for these loading effects. By doing so, the findings provide valuable guidance for efforts to determine sustainable groundwater extraction rates and preventing further depletion. Future research could combine the presented data-driven approach with physically based methods to further assess potential groundwater recharge and verify both local and regional water balances. In particular, obtaining reliable groundwater extraction data for use in comprehensive water balance analyses, together with a delta-wide hydro-geomechanically coupled groundwater model that accounts for surface loading, would provide a sound framework for more holistic recharge assessments. Such integrative approaches would enable a more detailed understanding of the spatial and temporal variability of actual recharge processes and surface loading, as well as the associated groundwater head dynamics, thereby supporting the validation of recharge estimates derived from the groundwater-loading analysis presented in this study.

By showing that groundwater recharge, if present at all, is insufficient to compensate for the current exploitation, the present study highlights the necessity for developing alternative freshwater resource strategies as essential concepts for the region's sustainable development.

## 5. Conclusion

Observed rising trends in groundwater heads within confined aquifers may not necessarily indicate groundwater recharge, as these trends can also be attributed to mechanical loading at the surface. A comprehensive analysis of loading signals in groundwater heads in the southern VMD between 2023 and 2024 strongly suggests that the seasonal rise in groundwater heads is primarily due to the accumulation of surface water in aquaculture ponds during the rainy season. By removing these loading signals from groundwater observation data, the trend reverses, showing declining groundwater heads. This implies that not only in dry but also in rainy season, groundwater recharge, if present, is incapable of compensating for the groundwater exploitation in the study area. Independent validation, for example using GRACE-based mass balances, is essential to substantiate the conclusions on persistent groundwater depletion.

The influence of temporary flooding and inundation on surface loading processes in deltas has been a persistent phenomenon over the course of geological history. However, the recent rapid expansion of aquaculture ponds has introduced a new infrastructure capable of holding substantial surface water volumes. This shift in land use and land cover is a common feature in many deltaic and coastal zones across Southeast Asia, highlighting the broader relevance of this study beyond the southern VMD. Future research would greatly benefit from monitoring data on pond water level dynamics, solid quantification of hydraulic properties, water balances and other factors contributing to seasonal surface loading to improve the correction of groundwater heads in confined aquifers beneath.

Data-driven approaches to groundwater recharge assessment must
be firmly grounded in a consistent physical framework. Purely statistical
correlations between surface and groundwater levels may lead to mis-
interpretations of recharge processes. The preconditions for applying
easy-to-use methods, such as the water table fluctuation method, must
be carefully evaluated, as these methods are not suitable for confined
aquifer systems. Within the framework of groundwater-loading response
in confined aquifers, the analysis of groundwater head responses to
distinct precipitation events represents a straightforward and effective
diagnostic indicator of loading effects. A sound quantification of surface
and, where present, unconfined aquifer water balances provide an
essential foundation for the (multi-)seasonal correction of groundwater
head time series for loading signals. Combined with a robust definition
of loading efficiency, or integrated within a multi-factor regression
deconvolution, these water balances constitute key input parameters for
isolating and removing loading signals from groundwater head time
series on (multi-)seasonal time scales. Furthermore, the reliability of
such deconvolution-based corrections strongly depends on the avail-
ability of long-term observation records and should be systematically
tested, for example through cross-validation, to ensure robust and
transferable conclusions.

The method of calculating an aquifer's loading efficiency from the groundwater head responses to atmospheric tides at the S₂ frequency has seen significant development over the past decade. Unlike earlier approaches limited to areas without river and ocean tidal influences, the current method incorporates these additional tidal drivers into the analysis, representing a novel advancement. This advancement enables the calculation of loading efficiency using S₂-based methods in deltaic groundwater systems elsewhere, where multiple tidal forcings influence the tidal constituents in groundwater heads.

These findings provide valuable insights for groundwater recharge
assessments and play a critical role in shaping sustainable water
resource management strategies in the southern Mekong Delta and other
confined deltaic aquifers. By evaluating loading signals within ground-
water time series, this study advances a more holistic understanding of
groundwater dynamics — especially the coupling between hydrological
and geomechanical processes — in confined deltaic aquifer systems.

CRediT authorship contribution statement

**Felix Dörr:** Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Jonas Bauer:** Writing – review & editing, Validation, Investigation, Formal analysis. **Gabriel C. Rau:** Writing – review & editing, Validation, Methodology, Conceptualization. **Elias Lewi:** Writing – review & editing, Methodology. **Viet Tran Hoan:** Data curation, Investigation, Writing – review & editing. **Le Thi Mai Van:** Writing – review & editing, Project administration, Data curation. **Remi Valois:** Writing – review & editing, Software. **Anke Steinel:** Writing – review & editing, Conceptualization. **Franz Nestmann:** Writing – review & editing, Supervision, Resources, Project administration, Funding acquisition. **Stefan Norra:** Writing – review & editing, Supervision, Resources, Project administration, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Acknowledgements

This research was conducted in the frame of the project “ViWaT –
Vietnam Water Technologies” funded by the German Federal Ministry of
Education and Research (funding reference: 02WCL1474A). Felix Dörr
thanks the Hector Fellow Academy for support. We gratefully thank the
Vietnamese Ministry of Science and Technology (MOST), the National

Center for Water Resources Planning and Investigation (NAWAPI) and
the Department of Natural Resources and Environment Ca Mau (DONRE
Ca Mau) as well as the Federal Institute for Geosciences and Natural
Resources (BGR) for the good cooperation and providing data. We
acknowledge support by the KIT-Publication Fund of the Karlsruhe
Institute of Technology, Germany. We acknowledge the developers of
the software tools PyGTide, HydroGeoSines, and Groundwater.app for
making their applications and methods publicly available, which sup-
ported parts of the data analysis in this study.

Data availability

The data used in this study are available at a Zenodo repository:
https://doi.org/10.5281/zenodo.17279956.

References

Acworth, R.I., Halloran, L.J.S., Rau, G.C., Cuthbert, M.O., Bernardi, T.L., 2016. An objective frequency domain method for quantifying confined aquifer compressible storage using Earth and atmospheric tides. Geophys. Res. Lett. 43, 11671. https://doi.org/10.1002/2016GL071328.

D.C. Agnew SPOTL: some programs for ocean-tide loading, SIO technical report, Scripps Institution of Oceanography 2012 UC San Diego, California 954322pg (last access: 22 May 2025.

R.G. Allen L.S. Pereira D. Raes M. Smith Crop evapotranspiration: guidelines for computing crop water requirements FAO Irrigation and Drainage Paper 56 1998 Rome, 300 pp.

Anderson, H. R., 1978: Hydrogeologic Reconnaissance of the Mekong Delta in South Vietnam and Cambodia. Contributions to The Hydrology of Asia and Oceania, Geological Survey Water-Supply Paper 1608-R.

Anh, P.T., Kroeze, C., Bush, S.R., Mol, A.P., 2010. Water pollution by intensive brackish
shrimp farming in south-east Vietnam: causes and options for control. Agric. Water
Manag. 97 (6), 872–882. https://doi.org/10.1016/j.agwat.2010.01.018.

Armstrong, D., Narayan, K., 1998. The basics of recharge and discharge. CSIRO, Collingwood, Victoria.

Bardsley, W.E., Campbell, D.I., 2000. Natural geological weighing lysimeters: calibration tools for satellite and ground surface gravity monitoring of subsurface water-mass change. Nat. Resour. Res. 9, 147–156. https://doi.org/10.1023/A:1010147527484.

Bauer, J., Börsig, N., Pham, V.C., Hoan, T.V., Nguyen, H.T., Norra, S., 2022. Geochemistry and evolution of groundwater resources in the context of salinization and freshening in the southernmost Mekong Delta. Vietnam. J. Hydrol. Reg. Stud. 40. https://doi.org/10.1016/j.ejrh.2022.101010.

Bauer, J., Dörr Vu, F., Duong, H.T., Schenk, A., Tran, H.V., Pham, V.C., Börsig, N., Van Der Linden, R., Nguyen, N.H., Eiche, E., Norra, S., 2025. Seawater intrusion in river delta systems. Inter-annual dynamics and drivers of salinity variations in the southern Mekong Delta Vietnam. J. Hydrol. 661, 133745. https://doi.org/10.1016/j.jhydrol.2025.133745.

Binh, D. Van, Tran, D.D., Thái Dương, V.H., Bauer, J., Park, E., Loc, H.H., 2025. Land use change in the Vietnamese Mekong Delta: Long-term impacts of drought and salinity intrusion using satellite and monitoring data. iScience 28. https://doi.org/10.1016/j.isci.2025.112723.

Burgess, W.G., Shamsudduma, H., Taylor, R.G., Zahid, A., Ahmed, K.M., Mukherjee, A., et al., 2017. Terrestrial water load and groundwater fluctuation in the Bengal Basin. Sci. Rep. 7 (1), 3872. https://doi.org/10.1038/s41598-017-04159-w.

Butler, J.J., Jin, W., Mohammed, G.A., Reboulet, E.C., 2011. New insights from well responses to fluctuations in barometric pressure. Ground Water 49 (4), 525–533. https://doi.org/10.1111/j.1745-6584.2010.00768.x.

Davydzenka, T., Tahmasebi, P., Shokri, N., 2024. Unveiling the global extent of land subsidence: the sinking crisis. Geophys. Res. Lett. 51, e2023GL104497. https://doi.org/10.1029/2023GL104497.

Domenico, P.A., Schwartz, F.W., 1998. Physical and chemical hydrogeology, 2nd ed.
Wiley, Chichester, USA.

Dörr, F., Bauer, J., Rau, G.C., Valois, R., Hoan, T.V., Pham, V.C., Van, L.T.M., Steinel, A., Nestmann, F., Norra, S., 2025. Passive subsurface characterization in subsiding deltas: assessing land subsidence mitigation potential with frequency analyses of groundwater heads and superposing harmonic drivers. J. Hydrol. https://doi.org/10.1016/j.jhydrol.2025.133844.

Dörr, N., Schenk, A., Hinz, S., 2024. Land subsidence in the mekong delta derived from advanced persistent scatterer interferometry with an infrastructural reference network. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 1–19. https://doi.org/10.1109/JSTARS.2024.3420130.

Dörr, F., Bauer, J., Tran, H.V., Norra, S., Nestmann, F., 2023. Vietnams Mekong-Delta - Landsenkung infolge von Grundwasserübernutzung. Forsch. | WASSERWIRTSCHAFT 64-68. (In German language) http://doi.org/10.1007/s35147-023-1922-3.

Duy, N.L., Nguyen, T.V.K., Nguyen, D.V., Tran, A.T., Nguyen, H.T., Heidbüchel, I., 2021. Groundwater dynamics in the Vietnamese Mekong Delta: trends, memory effects, and response times. J. Hydrol.: Reg. Stud. 33, 100746. https://doi.org/10.1016/j.ejrh.2020.100746.
Erban, L.E., Gorelick, S.M., Zebker, H.A., 2014. Groundwater extraction, land subsidence and sea-level rise in the Mekong Delta Vietnam. Environ. Res. Lett. 9. https://doi.org/10.1088/1748-9326/9/8/084010.

Furbish, D.J., 1991. The response of water level in a well to a time series of atmospheric loading under confined conditions. In Water Res. 27 (4), 557–568. https://doi.org/10.1029/90WR02775.

Haehnel, P., Rasmussen, T.C., Rau, G.C., 2024. Technical note: Removing dynamic sea-level influences from groundwater-level measurements. Hydrol. Earth Syst. 28 (12), 2767–2784. https://doi.org/10.5194/hess-28-2767-2024.

Healy, R.W., Cook, P.G., 2002. Using groundwater levels to estimate recharge. In Hydrogeology Journal 10 (1), 91–109. https://doi.org/10.1007/s10040-001-0178-0.

Herrera-García, G., Ezquerro, P., Tomás, R., Béjar-Pizarro, M., López-Vinielles, J., Rossi, M., et al., 2021. Mapping the global threat of land subsidence. Science 371 (6524), 34–36. https://doi.org/10.1126/science.abb8549.

Ho, T.Q., Do, H.L., Eggert, H., 2025. Shrimp farming industry in Vietnam: An aquaculture performance indicators approach. In Aquaculture Economics & Management, pp. 1–24. https://doi.org/10.1080/13657305.2024.2449410.

Hoan, T.V., Richter, K.-G., Börsig, N., Bauer, J., Ha, N.T., Norra, S., 2022. An improved groundwater model framework for aquifer structures of the quaternary-formed sediment body in the southernmost parts of the mekong delta Vietnam. Hydrology 9, 61. https://doi.org/10.3390/hydrology9040061.

Hoang, T.H., Steinel, A., 2021. Technical Note TN-IV-01. Groundwater level monitoring data at the U Minh well group September 2016 to October 2020. https://www.deutsche-rohstoffagentur.de/EN/Themen/Wasser/Projekte/abgeschlossen/TZ/Vietnam/tech note IV-01_en.pdf?blob=PublicationFile&v=4.

Hussein, M.E., Odling, N.E., Clark, R.A., 2013. Borehole water level response to barometric pressure as an indicator of aquifer vulnerability. Water Resour. Res. 49 (10), 7102–7119. https://doi.org/10.1002/2013WR014134.

Jacob, C.E., 1950. Flow of groundwater. Engineering Hydraulics 321–386.

Jan, C.-D., Chen, T.-H., Lo, W.-C., 2007. Effect of rainfall intensity and distribution on groundwater level fluctuations. J. Hydrol. 332 (3–4), 348–360. https://doi.org/10.1016/j.jhydrol.2006.07.010.

Jasechko, S., Seybold, H., Perrone, D., Fan, Y., Shamsuddha, M., Taylor, R.G., Fallatah, O., Kirchner, J.W., 2024. Rapid groundwater decline and some cases of recovery in aquifers globally. Nature 625, 715–721. https://doi.org/10.1038/s41586-023-06879-8.

Johnston, D., Trong, N.V., Tien, D.V., Xuan, T.T., 2000. Shrimp yields and harvest characteristics of mixed shrimp-mangrove forestry firms in southern Vietnam: factors affecting production. Aquaculture 188, 263–284. https://doi.org/10.1016/S0044-8486(00)00348-3.

Konikow, L.F., Kendy, E., 2005. Groundwater depletion: a global problem. Hydrolge. J. 13, 317–320. https://doi.org/10.1007/s10040-004-0411-8.

Ledieu, J., De Ridder, P., De Clerck, P., Dautredebe, S., 1986. A method of measuring soil moisture by time-domain reflectometry. J. Hydrol. 88, 319–328. https://doi.org/10.1016/0022-1694(86)90097-1.

Lee, J., Jayakumar, R., Shrestha, S., Han, Z., 2018. Assessment of transboundary aquifer resources in Asia: Status and progress towards sustainable groundwater management. J. Hydrol.: Reg. Stud. 20, 103–115. https://doi.org/10.1016/j.ejrh.2018.01.004.

Le, H.-A., Nguyen, T., Gratiot, N., Deleersnijder, E., Soarez-Frazão, S., 2023. The multi-channel system of the Vietnamese Mekong Delta: Impacts on the flow dynamics under relative sea-level rise scenarios. Water 15 (20), 3597. https://doi.org/10.3390/w15203597.

Luo, J., Sun, Z., Lu, L., Xiong, Z., Cui, L., Mao, Z., 2022. Rapid expansion of coastal aquaculture ponds in Southeast Asia: patterns, drivers and impacts. J. Environ. Manage. 315, 11510. https://doi.org/10.1016/j.jenvman.2022.11510.

Mandapaka, P.V., Qin, X., Lo, E.Y.M., 2017. Analysis of spatial patterns of daily precipitation and wet spell extremes in Southeast Asia. Int Jour of Climatology 37 (S1), 1161–1179. https://doi.org/10.1002/joc.5073.

McMillan, T.C., Rau, G.C., Timms, W.A., Andersen, M.S., 2019. Utilizing the impact of earth and atmospheric tides on groundwater systems: a review reveals the Future potential. Rev. Geophys. 57, 281–315. https://doi.org/10.1029/2018RG000630.

Minderhoud, P.S.J., Erkens, G., Pham, V.H., Bui, V.T., Erban, L., Kooi, H., Stouthamer, E., 2017. Impacts of 25 years of groundwater extraction on subsidence in the Mekong delta Vietnam. Environ. Res. Lett. 12 (6), 64006. https://doi.org/10.1088/1748-9326/aa7146.

Minderhoud, P.S.J., Coumou, L., Erkens, G., Middelkoop, H., Stouthamer, E., 2019. Mekong delta much lower than previously assumed in sea-level rise impact assessments. Nat. Commun. 10, 1–13. https://doi.org/10.1038/s41467-019-11602-1.

Osenbrück, K., Steinel, A., Montcoudiol, N., Manh, V., Bäumle, R., 2025. Geochemical evolution and flow of groundwater impacted by long-term abstraction in the Mekong Delta Vietnam. J. Hydrol. 132881. https://doi.org/10.1016/j.jhydrol.2025.132881.

Pechstein, A., Hoang, T.H., Orilski, J., Le, H.N., Le, V.M., 2018. Detailed investigations on the hydrogeological situation in Ca Mau Province, Mekong Delta. Technical Report No III-5. BGR. Ho Chi Minh. City. https://www.pebs-eu.de/EN/Themen/Wasser/Projekte/abgeschlossen/TZ/Vietnam/tech_repIII-5_en.pdf?blob=publicationFile&v=3.

Penman, H.L., 1948. Natural evaporation from open water, bare and grass. Proc. R. Soc. Lond. Ser. A 193, 120–145. https://doi.org/10.1098/rspa.1948.0037.

Pham, V.C., Bauer, J., Börsig, N., Ho, J., Vu, L., Tran, H., Dörr, F., Norra, S., 2023. Groundwater use habits and environmental awareness in Ca Mau Province, Vietnam: implications for sustainable water resource management. Environ. Challenges 13, 100742. https://doi.org/10.1016/j.envc.2023.100742.

Pham, H.V., Van Geer, F.C., Tran, V.B., Dubelaar, W., Oude Essink, G.H.P., 2019. Paleohydrogeological reconstruction of the fresh-saline groundwater distribution in the Vietnamese Mekong Delta since the late Pleistocene. J. Hydrol.: Reg. Stud. 23, 100594. https://doi.org/10.1016/j.ejrh.2019.100594.

Rasmussen, T.C., Crawford, L.A., 1997. Identifying and removing barometric pressure effects in confined and unconfined aquifers. Groundwater 35 (3), 502–511. https://doi.org/10.1111/j.1745-6584.1997.tb00111.x.

Rasmussen, T.C., Mote, T.L., 2007. Monitoring surface and subsurface water storage using confined aquifer water levels at the savannah river site, USA. Vadose Zone J. 6 (2), 327–335. https://doi.org/10.2136/vzj2006.0049.

Rau, G.C., Cuthbert, M.O., Acworth, R.I., Blum, P., 2020. Technical note: Disentangling the groundwater response to Earth and atmospheric tides to improve subsurface characterisation. Hydrol. Earth Syst. 24, 6033–6046. https://doi.org/10.5194/hess-24-6033-2020.

Rau, G.C., Eulenfeld T., Howe, D., Rietbroek R., Gosselin J.S., Staniewicz S., 2022. hydrogeoscience/pygtide: PyGTide v0.7.1, Zenodo [code].

Schulze, K.C., Kümpel, H.-J., Huenges, E., 2000. In-Situ Petrohydraulic Parameters from Tidal and Barometric Analysis of Fluid Level Variations in Deep Wells: Some Results From KTB. In Ingrid Stober, Kurt Bucher (Eds.): Hydrogeology of Crystalline Rocks. Dordrecht: Springer (Springer eBook Collection, 34), pp. 79–104 https://doi.org/10.1007/978-94-017-1816-5_4.

Schweizer, D., Ried, V., Rau, G.C., Tuck, J.E., Stoica, P., 2021. Comparing methods and defining practical requirements for extracting harmonic harmonic components from groundwater level measurements. Math. Geosci. 53, 1147–1169. https://doi.org/10.1007/s11004-020-09915-9.

Shukla, S., Sivakugan, N., Das, B., 2009. Methods for determination of the coefficient of consolidation and field observations of time rate of settlement — an overview. Int. J. Geotech. Eng. 3 (1), 89–108. https://doi.org/10.3328/IJGE.2009.03.01.89-108.

Smith, C., van der Kamp, G., Arnold, L., Schmidt, R., 2017. Measuring precipitation with a geolsimeter. Hydrol. Earth Syst. Sci. 21 (10). https://doi.org/10.5194/hess-2017-174.

Sophocleous, M., Bardsley, W.E., Healey, J., 2006. A rainfall loading response recorded at 300 m depth: implications for geological weighing lysimeters. J. Hydrol. 319, 237–244. https://doi.org/10.1016/j.jhydrol.2005.06.031.

Spane, F.A., 2002. Considering barometric pressure in groundwater flow investigations. Water Res. 38 (6). https://doi.org/10.1029/2001WR000701.

Spane, F.A., Mackley, R.D., 2011. Removal of river-stage fluctuations from well response using multiple regression. Ground Water 49 (6), 794–807. https://doi.org/10.1111/j.1745-6584.2010.00780.x.

Syvitski, J.P.M., Kettner, A.J., Overeem, I., Hutton, E.W.H., Hannon, M.T., Brakenridge, G.R., et al., 2009. Sinking deltas due to human activities. Nature Geosci. 2 (10), 681–686. https://doi.org/10.1038/ngeo6029.

Tapley, B.D., Bettadpur, S., Watkins, M., Reigber, C., 2004. The gravity recovery and climate experiment: mission overview and early results. Geophys. Res. Lett. 31, 2004GL019920. https://doi.org/10.1029/2004GL019920.

Terzaghi, K., 1943. Theoretical soil mechanics. Wiley, Chichester, UK.

Toll, N.J., Rasmussen, T.C., 2007. Removal of barometric pressure effects and earth tides from observed water levels. Ground Water 45 (1), 101–105. https://doi.org/10.1111/j.1745-6584.2006.00254.x.

Turnadge, C., Crosbie, R.S., Barron, O., Rau, G.C., 2019. Comparing methods of barometric efficiency characterization for specific storage estimation. Groundwater 57 (6), 844–859. https://doi.org/10.1111/gwat.12923.

Valiantzas, J.D., 2006. Simplified versions for the Penman evaporation equation using routine weather data. J. Hydro. 331 (3–4), 690–702. https://doi.org/10.1016/j.jhydrol.2006.06.012.

Valois, R., Rivière, A., Vouillamoz, J.M., Rau, G.C., 2024. Analytical solution for well water response to Earth tides in leaky aquifers with storage and compressibility in the aquitard. Hydrol. Earth Syst. 28 (4), 1041–1054. https://doi.org/10.5194/hess-28-1041-2024.

Valois, R., Derode, B., Vouillamoz, J.M., Valerie Kotchoni, D.O., Lawson, M.A., Rau, G.C., 2023. Use of atmospheric tides to estimate the hydraulic conductivity of confined and semi-confined aquifers. Hydrgeol. J. 31 (8), 2115–2128. https://doi.org/10.1007/s10040-023-02715-5.

Valois, R., Rau, G.C., Vouillamoz, J.M., Derode, B., 2022. Estimating hydraulic properties of the shallow subsurface using the groundwater response to Earth and atmospheric tides: a comparison with pumping tests. Water Resour. 58 (5), e2021WR031666. https://doi.org/10.1029/2021WR031666.

Van Tuan, P., Jiang, Y., Stigter, T., Zhou, Y., 2024. Understanding groundwater use and vulnerability of rural communities in the Mekong Delta: the case of Tra Vinh province Vietnam. Groundwater for Sustainable Development 25, 101095. https://doi.org/10.1016/j.gsd.2024.101095.

Van, L.T.T., Lertsirivorakul, R., Vuong, B.T., Thang, C.H., 2023. Groundwater recharge in Mekong River Delta: an application of the water-table fluctuation method in the Long Xuyen Quadrangle and the Ca Mau Peninsula. Songklanakarin J. Sci. Technol. 45 (1), 51–58.

van der Kamp, G., 1972. Tidal fluctuations in a confined aquifer extending under the sea, 24th International Geological Congress, Section 11 (pp. 101–106). Montreal, Que., Canada.

van der Kamp, G., Gale, J.E., 1983. Theory of earth tide and barometric effects in porous formations with compressible grains. Water Resour. Res. 19 (2), 538–544 ISSN 19447973.

van der Kamp, G., Schmidt, R., 2017. Review: Moisture loading—the hidden information in groundwater observation and well records. Hydrgeol. J. 25 (8), 2225–2233. https://doi.org/10.1007/s10040-017-1631-z.
van der Kämp, G., Maathuis, H., 1991. Annual fluctuations of groundwater levels as a result of loading by surface moisture. J. Hydrol. 127, 137–152. https://doi.org/10.1016/0022-1694(91)90112-U.

Vu, H.T.D., Tran, D.D., Schenk, A., Nguyen, C.P., Vu, H.L., Oberle, P., et al., 2022. Land use change in the Vietnamese Mekong Delta: new evidence from remote sensing. Sci. Total Environ. 813, 151918. https://doi.org/10.1016/j.scitotenv.2021.151918.

Wagner, F., Tran, V.B., Renaud, F.G., 2012. Chapter 7, Groundwater resources in the Mekong Delta: availability. In: Renaud Fabrice G., K.C. (Ed.), The Mekong Delta System - Interdisciplinary Analysis of a River Delta, 1st ed. Springer, Dordrecht, pp. 201–220. https://doi.org/10.1007/978-94-007-3962-87.

Wenzel, H.-G., 1996. The nanoGal software: Earth tide data processing package: ETERNA 3.3. Bulletin D'informations Des Marées Terrestres 124, 9425-9439.