Groundwater recharge processes in an Asian mega-delta: hydrometric evidence from Bangladesh

Groundwater is used intensively in Asian mega-deltas yet the processes by which groundwater is replenished in these deltaic systems remain inadequately understood. Drawing insight from hourly monitoring of groundwater levels and rainfall in two contrasting settings, comprising permeable surficial deposits of Holocene age and Plio-Pleistocene terrace deposits, together with longer-term, lower-frequency records of groundwater levels, river stage, and rainfall from the Bengal Basin, conceptual models of recharge processes in these two depositional environments are developed. The representivity of these conceptual models across the Bengal Basin in Bangladesh is explored by way of statistical cluster analysis of groundwater-level time series data. Observational records reveal that both diffuse and focused recharge processes occur in Holocene deposits, whereas recharge in Plio-Pleistocene deposits is dominated by indirect leakage from river channels where incision has enabled a direct hydraulic connection between river channels and the Plio-Pleistocene aquifer underlying surficial clays. Seasonal cycles of recharge and discharge including the onset of dry-season groundwater-fed irrigation are well characterised by compiled observational records. Groundwater depletion, evident from declining groundwater levels with a diminished seasonality, is pronounced in Plio-Pleistocene environments where direct recharge is inhibited by the surficial clays. In contrast, intensive shallow groundwater abstraction in Holocene environments can enhance direct and indirect recharge via a more permeable surface geology. The vital contributions of indirect recharge of shallow groundwater identified in both depositional settings in the Bengal Basin highlight the critical limitation of using models that exclude this process in the estimation of groundwater recharge in Asian mega-deltas.


Introduction
Thick sequences of unconsolidated sediments occur in 'megadeltas' throughout Asia and comprise the termini of major drainage basins that include the rivers Ganges-Brahmaputra-Meghna (GBM), Indus, Irrawaddy, Chao Phraya, Mekong, Red (Song Hong), Pearl (Zhujiang), Yangtze (Chiangjiang) and Yellow (Huanghe; Fig. 1). Coarser sediments within Quaternary fluvial, alluvial, and estuarine sequences form productive aquifers from which shallow groundwater is drawn intensely for dry-season irrigation (Alauddin and Quiggin 2008;Shamsudduha et al. 2011Shamsudduha et al. , 2012 as well as rural and urban domestic water supplies (Ravenscroft et al. 2005;Zahid and Ahmed 2006;Hoque et al. 2007;Naik et al. 2008;Hardoy et al. 2013). The sustainability of groundwater withdrawals and the processes by which groundwater is replenished in Asian mega-deltas, remain inadequately resolved. This knowledge gap constrains not only the modelling of groundwater recharge in response to global change but also our understanding of the vulnerability of shallow groundwater resources to contamination (e.g., Fendorf et al. 2010;Burgess et al. 2011;van Geen et al. 2013) and depletion (e.g., Shamsudduha et al. 2009Shamsudduha et al. , 2012Bui et al. 2012).
There remains a surprising dearth of studies examining the processes of groundwater recharge in Asian mega-deltas and deltaic environments more widely (Sawyer et al. 2015). Diffuse recharge resulting from the direct infiltration of rainfall at the land surface is commonly assumed in humid environments (de Vries and Simmers 2002;Healy and Scanlon 2010). Focused recharge resulting by way of leakage from surface-water bodies including ponds and ephemeral or perennial rivers is often considered to be the dominant process of groundwater recharge in drylands (Scanlon et al. 2006;Cuthbert et al. 2019). It has, however, also been traced to occur in the humid GBM delta around Dhaka City using Cl:Br ratios (Hoque et al. 2014), urban pollutants (Burgess et al. 2011), and stable isotope ratios of O and H (Darling et al. 2002). Exchanges between groundwater and surface water in deltaic environments have also been observed in response to hydrologic events including hurricanes in the lower Mississippi Delta (Li and Tsai 2020) and monsoonal flooding in the Okavango Delta (Wolski and Savenije 2006) and lower Amazon floodplain (Rudorff et al. 2014).
Substantial uncertainty remains, however, in the proportions of groundwater recharge in the GBM delta arising from diffuse and focused recharge pathways, and what factors control these processes. In distal areas of the GBM delta, seasonal oscillations in deep (>150 m below ground level, bgl) wells are attributed to the impact of seasonal surface water, and soilwater loading and unloading on groundwater levels rather than groundwater recharge and discharge (Burgess et al. 2017). Indeed, more recent work in distal areas of the GBM delta (Woodman et al. 2019) suggests that such poroelastic responses could explain seasonal groundwater-level oscillations at depths as shallow as 30 m bgl.  (Ravenscroft and Rahman 2003) and the major physiographic features (e.g. Brahmaputra floodplains, Madhupur Tract) along the transect line Hydrological impact of developing shallow groundwater Bengal Basin within the Ganges-Brahmaputra-Meghna (GBM) delta system, also known as the Bengal Mega-delta in Bangladesh and West Bengal (India), is the largest of the Asian mega-deltas covering an area of more than 100,000 km 2 -Department of Public Health Engineering (DPHE) 2001; Akter et al. 2016). Sedimentary sequences of Quaternary sands, silts, and clays are several kilometres thick (Ravenscroft et al. 2005). In Bangladesh, historic dependence upon surface waters switched to groundwater in the early 1980s with a rapid expansion in the use of shallow groundwater (<150 m bgl) through water wells (i.e., tubewells) to expand access to safe water and enable dry-season irrigation to improve food security. By 2007-2008, almost 80% of all water used in irrigation was supplied by groundwater of which about 80% is derived from shallow tubewells and 20% from deep (>150-300 mbgl) tubewells (Bangladesh Bureau of Statistics 2012).
The dramatic increase in shallow groundwater withdrawals in the Bengal Mega-delta has greatly influenced seasonal oscillations in groundwater levels and amplified previously low vertical hydraulic gradients that had persisted for millennia (Goodbred and Kuehl 2000;DPHE 2001). Where surface geologies are permeable, dry-season groundwater abstraction for irrigation has been observed to induce greater recharge during the subsequent monsoon (Shamsudduha et al. 2011(Shamsudduha et al. , 2015 through increased capture of surface-water flows (i.e. focused recharge) in a manner previously described as the "Ganges Water Machine" (Revelle and Lakshminarayana 1975) and increased vertical hydraulic gradients amplifying diffuse recharge. In contrast, where surface geologies are comparatively impermeable, intensive groundwater withdrawals are causing a net depletion in groundwater storage, which was recently estimated to be~1 km 3 year −1 in Bangladesh (Shamsudduha et al. 2012).

Quantification of recharge in the Bengal Basin
The quantification of groundwater recharge in the Bengal Basin has, to date, relied upon tools or models that assume recharge derives solely from the direct infiltration of rainfall-e.g., United Nations Development Programme (UNDP 1982); Master Plan Organization (MPO 1987). Indeed, the estimation of potential recharge using a lumped-parameter, soil-water balance model (MPO 1987) has been deemed conceptually appropriate (Saleh and Nishat 1989) and later revised by Water Resources Planning Organisation (2000a). Master Plan Organisation (1987) considered 'usable recharge' simply to be 75% of the 'potential recharge'. Quantified recharge via this approach has, however, been shown to differ substantially from observed recharge in Bangladesh computed using the water-table fluctuation method (Shamsudduha et al. 2011). This latter analysis, which is independent of recharge pathways or processes, estimates an average net groundwater recharge to bẽ 150 mm year −1 (σ = 140 mm year −1 ) for the period of 1985(Shamsudduha et al. 2011). Despite suggestions from previous studies that focused recharge can be an important recharge pathway in the Bengal Basin (DPHE 2001;Ravenscroft and Rahman 2003;Zahid et al. 2008;Burgess et al. 2011), locally developed recharge estimation tools and hydrological models (Kirby et al. 2013(Kirby et al. , 2015Ahmad et al. 2014) as well as large-scale models such as PCRaster GLOBal Water Balance model (Wada et al. 2010), are restricted to the estimation of diffuse recharge that results solely from the direct infiltration of precipitation at the soil surface.
This paper explores fundamental questions about the nature of monsoonal recharge in a mega-delta environment: What are the dominant processes by which recharge occurs? How do these processes vary under different depositional environments and soil lithologies? How has intensive groundwater withdrawals influenced recharge dynamics? This exploratory research exploits hydrometric evidence as the Bengal Megadelta of Bangladesh (~100,000 km 2 ) that features networks of dedicated monitoring stations for daily rainfall (~300), weekly groundwater-levels (~1,200), and subdaily to daily river stage (~300); groundwater-level records, managed by the government's agency, Bangladesh Water Development Board (BWDB), have been collected since early 1970s. Here, rare hourly in-situ measurements of colocated groundwater-level and rainfall records are examined together with a much larger body of lower-frequency (daily to weekly) measurements of rainfall, surface water, and groundwater levels. The latter involves regionalizing hydrologic phenomena through a cluster exercise that distinguishes physical patterns in groundwater time-series hydrographs. In-situ analyses build upon a limited number of past studies employing short-term but highfrequency groundwater-level records from localised observations at two locations in Bangladesh (Harvey et al. 2006;Stute et al. 2007;Aziz et al. 2008;Burgess et al. 2017).

Shallow groundwater regimes of the Bengal Basin
The surficial geology of Bangladesh ( Fig. 1) is primarily characterised by Holocene unconsolidated sediments that cov-er~80% of the land surface. The remainder comprise pre-Holocene sediments that occur primarily in terraces as well as hills located in the eastern part of Bangladesh (Fig. 1). Plio-Pleistocene terrace deposits, located in northwest (Barind Tract) and north-central (Madhupur Tract) Bangladesh (Alam et al. 1990;Reimann 1993), occupy 8% of the nation's land surface and feature annual groundwater withdrawals of ~9.7 km 3 , which is~30% of the total groundwater withdrawals (32 km 3 ) estimated for Bangladesh . Both Quaternary deposits consist of sands, silts, and clays that form aquifers throughout the basin. Groundwater generally occurs at shallow depth (<10 m bgl) over most of Bangladesh within Holocene alluviums, alluvial fan deposits, floodplains, and delta plains but at comparatively deeper depths (>10 m bgl) in Plio-Pleistocene fluvio-deltaic sediments of the Madhupur and Barind Tracts in the Bengal Aquifer System (Ahmed et al. 2004;Shamsudduha et al. 2011). These two primary hydrogeological regimes or typologies, namely the Holocene and Plio-Pleistocene, feature surficial covers with contrasting permeabilities, relatively higher and lower, respectively. Holocene aquifers are commonly unconfined or semiconfined with higher transmissivity (2,000-5,000 m 2 day −1 ), whereas Plio-Pleistocene aquifers are confined to semiconfined with lower transmissivity (300-3,000 m 2 day −1 ) (DPHE 2001;Ravenscroft and Rahman 2003). Plio-Pleistocene aquifers typically occur below a clay unit, stratigraphically known as the Madhupur Clay Formation (Fig. 1c), that varies in thickness from~8 to 45 m; the underlying aquifer is known as Dupi Tila (Ravenscroft et al. 2005;Shamsudduha and Uddin 2007). It is noteworthy that lithologies within Holocene sediments are highly variable both in horizontal and vertical directions producing localised (from several hundred metres up to tens of kilometres), highly heterogeneous aquifers at very shallow depths (<50 m bgl) throughout Bangladesh (DPHE 2001).

Sites of high-frequency (hourly) monitoring
High-frequency monitoring sites were established in central Bangladesh (Fig. 1b) as weekly groundwater-level observations are unable to represent diurnal processes and signals including pumping regimes, atmospheric/tidal effects, and evapotranspiration (Acworth et al. 2015). Monitoring stations were constructed in two contrasting surface-geological units: (1) the first site (latitude: 24.468°N and longitude: 89.875°E) in Bhuapur Upazila (subdistrict) of Tangail District located in the Holocene deposits (Fig. 2a), hereafter referred to as Bhuapur, and (2) the second site (latitude: 23.879°N and longitude: 90.274°E) in Savar Upazila (subdistrict) of Dhaka District located in the Plio-Pleistocene deposits (Fig.  2b), hereafter referred to as Savar. Bhuapur is located in a rural setting where the landscape is dominated by irrigated agricultural lands with a few rural settlements. Savar is proximate to Dhaka (20 km) and the periurban town of Savar (4 km). The landscape surrounding Savar is dominated by manufacturing industries, garment factories, business facilities, and settlements. At Bhuapur, groundwater is predominantly used for dry-season rice cultivation, and domestic use, and supplied respectively through shallow irrigation pumps and handoperated tubewells. Boro (dry-season: December to April) and transplanted Aman (wet-season: July to November) rice are produced in >80% of the total cultivable land in Bhuapur Upazila. In contrast, groundwater in and around the site at Savar is also widely used for industrial, municipal, and domestic water supplies throughout the year.
Rainfall monitoring stations of the Bangladesh Meteorological Department (BMD) are located at distances of 15 km (BMD Headquarters in Dhaka) and 25 km (BMD Station in Tangail Town) from Savar and Bhuapur, respectively (Fig. 2). Summary statistics of long-term hydro-meteorological observations (i.e. groundwater levels, river-levels, rainfall) are presented in Table 1 and provide a broad hydrogeological context for high-frequency records collected in Bhuapur and Savar over a period of two hydrological years (2009)(2010)(2011). In addition, lithological logs recorded during monitoring-well (i.e. borehole) installation-see Fig. S1 in the electronic supplementary material (ESM)-illustrate the contrasting hydrogeological conditions at these two sites. At Bhuapur, located within Holocene deposits, groundwater levels are commonly shallow (<5 m bgl), whereas at Savar, groundwater levels within the Plio-Pleistocene deposits are typically deeper (>10 m bgl) and overlain by a thick clay layer (i.e. Madhupur Clay Formation).
At Bhuapur, the surficial geology comprises mainly Young and Old Brahmaputra Floodplain deposits according to the National Geological Map of Bangladesh (Alam et al. 1990). The Brahmaputra floodplain sediments feature a complex relief of broad and narrow ridges, interridge depressions, partially filled cut-off channels, and localised depressions. Bhuapur is covered by permeable silt loam to silty-clay loam soils on ridges and impermeable clays at low-land depressions. Soils predominantly comprise grey (unoxidised) floodplain soils with a generally thin and permeable surficial silt and clay unit (Water Resources Planning Organisation 2000a). At Bhuapur, near-surface lithology is predominantly light-brown sand with some silt having a thickness of~3 m (Fig. S1a of the ESM). Below the surficial sand layer there is a thin (<1 m) silty-clay layer and the shallow aquifer occurs just below that layer. The shallow aquifer comprises medium gray sand with little fine sand. Shallow groundwater at Bhuapur is generally unconfined with dry-season groundwater levels occurring at <8 m bgl and wet-season levels at~1 m bgl; the screen depth interval is 15-18 m bgl.
The site in Savar is located on the southern side of the Madhupur Tract, which is slightly more elevated than the surrounding floodplains and less prone to seasonal flooding. The main soil type of the Madhupur Tract is red-brown coloured lateritic terrace soil rich in clay (Fig. S1b of the ESM). Consequently, the vertical hydraulic conductivity of the Madhupur Clay is low, ranging from~10 − 3 to 10 −2 m day −1 (Michael and Voss 2009;Hassan and Zahid 2011); direct rain-fed recharge to the underlying aquifer (i.e. Dupi Tila sand) is restricted. The National Hydrochemical

Hydrometric observations
High-frequency observations of groundwater levels and rainfall were established in 2009. Both sites were equipped with automatic data loggers that include a groundwater-level logger (Solinst Gold LT M30 at Bhuapur and LT M100 at Savar), a barometric-level logger (Solinst Gold Baro-logger M1.5), and an automated logger (AR-DT2) installed with a tippingbucket rain-gauge (ARG100); rain-loggers at both sites failed in June 2010. Half-hourly groundwater-level data for the period of February 2009 to April 2011 were recorded at both Bhuapur and Savar and converted to hourly data. Proximate to both locations, daily rainfall records for the period of 1987 to 2014 were collated from two BMD stations (Tangail and Dhaka) along with hourly surface-water level (SWL) records from a Bangladesh Inland Water Transport Authority (BIWTA) gauging station in Bhuapur (BIWTA-2730). Additionally, daily records from three BWDB river-gauging stations (SW342 on River Jhenai near Bhuapur; SW14.5 on River Bangshi and SW69 on River Dhaleshwari, the most proximate gauges to Savar) and weekly time-series records (1987 to 2014) at four BWDB boreholes (TA013, TA014, DH073, and DH115) were collated to investigate long-term trends and seasonality and to compare to high-frequency observations (Fig. S2 of the ESM) at both monitoring stations. Fourier frequency analysis in R environment platform (R Core Team 2016) was conducted on the high-resolution records of barometric pressure, groundwater levels, and riverstage levels in order to identify any dominant signals in the time-series data. Commonly used extreme rainfall indices used in Bangladesh by several authors Nowreen et al. 2015) were computed by using R package called 'Rclimdex' (Zhang and Yang 2004).
Atmospheric or barometric pressure (P atm ) was removed from measured total pressure (P tot ) by level-loggers in order to derive observed water levels (W obs = P tot -P atm ) in monitoring boreholes. Despite heavy influence of pumping on groundwater-level observations at both sites, barometric efficiency (BE) was calculated from the slope of a linear trend line fitted through changes in measured water levels (ΔW obs ) and changes in concurrent atmospheric pressure (ΔP atm ) over the same period of time (Burgess et al. 2017). BE is used here simply as an indicator of aquifer conditions (e.g., unconfined versus confined) and to consider potential poroelastic responses in measured groundwater levels.

Statistical clustering of national groundwater-level observations
The representivity of observations at Bhuapur and Savar was explored using evidence from a cluster analysis of groundwater-level records across the Bengal Basin in Bangladesh. These national-scale datasets were thoroughly checked for spurious records; all statistical analyses were conducted in R programming language platform (R Core Team 2016). Initially, quality assured longer records of 786 monitored wells were picked out of a total of 1265 wells monitored by the BWDB. Station records with less than 10% missing values were selected for cluster analyses where no well record had two consecutive missing years. This process resulted in a time series (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) of 464 stations available for analysis.
Statistical clustering is an unsupervised machinelearning method for partitioning datasets into a set of groups or "clusters". Unlike feature-based groupings of groundwater-level times series (e.g., Vernieuwe et al. 2007;Bloomfield et al. 2015), this study applies instance-based time series groupings that characterise the amplitude of trends and fluctuations in groundwater-level times series. Clustering algorithms can be broadly classified into non hierarchical (partitioning method) k-means and hierarchical (i.e. agglomerative or divisive) methods. The k-means clustering procedure splits a set of objects into a selected number of groups by maximizing between-cluster variation and minimizing within-cluster variation. This method is efficient for large datasets. In hierarchical clustering, initially, each object (i.e., case or variable) is considered as a separate cluster. Clusters are combined based on their "closeness" using one of many possible definitions until all objects are combined into one single cluster. Here, 'closeness' was computed using linkage methods in the NbClust package of R platform as described by Charrad et al. (2014). Finally, three variables-groundwater abstraction, lithology of the soil (i.e., upper silt-clay thickness), and local elevation represented by inundation land types (i.e., flood phases as described in Brammer 1988)-are considered collectively to present a process-based description of the clusters. Three ranges (low, moderate, high) for these variables were considered in the interpretation of clusters.

Long-term groundwater-level hydrographs in Bhuapur and Savar
Weekly time-series (1987-2014) records illustrate long-term patterns in BWDB groundwater-level (GWL) fluctuations at Bhuapur and Savar (Fig. 3); details of these BWDB monitoring boreholes are provided in Table 1 and Fig. S1 of the ESM. Two hydrographs (TA013 and TA014) within Holocene deposits in Bhuapur show no discernible change in seasonality and trend in long-term mean GWL fluctuations (Fig. 3a). Both hydrographs reveal interannual variations with some years being wetter (e.g., 1998, 2007) or drier (e.g., 1994, 2006, 2012) compared to mean GWL fluctuations. Median  annual GWL fluctuations (difference between annual minima and maxima) in Bhuapur and Savar are~5.5 and 5.0 m, respectively. Weekly monitored GWLs (TA013 and TA014) and hourly (BIWTA-2730) river levels at Bhuapur show fluctuations of similar magnitude throughout the period of monitoring.
At Savar, two hydrographs (DH073 and DH115), located within the Pleistocene Madhupur Tract, show a dynamic temporal response over nearly three decades (1987Fig. 3c). GWLs at both sites were dominated by strong seasonality until early 2004-2005 when GWLs in borehole DH115 start to decline with progressively reduced seasonality; a decline in GWLs at borehole DH073 occurred later (since 2009-2010) and more gradually. Mean annual GWL fluctuations in DH073 and DH115 from 1990 to 2000 are~4.5 and~7.0 m, Mean monthly records in Fig. 4 depict seasonal regimes in GWLs, rainfall, and river stage (SWLs) at Bhuapur and Savar. At both Bhuapur and Savar, GWLs are lowest in March-April at the end of dry-season during Boro rice cultivation, when seasonal rainfall is low (~100 mm month −1 ), whereas GWLs rise rapidly in May and June, reaching peak levels in July (TA014) and August (TA013) at Bhuapur, and in September (DH115) and October (DH073) at Savar. SWLs rise earlier and reach peak stages in July at Bhuapur and August at Savar, suggesting an approximate 1 month lag between peak GWLs and SWLs in both areas.

Water-level dynamics in Holocene deposits: Bhuapur
Hourly monitoring of shallow GWLs at Bhuapur covers two complete annual cycles (February 2009 to April 2011), showing GWL responses to heavy rainfall events and large SWL fluctuations (Fig. 5). An analysis of daily rainfall records (1987-2014) from a BMD station (Tangail Town) near Bhuapur reveals that the 95th percentile and 50th percentile of 1-day maximum rainfall vary from 222 to 120 mm day −1 respectively (Table S1 of Fig. S3 of the ESM. These individual, daily responses are partially masked during the monsoon period by the seasonal rise in GWLs. It is noteworthy that variations in peak GWLs approximately coincide with peak river stages in the adjoining River Brahmaputra (Fig. 5a). High-frequency monitoring of individual peaks in SWLs (BIWTA-2730) and GWLs over the monsoon season reveals, however, that these are not synchronous (Fig. 5) but feature lags of 0.5-7.5 days.
Following the end of the monsoon in October, river stage (BIWTA-2730) initially recedes more rapidly than the GWL at Bhuapur. A sharp increase in the GWL recession occurs at the end of the year (December) and is associated with a rise in diurnal GWL oscillations that continue until the start of the following monsoon (Fig. S4 of the ESM). Fourier frequency analysis identifies responses to four atmospheric signals (Fig.  S5 of the ESM): (1) once daily (~24 hourly) S 1 or solar diurnal, (2) twice daily (~12 hourly) S 2 or solar semidiurnal, (3) thrice daily (~8 hourly) S 3 or solar terdiurnal, and (4) S 4 or solar quarter-diurnal. Of note is the contrast between the higher amplitude S 2 signal in barometric pressure relative to the hydrostatic pressure. The computed BE applying the linear regression between standardised hydrostatic and barometric pressure (Burgess et al. 2017) at Bhuapur ranges from 13 to 32% (Fig. S6 of the ESM). Frequency analysis of surfacewater levels in rivers Brahmaputra and Futikjani reveals no discernible dominant frequencies at repeated intervals (Fig. S7 of the ESM). Groundwater recharge estimated using a watertable fluctuation method that does not consider drainage (groundwater-level recessions) is 485 mm in 2009 and 540 mm in 2010, based on high-frequency observations, and 568 mm (σ = 82 mm year −1 ) and 589 mm (σ = 85 mm year −1 ) based on lower-frequency observations over period of 1987-2014 at sites TA013 and TA014 (Fig. 3a), respectively.

Water-level dynamics in Plio-Pleistocene deposits: Savar
At Savar, hourly monitoring (February 2009 to April 2011) within the Plio-Pleistocene environment reveals diurnal to seasonal fluctuation patterns in GWLs (Fig. 6). The amplitude Fig. 5 a Bhuapur high-resolution (half-hourly) field observations of normalized (i.e. subtracting the mean) groundwater levels (red) and hourly surface-water levels (black) in River Brahmaputra; b daily (i.e. sum of half-hourly records over 24 h) rainfall observations (blue) and daily rainfall (pink) from nearby BMD (Bangladesh Meteorological Department) station located at 25 km away from Bhuapur to fill in gaps in field monitoring records; c hourly fluctuations in groundwater levels, calculated by subtracting a 1-h moving average from the observational records. The inset on the top panel (d) is shown in Fig. S4 of the ESM of seasonal groundwater levels at Savar during the monsoon decreases in 2010 (0.8 m), relative to 2009 (1.5 m); this decline occurs alongside an overall fall in groundwater levels of 2 (~1 m year -1 ). Of note is that GWLs at Savar do not respond to extreme rainfalls. Standardised GWLs in Fig. 6e (and Fig. S8 of the ESM) show no response in the piezometric level to an extreme rainfall event with a peak intensity of 90 mm h −1 (day total of~320 mm) that occurred on the 28th of July in 2009. On the same day in Dhaka, a rainfall event of 333 mm was recorded and is the second highest oneday rainfall in the past 60 years period (Ahammed et al. 2014). Fourier frequency analysis of groundwater levels reveals responses to primarily three atmospheric signals (Fig. S9 of the ESM): S 1 , S 2, and S 3 ; S 4 is detectable but weak. In contrast to Bhuapur, the amplitudes of the solar signals in hydrostatic pressure at Savar vary from wet to dry season and curiously exceed barometric pressure during the dry-season. The computed BE at Savar (Fig. S10 of the ESM) is substantially higher than Bhuapur and ranges from 64 to 81%.
Surface-water levels in the nearest rivers (River Dhaleshwari, BWDB station number SW69; River Bangshi, BWDB station number SW14.5) shown in Fig. 2bb are influenced by tides, particularly during the dry-season low-flow condition (Fig. S11 of the ESM). Specifically, 14.76 daily cycles correspond to a lunisolar synodic fortnightly (Msf) signal that is associated with spring and neap tides for both rivers. Standardised GWLs (i.e. zero mean head oscillations) in Fig. 6c do not show any significant variations over the seasons or years. There are few episodic spikes in piezometric levels but their origin remains unclear. Groundwater recharge estimated using a water-table fluctuation method (as previously discussed at Bhuapur) is 44 mm in 2009 and 34 mm in 2010, based on high-frequency observations, and 179 mm (σ = 35 mm year −1 ) and 229 mm (σ = 38 mm year −1 ) based on lower-frequency observations over the period of 1987 to 2004 at sites DH073 and DH115 (Fig. 3c), respectively.

National-scale representivity of hydrographs at Bhuapur and Savar
The representivity of GWL regimes observed at Bhuapur and Savar was evaluated across the Bengal Basin of Bangladesh using a hierarchical cluster analysis of 464 Fig. 6 a Savar high-resolution (half-hourly) field observations of normalized (i.e. subtracting the mean) groundwater levels (red) and hourly surface-water levels (black) in River Dhaleswari; b daily (i.e. sum of half-hourly records over 24 h) rainfall observations (blue) and daily rainfall (pink) from nearby BMD (Bangladesh Meteorological Department) station located at 15 km away from Savar to fill in gaps in field monitoring records; c hourly fluctuations in groundwater levels, calculated by subtracting a 1-h moving average from the observational records. The inset on the top panel (d) is shown in Fig. S10 of the ESM, e highlights an extreme rainfall event time series records of weekly groundwater levels between 1994 and 2013. The results of the finalised hierarchical clustering are displayed in Fig. S12 of the ESM in the form of a dendrogram using Canberra distance as the measure of dissimilarity and Ward.D2 method as the primary clustering algorithm with the cluster numbers fitted at k = 5 (Fig. S13 of the ESM). For each of the five clusters identified in the Bengal Basin of Bangladesh, groundwaterlevel time series are plotted against the long-term (1994-2013) mean in Fig. S14 of the ESM. Table 2 reports the average of the time-series records for each group and their corresponding statistics.
The characteristics of each cluster are outlined in the following. Cluster 1 (CL1) is characterised by strong seasonal fluctuations in GWLs that increase in their amplitude with a very slight decline in the mean (Fig. S14 of the ESM); this cluster includes monitoring well TA014 in Bhuapur and represents 34% (165 wells) of national subset of monitoring wells (Fig. S15a of the ESM). Cluster 2 (CL2) also shows strong seasonality but without a noticeable increase in their amplitude as observed in CL1; a slight decline in GWLs occurs similar to CL1. Cluster 3 (CL3) demonstrates diminishing seasonality and much greater declining trends relative to the other four clusters. This cluster includes DH115 of Savar and represents 12% (60 wells) of the subset of national monitoring boreholes (Fig. S15a of the ESM). Cluster 4 (CL4) shows a somewhat suppressed seasonality, though much less than CL3, with moderate declining trends in both wet and dry seasons. Cluster 5 (CL5) is similar to CL1 showing slightly declining trends but with a slight dampening of seasonality relative to CL1. Notably, Cluster 5 exhibits high annual fluctuations throughout the time series compared to other clusters.

Groundwater recharge in Holocene sedimentary environments
The annual hydrograph at Bhuapur (Fig. 5a) can be broadly dissected into three distinctive phases: (1) during May to September, GWL rise indicates recharge in excess of discharge during the monsoon; (2) during September to early November, GWL declines as groundwater discharge, which occurs naturally to surface-water bodies (i.e. streams, canals, wetlands), via evapotranspiration by phreatophytic plants and as a result of sporadic pumping for supplemental irrigation around the Aman harvesting period (late November-late December), exceeds recharge following the monsoon; and (3) during late November to early May, a sharp change in the recessionary trend in the hydrographs following the onset of groundwater abstraction for dry-season Boro rice cultivation. During this third stage at Bhuapur, the amplitude of diurnal hydraulic head oscillations increases at the onset of irrigation pumping (on/off) cycles for Boro rice cultivation from mid-January to mid-May (Fig. S4 of the ESM) and has similarly been observed by Harvey et al. (2006). During the peak of the irrigation season (March/April) in 2012-2013, field surveys recorded 2638 shallow (35-55 m) pumping tubewells and 18 deep (90-150 m) tubewells operating in Bhuapur Upazila.
High-frequency GWL observations at Bhuapur show responses to observed heavy rainfall events (Fig. S3 of the ESM) that are consistent with the occurrence of direct recharge; the possibility that these piezometric responses stem from poroelastic deformation is challenged by the very shallow depth of groundwater (mean depth to groundwater was 6.25 m from 2009 to 2011) and absence of compressible sediments within the shallow surface geology comprising light brown sand (Fig. S1a of the ESM). Fourier analysis of highfrequency GWL observations at Bhuapur indicates that the amplitude of the solar atmospheric signal S2 exceeds that of other solar signals (S1, S3, S4) as predicted theoretically (Acworth et al. 2015). The muted response in hydrostatic pressure at Bhuapur (Fig. S5 of the ESM), relative to Savar (Fig. S9 of the ESM), is expected from the low computed BE (13-32%), consistent with unconfined aquifer conditions suggested from lithological logs (Fig. S1 of the ESM). Time lags are evident between monsoonal peaks in SWLs and GWLs (Fig. 5). The Pearson correlation coefficient between hourly GWL and SWL (BIWTA-2730) at Bhuapur from 2009 to 2011 improves from 0.65 to 0.86 with the incorporation of a time lag of 46 days. Using the long-term  weekly time series data of GWLs (TA014) and SWLs (BIWTA-2730), the correlation improves from 0.87 to 0.91 with a time lag of 2 weeks. Further, the hydraulic gradient is from the river to the groundwater during the monsoon (Fig. S16 of the ESM) but then reverses during the dry season. Overall, the monitoring evidence indicates a strong coupling of GWLs and SWLs in which groundwater recharge can occur not only directly (i.e. diffuse recharge) but also indirectly via leakage from adjacent surface waters. Closely corresponding but lagged individual peaks in SWLs and GWLs during the monsoon (June to October) have been described as "lateral pressure pulses" elsewhere in the Bengal Basin and attributed to recharge (Stute et al. 2007) yet may alternatively comprise poroelastic responses as suggested by recent modelling (Woodman et al. 2019). The conclusion of strong coupling between SWLs and GWLs is consistent with Ravenscroft Fig. 7 Conceptual diagrams summarizing shallow groundwater recharge processes in the a Holocene and b Plio-Pleistocene environments in the Bengal Basin of Bangladesh; blue and red horizontal lines indicate general position of groundwater levels during the monsoon and dry seasons, respectively. Note that the lithological logs from Bhuapur and Savar provide a legend for the employed shading in the adjacent block diagrams and Rahman 2003) who assert that stream-bed sediments of the River Brahmaputra are hydraulically connected to adjacent shallow aquifers.
The GWL regime observed at Bhuapur and represented by monitoring well TA014 (Fig. 5a), occurs within Holocene sediments (Fig. S15a of the ESM) and is characterised by CL1. The increasing amplitude of seasonal fluctuations that defines this cluster is consistent with the observation of increasing seasonal groundwater abstraction replenished by recharge, be it diffuse and/or focused (Fig. 7a). Such dynamics have similarly been noted in Bangladesh by Shamsudduha et al. (2009Shamsudduha et al. ( , 2011 and illustrate the Ganges Water Machine originally proposed by Revelle and Lakshminarayana (1975); the steady increase in seasonality for CL1 does not reflect mass loading observed in coastal areas (Burgess et al. 2017). Manual inspection of 55 hydrographs belonging to CL1 (33% of total) prior to the 1994 indicates amplification of seasonal oscillations associated with pumping which occurred as early as the late 1980s.
Monitoring-well hydrographs in CL1 are largely absent from the coastal zone and show a bias to areas where upper silty-clay thickness is low (Fig. S15 of the ESM) facilitating direct recharge. No clear association is evident between CL1 sites and the magnitude of groundwater abstraction for irrigated agriculture (Fig. S15b of the ESM). CL1 sites also show no clear bias to annual flood inundation depths (Fig. S15d of the ESM). In northwestern Bangladesh where annual flood inundation depths are low (<0.9 m), groundwater replenishment is expected to derive predominantly from direct, diffuse recharge. In south-central and northeastern Bangladesh as well as areas proximate to river channels where annual flood inundation depths are higher (>0.9 m), focused groundwater recharge pathways are expected to become more prominent (Fig. 7a). Focused groundwater recharge ultimately depends, however, upon favourable hydraulic gradients between flood depths and underlying groundwater levels. This evidence from statistical clustering, combined with site observations at Bhuapur, challenge assumptions (e.g., Water Resources Planning Organisation 2000a) that groundwater recharge to shallow aquifers in Holocene sediments occurs solely via diffuse recharge.

Groundwater recharge in Plio-Pleistocene sedimentary environments
In the Plio-Pleistocene environment of Savar, more gradual and smoother seasonal rises in GWLs (0.5-1 m) are observed relative to Bhuapur. GWLs lag the seasonal rise in adjacent SWLs by one and a half months (Fig. 6b). Over two consecutive years in Fig. 6a, the second peak in hydraulic head during 2010 is lower (by~1.5 m) and delayed by 1 month (from October to November) relative to 2009. A consistent 1-month lag between peak levels of shallow GWLs and SWLs is commonly observed in the monitoring boreholes in Madhupur and Barind Tracts where direct rain-fed recharge to underlying Plio-Pleistocene aquifers is constrained by a low-permeable surface geology (Fig. S8 of the ESM; Shamsudduha et al. 2011).
Confined aquifer conditions suggested by the high BE (64-81%) correspond to lithological data (Fig. S1 of the ESM) showing the presence of a low-permeability confining layer (Madhupur Clay) to a depth of between 13 and 14 m below ground level. In practice, the degree of confinement is expected to change slowly with depth in association with a fining upwards sequence from sands to clay that is characteristic of Plio-Pleistocene depositional environments in the Bengal Basin. Seasonality in the GWL regime at Savar reduces dramatically in the deeper borehole (DH115) after 2003-2004 suggesting a gradual change in hydraulic conditions that likely comprises a transition from a confined to unconfined conditions as hydraulic head declines. In Savar, sharp increases in groundwater abstraction for industry began approximately in 2000 and have impacted wells locally including DH115. This transition is commonly observed in wells across the Plio-Pleistocene terraces including the Madhupur and Barind Tracts where GWLs drop below the upper silt and clay layer due to intensive abstraction (Water Resources Planning Organisation 2000b; Shamsudduha et al. 2011;Nowreen 2017).
Seasonal rises in GWLs are consequently attributed to focused recharge occurring via incised drainage channels that have cut through the overlying surface clays into the underlying sand aquifer. Lateral discontinuity in surface clays may also exist and enable hydraulic connections between river channels and shallow aquifers. Such indirect recharge pathways to the Plio-Pleistocene aquifer have previously been proposed by Zahid et al. (2008) and Burgess et al. (2011). The reduction in seasonality that characterises hydrographs in CL3, combined with noted lags between GWLs and SWLs, do not suggest that seasonal fluctuations result from poroelastic responses to surface loading and unloading.
The observed GWL regime at Savar, as represented by DH115, is characterised as CL3 and associated predominantly with Pleistocene terrace areas of Barind and Madhupur Tracts (Fig. S15a,c of the ESM) where confining and semiconfining aquifer conditions are common (Ravenscroft and Rahman 2003). This cluster is also associated with moderate to high abstraction (as mentioned in Table S2 of the ESM) and bias to areas where annual flood inundation is low, primarily <0.3 m (Fig.  S15b,d and Table S3 of the ESM). Consequently, opportunities for localised, focused recharge from seasonal floodwaters are reduced not only by lower surface-water levels but more importantly by much less permeable surface soils. Indeed, in these areas with a substantial thickness of silty-clay cover, diffuse recharge is often insufficient to sustain high rates of groundwater abstraction for dry-season irrigation, leading to groundwaterlevel declines and substantial reductions in seasonality associated with shifts from confined to unconfined storage conditions.

Concluding discussion
High-frequency and long-term groundwater-level observations, combined with records of rainfall and river stage, reveal contrasting hydrogeological responses to the Asian monsoon in the two primary depositional environments in the Bengal Mega-delta of Bangladesh. These contrasting hydrological responses reflect different recharge processes. In Holocene sediments, rain-fed recharge is observed to occur directly and rapidly via a transmissive surface geology; indirect recharge via adjacent river channels is also strongly indicated. Despite substantial increases in abstraction, declining trends in groundwater levels are largely absent and point to increased recharge capture (i.e. induced recharge) that is consistent with the concept of the Ganges Water Machine proposed by Revelle and Lakshminarayana (1975) and observed previously in Bangladesh by Shamsudduha et al. (2011). In Plio-Pleistocene sediments, direct groundwater recharge is inhibited, even under extreme rainfall, by a relatively impermeable clay layer (Madhupur Clay) that produces confining conditions. There is, however, evidence that intensive abstraction in the noted absence of direct recharge has lowered groundwater levels below the confining layer leading to a transition from confined to unconfined aquifer conditions. Seasonal recharge in Plio-Pleistocene sediments can occur via leakage from hydraulically connected distributaries during the monsoon. Steady changes in the seasonality of groundwater-level time series (i.e. amplified seasonality for the Bhurapur cluster 1, diminished seasonality for Savar cluster 3), combined with time-lags between observed oscillations in surface water and groundwater levels at Bhuapur and Savar, suggest that seasonal groundwater-level rises do not derive substantially from poroelastic effects, noted recently in coastal Bangladesh by Burgess et al. (2017). Poroelastic effects associated with seasonal mass loading may, however, also contribute to groundwater-level responses in Holocene sediments observed at Bhuapur.
The evidence of contrasting recharge regimes in the two primary depositional environments of the Bengal Mega-delta has important implications for water supply. Groundwater depletion is pronounced in Plio-Pleistocene deposits where the low permeability of the surficial clays inhibits rain-fed recharge. The more permeable surface geology in Holocene deposits enables direct and indirect recharge, and facilitates recharge induced by intensive groundwater abstraction during the dry-season for irrigation. The current analysis is unable, however, to quantify the relative proportions of direct and indirect recharge. Nevertheless, the observation of both direct and indirect recharge processes in the Bengal Mega-delta highlights the potential limitations of recharge assessments which employ tools and models that exclude indirect recharge. As the hydrological and depositional environments of other Asian mega-deltas mirror those of the Bengal Basin, assessments of groundwater recharge and the sustainability of groundwater use will need to recognise shallow groundwater and surface water as coupled systems in which changes in one system (e.g., river discharge, groundwater abstraction) influence the other.