Buried deep freshwater reserves beneath salinity-stressed coastal Bangladesh

Domain: earth_science
Published: 2025-11-28
Editorial outcome: major_revision

Abstract
Aquifer overexploitation and saltwater intrusion threaten freshwater resources in coastal regions worldwide. In Bangladesh, arsenic contamination further reduces shallow freshwater availability. In the coastal zone, much of the shallow groundwater is saline, while the availability of deeper fresh groundwater remains poorly understood. Here, we utilize deep-sensing magnetotelluric soundings to image contrasts in electrical resistivity between fresh and saline groundwater along the Pusur River in the Ganges-Brahmaputra Delta. Our data reveal two distinct deep freshwater bodies separated by a high-salinity zone. We propose that these aquifers formed during the Last Glacial Maximum sea-level lowstand and are protected by overlying fine-grained sediments, whereas the Ganges paleovalley incision, followed by marine transgression and deposition, created the saline gap. Our work maps potential resources for this water-stressed region and suggests that the interplay between past sea-level cycles, sedimentation, and hydrogeological processes demonstrated here may also control the distribution of fresh groundwater in other deltas.

Introduction
Freshwater scarcity, exacerbated by climate change, is a severe problem in coastal regions worldwide1–5. In coastal Bangladesh, particularly in the southwestern (SW) region, high population density and saline groundwater place substantial pressure on freshwater resources6,7. Relative sea-level rise, enhanced by embankments8 and shrimpfarming9, has intensified saline intrusion and further strained freshwater availability10,11. Farther inland, shallow aquifers are often contaminated with arsenic (As) and other pollutants12–14, further limiting potable and irrigation water. It is crucial to identify sustainable management strategies and new freshwater sources to address this pressing water crisis. Deep groundwater (>150 m) is one of the vital water supplies, as it has low salinity and As levels15–17. Multiple studies have suggested the presence of deep groundwater at a regional scale6,15 and evaluated the susceptibility of deep aquifers to As contamination18,19, but little is known about their extent and the underlying hydrogeological processes.

Electromagnetic (EM) geophysical methods, sensitive to pore fluid salinity, effectively map resistivity variations in groundwater systems20–22. Freshwater-saturated sediments exhibit significantly higher resistivity than their saline counterparts, enabling EM methods to delineate freshwater resources within conductive host environments, such as saltwater-dominated zones in sedimentary basins. While near-surface EM methods have been used to map shallow groundwater and Holocene lithology in other parts of Bangladesh23,24, these techniques are insufficient to delineate deep targets. To image the deep, large-scale aquifers in coastal Bangladesh, we employed broadband magnetotelluric (MT) soundings, a non-invasive geophysical method capable of sensing electrical resistivity from depths of meters to several kilometers and beyond (see Methods). We collected MT data at 25 sites along riverbanks and fields adjacent to the Pusur River, a Ganges River distributary in SW Bangladesh, creating a 120 km-long MT transect from Khulna City to the coastline at the southern edge of the Sundarbans mangrove forest (Fig. 1c). With an average site spacing of about 5 km, the transect captures the large-scale spatial distribution of fresh and saline groundwater in the study area.

The Ganges-Brahmaputra Delta (GBD), encompassing over 100,000 km2 and formed by the convergence of these two eponymous Himalayan rivers, is the largest delta in the world (Fig. 1). Two-thirds of the GBD lies within Bangladesh, bordered by the Indian Shield to the west, the Shillong Massif to the northeast, and the Indo-Burma Ranges to the east25. Annually, the Ganges and Brahmaputra Rivers transport an estimated 1 billion tons of sediment eroded from the Himalayas to the delta26. The sediment load, ranging from coarse sand to clayey silt27,28, consists of medium to coarse sand (28%), very fine to fine sand (38%), and muds (34%)—proportions that are roughly equal within the delta29. These sediments are coarsest upstream and progressively fine downstream30. Pleistocene terrace uplands, such as the Barind Tract and Madhupur Terrace (Fig. 1b), act as natural barriers controlling sediment dispersal and channel migration28,31,32. Meanwhile, the Sundarbans mangrove forest, located in SW GBD, functions as a sediment storage area and a fresh-to-saltwater transition zone for surface water. High sedimentation rates33,34 drive basin subsidence due to compaction and flexural isostatic loading, increasing toward the coastline, and simultaneously contribute to delta-front seaward progradation35–38. This substantial sediment budget is primarily transported during the wet summer monsoon season, when high precipitation and river discharge facilitate sediment delivery to the lower delta plain39.

Sea-level fluctuations, coupled with sedimentary processes, have been dominant factors controlling the formation of deep aquifers in coastal GBD40,41. Marine transgressions and regressions significantly influenced stratigraphy and depositional environments, affecting the emplacement of both fresh and saline water. During the Pleistocene sea-level lowstand, the Bengal delta shoreline regressed 150–200 km offshore42,43, exposing the continental shelf and facilitating freshwater recharge over vast areas. Conversely, the rapid sea-level rise during the early Holocene caused the tidal zone and shoreline to migrate ~100 km inland36,40 from their present-day positions. The stronger monsoon during the early Holocene increased river discharge and sediment supply, initiating seaward progradation of the GBD30,31. This dynamic interplay between sea-level cycles and sediment supply created conditions conducive to forming deep freshwater aquifers overlain by shallower saline groundwater. Delineating the freshwater distribution and dimensions of these coastal aquifers is a primary objective of this study, as their extent remains largely uncharted.

Results
Electrical resistivity of deep groundwater systems
The regularized 2D inversion of the MT data (see Methods) yielded a resistivity model characterized by two prominent resistive zones, R1 and R2, separated by a conductive zone, C1 (Fig. 2a). R1 extends down to 800 m depth along a 40 km span in the northern part of the transect, with a slight southward dip in its upper boundary. R1 is likely to extend farther north beyond the survey area. R2, located within the central Sundarbans, spans 40 km and reaches a depth of 250 m. Notably, R1 exhibits a broader zone of resistivity exceeding 10 Ωm compared to R2, with a maximum resistivity of ~20 Ωm, double that of R2. The transition between resistive and conductive zones is more gradual along the boundary of R1 than that of R2, particularly the lower boundary.

The conductive zone C1, situated in the northern Sundarbans, forms a ~20 km wide discontinuity between resistive zones R1 and R2, potentially connecting vertically with a deeper conductor, C2, extending from 700 to 1800 m depth. A shallower conductive zone, C3, confined within the upper 70 m near the coastline in the southern part of the Sundarbans, exhibits the lowest resistivity (0.55 Ωm) encountered along the entire transect. This shallow layer gradually increases in resistivity northward, away from the coastline.

Salinity interpretation
We estimated groundwater salinity by converting bulk resistivity to pore fluid salinity, using lithological information from well logs and a model of water conductivity as a function of salinity (see Methods). The calculations incorporated variables, such as temperature, pressure, sediment porosity, compaction, and cementation factor. Our preferred salinity model (Fig. 2b) defines freshwater (<3 psu) and saltwater (>30 psu) zones, with the 10 psu contour marking the transition between low- and high-salinity brackish water. This model reveals a clear correlation between resistive zones (R1, R2) and low-salinity aquifers, as well as between conductive zones (C1, C2, C3) and high-salinity groundwater.

Freshwater reserves
Although our results provide first-order details of deep fresh aquifers along the Pusur River, explaining the origin of these aquifers and the factors controlling their distribution using resistivity models alone is challenging. Therefore, we incorporated additional nearby data—including stratigraphic information, as well as groundwater age and salinity measurements from boreholes—to further investigate aquifer characteristics (Fig. 3). Lithological and isotopic data from sediments in boreholes up to 90 m depth in the northern part of our profile28,30,44 reveal spatial variations of the Holocene-Pleistocene sediment boundary (Fig. 3). Shallow Holocene aquifers are delineated in the upper 90 m7, whereas the deep aquifers (>100 m) identified in this study are likely associated with Pleistocene or older strata45.

We propose that resistive zones R1 and R2 represent primarily freshwater reservoirs formed due to the low sea-level stage in the Last Glacial Maximum (LGM) and subsequently sealed by a deposition of finer-grained, low-permeability layer44. During the LGM lowstand, when sea level was as much as 120 m lower than today, and the major rivers incised deep, wide paleochannels across the GBD32,40,44,46, converging in the lower delta (Figs. 3 & 4). At that time, the shoreline was near the present continental shelf edge41, exposing the entire shelf to freshwater recharge45 (Fig. 4a). Weathering of the Pleistocene interfluves exposed during glacial periods (Supplementary Fig. 7) created a stiff clay layer known as the Last Glacial Maximum Paleosol (LGMP)44. The subsequent rapid sea-level rise flooded the coastal GBD45,47, trapping freshwater beneath fine-grained transgressive sediments and the LGMP31,44,47. The LGMP acts as a barrier inhibiting vertical mixing between dense saltwater and lighter freshwater below44,47. However, marine sediments deposited within paleovalleys dissect the lateral continuity of the LGMP44,46,47, such as the NW-SE oriented Ganges paleovalley that intersects the MT profile in the northern Sundarbans (Figs. 3 & 4b) and forms conductor C1. Since the LGMP is crucial for protecting fresh groundwater, its lateral discontinuity across the lower delta43,47–49 strongly influences the distribution of these aquifers.

Our imaging results suggest that resistor R1 is potentially part of a larger deep freshwater system extending tens of kilometers or more to the north, as R1 maintains its full thickness northwards (Fig. 2a). The 10 Ωm contour at the top of R1 (Figs. 2 & 4c) might delineate the LGMP, which separates the thick resistor R1 from the overlying thin, conductive Holocene sediments that dip southward. The LGMP layer (subcrop in Fig. 3) inhibits local vertical recharge and prevents saltwater infiltration from shallow Holocene aquifers to R1. Therefore, if active recharge to R1 exists, it likely occurs farther north via lateral flow. The alternating layers of fine- and coarse-grained sediments reduce effective vertical hydraulic conductivity, facilitating the regional lateral flow19. Due to the minimal topographic relief of coastal Bangladesh and the low hydraulic gradients, natural horizontal recharge to deep aquifers at regional scales would require thousands of years or more18,19,45,50. Although the LGMP inhibits vertical recharge, Holocene-Pleistocene groundwater mixing may still occur near the southern edge of freshwater R1, where the LGMP is absent. In our models (Fig. 2), this mixing potentially occurs within the gradual vertical transition zone between R1 and C1 via vertical salt diffusion51. In addition to geophysical imaging and mapped stratigraphy, regional isotopic data provide insights into the origin and flow patterns of coastal aquifers in Bangladesh. A decrease in deep groundwater 14C activity from northern Bangladesh toward the coastline, along with δ18O enrichment in deep southern aquifers52, suggests slow seaward recharge through both vertical and horizontal flow pathways and long residence times of up to ~25 kyrs. Since R1 is likely the southern, distal end of a larger coastal aquifer located below the mapped interfluvial LGMP subcrop44 (Fig. 3) and lies in an area of low 14C activity52, we infer that R1 may partially to mostly contain Pleistocene-age freshwater.

Groundwater 14C ages, based on dissolved inorganic carbon in the interfluve and paleochannel areas that go across the Ganges paleovalley located to the northwest of our MT transect (Fig. 3) provide further constraints on the origin and dynamics of deep freshwater reserves13,16. Within the Ganges paleovalley area, without the LGMP layer and where the Ganges paleovalley fill is sandy (Fig. 3 transect at 23.2°N), reported uncorrected 14C ages that cluster around 2.0 kyrs at depths of 150–250 m, indicating the downward flow of Holocene As-contaminated groundwater13. In contrast, in the interfluve area associated with the LGMP farther south where the paleovalley is muddy (south of 23°N in Fig. 3), most groundwater samples have 14C ages approximately 10-12 kyrs, with two low-As wells exhibiting Pleistocene ages of 15.1 kyrs at 162 m and 29.6 kyrs at 249 m depth13 (Fig. 3). The older ages in interfluve areas reflect minimal vertical As-contaminated groundwater transport. We reported the published uncorrected 14C ages instead of corrected ages because of the high uncertainty in correcting for the initial 14C content13, however, these 14C ages delineate a clear spatial pattern in the distribution of deep groundwater ages. The older groundwater in interfluve areas suggests that deep confined systems protected by the LGMP, like the R1 reservoir, are likely to preserve paleo water at least from latest Pleistocene with limited slow regional recharge during the Holocene.

The resistive body R2, located within the Sundarbans mangrove forest, appears to be an isolated aquifer disconnected from the deep regional groundwater flow system due to the incision of the Ganges paleovalley (Figs. 2 & 3). This apparent isolation suggests that R2 are primarily relict Pleistocene groundwater with minimal connection to modern active recharge. The smaller thickness and higher salinity of R2 compared to R1 support this interpretation. An alternative explanation revealed by LGMP and isopach mapping (Fig. 3) is that R2 could be connected to flow paths originating in West Bengal India. In this scenario, R2, similar to R1, would contain paleo-groundwater emplaced prior to the Holocene, potentially connected to distal upland recharge farther west at extremely slow rates due to low topographic gradient. The LGMP surface may end or become deeper than the maximum depth of tube wells at the southern limit of data44. In the Sundarban Mangrove Forest, the LGMP cannot be mapped due to the lack of wells. If the LGMP did not develop as far south as R2, it would reduce R2’s isolation from saline Holocene groundwater and regional mixing. Additionally, the increasing proportion of clayey silt toward the coastline may lower the bulk resistivity of R2, due to enhanced surface conduction effects in clay mixtures53. We neglected the effects of clay content on formation resistivity since saltwater conductivity is the dominant factor.

High-salinity zones
The 20 km-wide conductive zone C1, representing a saline gap between freshwater aquifers R1 and R2, corresponds with the NW-SE oriented Ganges paleochannel incised through the northern Sundarbans, as indicated by tube well data52,54 (Figs. 2 & 3). This paleochannel formed during the Pleistocene lowstand, when the Ganges River carved a deep valley and transported sediments seaward (Fig. 4b). During that lowstand, both the Ganges and Brahmaputra Rivers emptied into the shelf-indenting Swatch of No Ground (SoNG) canyon (Fig. 1b). As sea levels rose in the Holocene, saltwater inundated the paleochannels, depositing marine or estuarine sediments and creating the observed saline gap. Although tube well data only extend to a depth of 91 m, the paleochannels once connected to the present-day SoNG canyon, itself incised to depths greater than 500 m during the Pleistocene42 (Figs. 1 & 3). We speculate that the paleochannel could extend up to a few hundred meters deep at our MT transect (Fig. 4b, c), providing abundant accommodation space for saltwater-filled sediments during marine transgression, thereby giving rise to feature C1.

Our data reveal a thin, relatively conductive layer overlying freshwater aquifer R1, which correlates well with the south-dipping Holocene deposits (30 to >90 m thick) documented in tube well data29. Increased sediment supply during the early Holocene, driven by a stronger Indian monsoon, resulted in the progradation of a broad tidal zone over the marine Holocene deposits27,38,55,56. Shallow salinity decreases northward from 30 psu at conductor C3 near the coast to 3 psu at the northern end of the transect, 120 km inland, coinciding with the maximum Holocene transgression40. This gradient mirrors the high surface water salinity observed during the dry season, extending inland to the topographic break between the fluvial and fluvial-tidal delta50. These findings confirm extensive shallow saltwater intrusion within the coastal GBD, contributing to freshwater scarcity during the dry season.

The broad conductive region below two freshwater aquifers, between depths of 1 and 2 km (Fig. 2), could represent marine clay extending across the Bengal Basin. The offshore Sangu-1 well (Fig. 1) yielded sands down to ~600 m depth, overlying sandy shale57 within the Pleistocene Megasequence 3 (base 0.4–1.9 Ma58). Since conductor C2 exhibits resistivity values lower than its surrounding and salinity levels close to those of seawater, it may represent offshore shelf deposits seaward of the coast. Based on seismic data43, we hypothesize that this deep saline unit formed during marine deposition prior to large-amplitude sea-level cycles that extended the shoreline and deposited coarser-grained sediments near the shelf edge. These processes subsequently facilitated the emplacement of freshwater bodies R1 and R2.

Sea-level change and the formation of deep freshwater resources
Our MT survey along the Pusur River reveal two distinct freshwater bodies: a seaward-dipping wedge extending to a depth of 600 m in the northern part of the transect, and a smaller, shallower fresh-to-brackish zone located closer to the river mouth. Corroborating evidence, including Holocene sediment thickness, LGMP distribution associated with interfluves and paleovalleys, and nearby groundwater 14C ages, supports the conclusion that the changes in sea-level cycles and subsequent stratigraphic deposition have governed the formation and distribution of the deep freshwater reserves, R1 and R2. A high-salinity zone between the two freshwater bodies correlates with the Ganges paleovalley, which was incised during the LGM sea-level lowstand and subsequently infilled with conductive muddy sediments during the marine transgression. Our results highlight the critical role of paleo-incised valleys in controlling the spatial distribution of fresh and saline groundwater as they are geological features commonly formed in coastal plains worldwide during periods of low sea-level51,59,60. Furthermore, paleovalley areas also serve as potential conduits for saline water intrusion into deeper flow systems, thereby contributing to coastal aquifer salinization. This onshore groundwater system mirrors recently studied offshore systems61,62, where interactions between sea-level fluctuations, sediment transport, groundwater flow, and geological processes shape freshwater formation and distribution over geologic time scales. Our findings suggest that the two identified aquifers represent notable freshwater reserves requiring millennia to develop. These deep aquifers could serve as critical freshwater sources for saline-affected coastal regions in Bangladesh6,15. However, their slow recharge rates necessitate careful integration into long-term water management strategies. Compared to other near-surface EM techniques used in groundwater exploration, MT soundings can image the full depth extent of deep groundwater systems beyond just the water table or saltwater-freshwater interface. Understanding the dimensions of aquifers allows for volumetric estimation, which is essential for sustainable water management. While deep freshwater has been locally exploited in coastal Bangladesh, our MT-derived salinity model, combined with lithological data, provides the framework for understanding the distribution of fresh and saline deep groundwater systems. This framework is valuable for strategically locating future deep extraction wells and determining pumping rates. Similar interactions of sea level cycles and fluvial incision may shape deep fresh aquifers in other coastal deltas and continental margins, as sea-level lowstands were widespread globally during the Pleistocene58–61.

Methods
MT data collection and processing
The MT method is a passive EM geophysical technique that measures the temporal variations of electric and magnetic fields to image subsurface resistivity structures. The primary source fields for MT data at lower frequencies (typically below 1 Hz) originate from the interaction between charged particles from the solar wind with Earth’s magnetosphere and ionosphere, generating fluctuating ionospheric electric currents, which in turn induce EM fields that propagate downward through the atmosphere and into the subsurface. At higher frequencies (typically above 1 Hz), global lightning discharges are the dominant source of EM fields for MT measurements. As the Earth’s interior is more conductive than the atmosphere, the penetrating EM fields interact with the resistivity structures and induce EM signals that are measured across six decades of frequencies.

In March 2022, we collected broadband MT data at 25 sites along riverbanks and agricultural fields adjacent to the Pusur River, a distributary of the Ganges River in SW Bangladesh. Our survey extended from Khulna City to the coastline (Fig. 1c). We conducted overnight MT deployments with roughly 5 km spacing, yielding MT responses from 10-4 to 2×103 s periods. This wide bandwidth allowed for sensitivity at multiple depth ranges. Cultural noise near Khulna City, such as 50 Hz electrical power lines, brick factories, farming, and livestock, degraded the MT signals at five sites in the northern part of the transect. However, we observed a significant increase in the signal-to-noise ratio at sites within the Sundarbans Mangrove Forest. Approaching the coastline, the MT responses at periods greater than 2 s become noisy, likely due to seawater motional induction noise generated by waves and strong water currents. We used remote referencing and robust processing techniques to estimate broadband MT impedance responses63. For most sites, we excluded bad data with large error bars and outliers and set an error floor of 5%. We discarded two northernmost sites (P25 and P26) because of their high noise level and set an error floor of 10% for P01, P02, and P03. These processes generally yielded high-quality MT responses in the period range of 10-3 to 103 s (Supplementary Fig. 1).

Electrical resistivity inversion methods
To obtain electrical resistivity models, we applied two complementary inversion techniques to the MT data: conventional regularized inversion and a Bayesian sampling method. The non-linear regularized Occam’s inversion approach64 generates a smooth model that fits the data within a specified tolerance, using a model roughness penalty to minimize the introduction of unnecessary structures. This approach is computationally efficient for identifying prominent resistivity features but does not provide quantitative estimates of model parameter uncertainty. Recognizing that model uncertainty is inherent in geophysical inverse problems due to data noise, model simplification, and computational limitations65, we also employed a Bayesian sampling method to quantify the range of plausible models that fit the observed data.

For the regularized inversion, we employed MARE2DEM66, a non-linear regularized Occam’s inversion code, to invert log10 apparent resistivity and phase of MT responses into a 2D bulk resistivity model. The inversion domain was parameterized using a dense quadrilateral grid in the upper 2 km to resolve near-surface groundwater systems, transitioning to a sparser grid from 2 to 10 km depth to accommodate potential deep structures influencing the longer-period data. A horizontal-to-vertical roughness penalty ratio of 200 was assigned to reflect both the large length-to-depth ratio (120 km: 10 km) of the survey area and to favor models with laterally continuous features, given the approximately 5 km station spacing, which is much wider than our 0 to 2 km depth of interest for groundwater mapping. To account for the conductive effect of the Bay of Bengal at the southern end of the transect, we incorporated an adjacent ocean region with a fixed resistivity of 0.3 Ωm. Our preferred 2D resistivity model (Fig. 2a) achieved a root-mean-square (RMS) misfit of 1.005 after 13 iterations.

In the Bayesian sampling approach, the log10 resistivity model parameters are treated as random variables. We compute their probability density function (PDF) as a function of depth using a parallel tempering reversible jump Markov chain Monte Carlo algorithm65,67. By applying parallel tempering, the MCMC algorithm can more effectively explore the entire model space and avoid getting trapped in the local maxima of the probability distribution, thereby reducing potential biases in the sampling results. Additionally, we treat the model parameters, such as the number of resistivity layers, their thicknesses, and depths, as unknowns and include them in the inversion process. The algorithm applies the concept of Bayes’ rule:

p(m|d) / p(d|m) × p(m) (1)

where m is a vector of model parameters describing subsurface resistivity structure, and d is the observed data vector. The p(d|m) term is the model likelihood, which measures the probability of observing the data given a specific model. The p(m) term is the assumed prior probability distribution of the model parameters and is independent of the observed data, reflecting prior knowledge or assumptions about the subsurface. On the left-hand side, our target is the posterior model probability distribution p(m|d), which quantifies our updated knowledge about the model parameters after considering the observed data. The accuracy of the posterior distribution estimation relies on the careful selection of the likelihood function and prior distribution.

Supplementary Fig. 2 shows representative Bayesian 1D inversion models at four sites, overlaid with profiles extracted from the 2D regularized inversion results. The agreement between the regularized inversion profile and the high-probability density (PDF) regions in the Bayesian models at each site confirm the interpreted distributions of resistive freshwater and conductive saltwater zones at depth. Furthermore, the region between the 5th and 95th percentiles of the posterior PDF represents the range of resistivity values compatible with the data at a given depth, as 90% of the models in the Bayesian ensemble fall within this range.

Although both Bayesian and regularized inversions resolve three resistivity layers in the upper 2 km, the Bayesian models reveal that conductors are better constrained than resistors, highlighting the inherent higher uncertainty associated with resistors in MT inversion. This observation is consistent with the MT method’s sensitivity to the conductivity-thickness product; MT is inherently more sensitive to conductors than resistors and thick resistors than thin ones. Given the significant number of Bayesian models indicating the resistivity values for the fresh aquifer vary from 10 to 500 Ωm, we postulate that the 10 Ωm estimate provides a reliable lower bound for the freshwater layer. Additionally, the Bayesian method offers improved resolution of the transition depth between resistive and conductive zones. While the regularized models show a smooth change, the high probability regions in Bayesian models clearly delineate the possibility of a sharp boundary between these two anomalies (Supplementary Fig. 2).

Salinity estimation
From our resistivity model, we calculated the pore fluid resistivity ρf using the empirical Archie’s law68:

ρf = ρb × ϕm (2)

where ρb is the bulk resistivity obtained from our 2D inversion results, ϕ is the porosity, and m is the cementation exponent. We estimated porosity as a function of depth ϕ(z) using exponentially dependent Athy’s law69:

ϕ(z) = ϕo × e-cz + ϕmin (3)

where ϕo + ϕmin is the initial surface porosity, ϕmin is the minimum porosity at depth, and c is the compaction factor. Based on constraints from seismic data (see Supplementary Fig. 6), we set ϕo = 40%, the minimum porosity ϕmin = 15%, corresponding to a surface porosity of 55%, and the compaction factor to 0.4 km-1 for the upper 2 km. The cementation exponent is poorly constrained and can vary significantly with lithology. Therefore, we conducted a sensitivity analysis using a range of cementation values from 1.0 to 2.5, aiming to achieve an estimated pore fluid salinity near the end of the profile that is at or less than the seawater value of 35 psu.

We converted pore fluid resistivity to salinity using the Practical Salinity Scale 1978 (PSS-78)70. The PSS-78 is a dimensionless scale that defines salinity based on the electrical conductivity of seawater, incorporating the effects of temperature and pressure. This approach is crucial for accurately interpreting geophysical data in marine and coastal environments. Considering surface water temperature varies between 25 °C and 28 °C in tropical regions, we incorporated a geothermal gradient of 28.5 °C/km71 in the conversion to account for increased groundwater temperature with depth (Supplementary Fig. 5). Since temperature affects fluid conductivity, and porosity governs the volume of groundwater, these two factors are crucial for salinity estimation (Supplementary Figs. 3 and 4). In our study area, the clay content is significantly lower than the proportion of sand and silt sediments, particularly for the deeper Holocene. Moreover, since the conductivity of saltwater is significantly higher than the conductivity of clay minerals; therefore, we can safely neglect clay conduction effects when applying Archie’s law to estimate fluid resistivity.

While the regularized 2D inversions used for salinity estimation are inherently non-unique, our Bayesian models provided valuable insights into the uncertainty associated with these salinity models. It is well-known that MT data are unable to resolve the magnitude of highly thin resistive layers due to response saturation. For example, the Bayesian 95th percentiles run up close to the prior upper bound of 1000 Ωm, whereas the 5th percentiles are more clearly resolved by MT data (Supplementary Fig. 2) as they are well above the prior low bound. Therefore, the 10 Ωm value, corresponding to the 3 psu contour in the salinity model, represent the lower bound of the resistivity PDF for the freshwater zones in the Bayesian models (Supplementary Fig. 2). Given the inverse relationship between resistivity and salinity, this suggests that the actual salinity of the fresh aquifers could be even lower than 1 psu—the threshold commonly used to define potable freshwater in hydrology. Therefore, we used a combination of 10 Ωm and 3 psu contours as the bounded values to delineate the freshwater zones. Although further validation with deep boreholes and groundwater samples is needed, our estimated salinity model provides a robust initial picture of onshore deep fresh and saline groundwater distribution in coastal Bangladesh. This highlights the role of glacial sea-level fluctuations in emplacing and preserving deep fresh aquifers in coastal zones, opening avenues for further investigation.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript by Le et al. presents the results of an interesting study to delineate deep, freshwater zones in the subsurface in Bangladesh. There is a desperate need for new sources of freshwater for the population as the easier to access to groundwater is from sources contaminated with arsenic and/or affected by increasing salinity. The use of magnetotelluric soundings is a good for identifying broad areas of resistivity contrasts in the subsurface to depths much greater than near-surface methods. Therefore, it is an excellent method for identifying larger bodies of freshwater with relatively minimal effort of data collection over a relatively large area.\nThis investigation was successful in identifying 2 zones of freshwater extending to depths of 600 m surrounded by regions of saline water. The origin and location of these two zones were placed in the geological context of the Ganges Delta region, relating to the influence of sea-level changes during the glacial and post-glacial periods. Therefore, the freshwater zones are interpreted as being Pleistocene in age. The success was due in part to the processing of the MT data and using methods that could significantly reduce error to allow more confidence in processing results, which therefore strengthens the interpretations.\nThe success of this study has several implications. First, they were able to demonstrate that the use of MT can identify large areas of freshwater in the subsurface in an area that has desperate needs. The two identified zones have potential as sources of freshwater for decades or more, pending proper water management. Secondly, they demonstrated that MT could be an important method for similar studies in other coastal aquifers throughout the world,not only identifying locations of freshwater but also for improving interpretations of the geologic history of an area as the presence of freshwater zones needs to be accounted for by geologic processes.\nI find the extended data/tables/graphs very helpful.\nI do have some concerns. The age of the northern deeper freshwater zone are taken from 14C data collected from a depth of 90 m, while these zones extend to depth of nearly 600 m. Also, uncorrected ages were reported which can sometimes lead to significant differences once corrected. However, the 30,000 yr age could put a lower limit to the age. Though not available, an estimate of age of the seaward zone would be beneficial. I would like to see a bit more on the methodology of data collection and the instrument that was used to collect data. This would allow the reader to evaluate if the methodology/equipment are appropriate for this purpose.\nHere are some specific comments:\nIn line 28, embankments are mentioned. It might be nice to be a bit more specific as to what type, even if it is in the reference.\nAround line 43, it would be helpful to refer to figure 1.\nIn a few places, word deep is used. This is a relative term and it is used to describe both the 90 m depth upper aquifer and the water zones that extend to 600 m depth. R1 and R2 are referred to as deep and so is the depth of 90 m that the 14C ages were determined.\nIn Figure 1b, the location of the seismic lines are very difficult to see, as well as the legend. I am not sure this figure will be larger if published.\nFigures 2 are also difficult to read, especially 2c. Personally, I would have liked to see in the caption that the distances are relative to the location of P11.\nI am not sure where Table 1 applies.", "reviewer_2": "The manuscript addresses deep groundwater resources in the Ganges-Brahmaputra Delta, Bangladesh, applying magnotelluric (MT) soundings to delineate distribution and dimensions of deep fresh and saline groundwater zones in this coastal system. The authors recorded data at 25 sites along a 120 km transect adjacent to the Pursur River in southwestern Bangladesh with a vertical extension in a km-scale and identified two distinct deep freshwater bodies, which they attribute to being remnants of preserved Pleistocence groundwater. The authors base their findings on deposition environments /sedimentary processes and pleistocence sea-level fluctuations forming the deep groundwater systems and recorded electrical resistivities of groundwater.\nThe manuscript appears technically sound and is very well-written. It applies established up-to-date methods. The topic is of scientific interest and well-worth publishing. However, the authors should make the global relevance of their shown findings more clear and identify the their study’s unique characteristics in order to distinguish it better from pure regional applicability. In addition, the authors should highlight the potential of this study relating to implications on other (global/overregional) findings and set this into a broader context.\nL27: add more recent references\nL29-31: provide references\nL46/47: stress significance (->large scale)\nL 48-62: This section is rather detailed compared to the subsequent section and may be condensed, (and some hydrogeological aspects are suggested to be included)\nL 64: Define deep in your context\nL72/73: As indicated in the main comments, this is one point that could be used for highlighting the relevance of the presented study\nL93: provide short information why it is the preferred salinity model\nL93: define why 3 psu is chosen (and not 1)\nL98: specify which isotopes or at least state if sediment or water is meant here\nL98/99: In Fig. 3 150-250m are given (instead of: 90m), or is another figure meant here?\nL101-103: connection to onshore aquifers?\nL12: Fig. 4 a and b cover depths up to 150 m, R1 and R2 fresh water extend (much) deeper and only the topmost parts overlap Fig 4a and b, this should be mentioned more clearly when deriving analogies and processes\nL112: why are the deductions made from uncorrected Radiocarbon data? (especially as this (published) data is used as a strong argument here for fossil water in the investigated system)\nL106-118 preserved Pleistocence groundwater vs. resemblance to groundwater system offshore New Jersey: It should be clarified, if only the morphology or the (emplacement) processes are expected to be similar\nL127: refer to Fig. 2\nL152: Sea-Level Change\nL302: refer to Fig. 1c\nL456: add link\nFigures and Tables:\nAre Figures 4a and 4c and Table 1 mentioned in the manuscript text?\nL 470: please delete this line\nFig. 1: legends in a and b are too small, an in b the red entries are difficult to see\nL475: not clear how subaqueous delta depth an SoNG can be distinguished\nL481 not clear where shrimp farms are (dark blue to purple appears a bit masked)\nFig. 2: The axis gives distances (+-60 km), where P11 appears to be the center. Why did you choose this scale (and the center at P11) instead of 0 to 120 km? And why are Sundarbarns explicitly mentioned here?\nFig. 3: connection between Fig. 3 and Fig. 4 somewhat difficult\nFig. 4: maybe insert R11 and R2 into Fig. 4c", "reviewer_3": "Thank you for the opportunity to review your manuscript titled \"Deep fossil freshwater in coastal Bangladesh\". The authors have a fascinating dataset with their N-S MT transect across a portion of the southern Ganges delta. Admittedly, the geophysics underlying it are outside my field of expertise and I accepted the profiles derived from the MT data as correct during my review. The authors demonstrate two zones of low resistance groundwater which is justifiably interpreted as freshwater at depths up to 600m-800m (R1) and up to 200m-300m (R2). Utilizing data from other studies in the Ganges delta and further afield, they construct a conceptual model to explain their findings. Unfortunately, the conceptual model is plausible but not well-supported with the data as presented. In particular:\n(1) The carbon-14 data used to support the conceptual model lacks geographic specificity, which leads to apparent overinterpretation of the data. The data comes from one study (ref #5, Kahn et al. 2019). Most of the carbon-14 ages in Kahn are Holocene and the subset that are Pleistocene are not identified and located relative to the current study area. Some of the points in Kahn are more than 100km from the closest boundary of the study area. Extrapolating over such a distance is unreasonable in most contexts and certainly in a major delta complex such as the Ganges. Further, the authors argue for the water in R1 and R2 to have been emplaced during the last glacial maximum -- only 7 of the carbon-14 dates in Kahn fall in the range of the LGM, which makes understanding the location of Kahn's data relative to the current study even more crucial. Confidence in Pleistocene ages for the vast majority of Kahn's samples nearest to the subject study area is key. If they're Holocene, the entire conceptual model falls apart. To invoke Pleistocene emplacement of R1 and/or R2, the authors must unequivocally demonstrate the presence of Pleistocene water in close proximity to their transect.\n(2) In the development of their conceptual model, the authors argue that the fresh water identified in R1 and R2 is relict, connate water emplaced during the LGM, as opposed to the terminal portion of an active groundwater system. At several times they contradict this conceptualization by referencing upgradient connectivity and making comparisons to the Ganges system to an offshore system in New Jersey (specifics noted in comments in PDF file). The contradictory statements in the text undermine their conceptual model. Further, the northern extent of the R1 freshwater region extends an unknown distance outside the study area. There is no evidence presented for the reader not to envision this lobe of freshwater continuing upgradient as part of a continuous, active groundwater flow system.\nThe deficiencies in this study principally stem from its heavy reliance on data from existing published work in areas that do not exactly overlap with theirs to support their interpretations of their interesting MT transect results. Without site specific data to support their interpretations, the authors often find themselves at odds with the data and with their proposed conceptual model. The original work (the MT transect) appears sound and compelling; its interpretation much less so."}, "author_response": "Point-by-point response to referees’ comments\nDate: February 20, 2025\nManuscript Number: NCOMMS-24-67705\nTitle of Article: Deep fossil freshwater in coastal Bangladesh\nName of the Corresponding Author: Huy Le\nEmail Address of the Corresponding Author: hdl2115@columbia.edu\nDear Editor,\nWe thank all three reviewers for their valuable time and insightful feedback that has helped improve the\nmanuscript. We have incorporated their suggestions into the main text of our revised manuscript. Below\nwe provide the reviewer comments in black text and our replies in blue. We have also included responses\nto review feedback that was given in the manuscript PDF.\nBest regards,\nHuy Le and coauthors\nSpecific Responses:\nResponse to Reviewer #1:\nComment 1: In line 28, embankments are mentioned. It might be nice to be a bit more specific as to what\ntype, even if it is in the reference.\nReply 1: The embankments mentioned in the text are river embankments. We did not extend the context\nof the embankments further since the main idea of this section is introducing water problems in\nBangladesh.\nComment 2: Around line 43, it would be helpful to refer to figure 1.\nReply 2: We added “(Fig. 1)” in line 45 to clarify the survey profile with figure 1.\nComment 3: In a few places, word deep is used. This is a relative term and it is used to describe both the\n90 m depth upper aquifer and the water zones that extend to 600 m depth. R1 and R2 are referred to as\ndeep and so is the depth of 90 m that the 14C ages were determined.\nReply 3: The word “deep” is used in the main text to separate the MT-imaged aquifers from shallow\nHolocene aquifer. The deep Pleistocene fresh groundwater in our study area is primarily protected by the\nLast Glacial Maximum Paleosol (LGMP), whose spatial distribution and thickness vary across the delta\n(Fig. 3). Ravenscroft et al. (2018) Sarker et al. (2021) and other groundwater studies in Bangladesh\ndefined deep groundwater systems differently. Previously, we used the word “deep” and did not specify\nthe fixed value to consider the complexity of variations of groundwater systems. In the revised edition,\nWe added “(>100 m)” on line 106 to clarify the depth extension of deep aquifers. Note that the upper part\nof R2 is shallower than 100 m and might not fit to the defined depth range.\nComment 4: In Figure 1b, the location of the seismic lines are very difficult to see, as well as the legend,\nI am not sure this figure will be larger if published.\nPoint-by-point response to referees’ comments\nReply 4: The resolution in the word document is not great. We hope the resolution in the published\nversion is better to see the seismic lines. The legend shows elevation in meters in color scale.\nWe enlarged the legends for Fig. 1a and b.\nComment 5: Figures 2 are also difficult to read, especially 2c. Personally, I would have liked to see in\nthe caption that the distances are relative to the location of P11.\nReply 5: Again, we hope the resolution gets better in the published version. Another way to improve the\nreadability of figure 2 is rotating the figure 90 degrees and enlarging Fig. 2c. Fig. 2c is important because\nit links bulk resistivity in Fig. 2a with salinity in Fig. 2b via temperature, cementation factor, and porosity.\nThe zero marker of the distance is chosen to be the center of the profile at the beginning of the inversion.\nIt will add unnecessary workload to rescale the distances relative to the location of P11.\nWe moved the legend of Fig. 2c to the space below it.\nComment 6: I am not sure where Table 1 applies\nReply 5: Table 1 went with Fig. 2 because it lists the bulk resistivity and estimated salinity of features in\nFig. 2a and 2b. Although, the color scales in Fig. 2a and 2b indicate clearly the values. Since the table is\nredundant, we removed it from the revised manuscript.\nResponse to Reviewer #2:\nComment 1: L27: add more recent references\nReply 1: We added He, C. et al. (2021). We keep the number of citations close to 50 references to match\nNature Communications suggested limit.\nComment 2: L29-31: provide references\nReply 2: We added van Geen et al. (2008) and Rahman et al. (2020) papers.\nComment 3: L46/47: (-> large scale)\nReply 3: L45-47: “With an average site spacing of about 5 km, the profile captures the large-scale spatial\ndistribution of fresh and saline groundwater in the study area.”\nComment 4: L48-62: This section is rather detailed compared to the subsequent section and may be\ncondensed, (and some hydrological aspects are suggested to be included)\nReply 4: This paragraph provides geological and lithological background relevant to groundwater system.\nLithology, as shown later in the text, controls freshwater distribution. Adding hydrological background\nmay create duplication in the discussion, such as slow horizontal recharge due to low topography and\nvertical mixing prohibited due to LGMP.\nWe added hydrology background in the first paragraph.\nComment 5: L64: Define deep in your context\nReply 5: We defined deep in Fossil fresh aquifers section on line 104-107. We added more details of\nhydrology background in the first paragraph.\nComment 6: L72/73: As indicated in the main comments, this is one point that could be used for\nhighlighting the relevance of the presented study\nPoint-by-point response to referees’ comments\nReply 6: We did not see any comments relevant to line 72/73 in the main text. We agree that the point of\nthe sentence on line 72/73 is the goal of our study that we emphasized on the first and last paragraphs of\nthe main text.\nComment 7: L93: provide short information why it is the preferred salinity model\nReply 7: In geophysics, the non-uniqueness in inverse problems is well-known, in which an infinite\nnumber of models can fit a specific dataset. Thus, we applied regularized and Bayesian inversion methods\nto quantitatively estimate the non-uniqueness in the resistivity models. Thus, the estimated salinity\nderived from the resistivity inherits non-uniqueness from bulk resistivity models. The second reason is\nthat we used a range of values for porosity, temperature, and cementation factor to estimate salinity and\nonly showed a model with one set of variables which we think are the best available constraints.\nTherefore, we used the term “preferred resistivity and salinity models” in the main text. We provided\nvariations of these variables as extended data figures in the Supplementary materials.\nComment 8: L93: define why 3 psu is chosen (and not 1)\nReply 8: on lines 473-480, we explained why 3 psu is chosen. We used a combination of 10 Ωm and 3 psu\nas references to delineate freshwater zones. Converting bulk resistivity to salinity by using Archie’s law\n(1942) involves temperature, pressure, porosity, and cementation factor. These variables can vary in\npractice, so we only provide the estimation of salinity from our resistivity models.\nAlso, we notice that MT data are unable to resolve the magnitude of highly resistive layers due to\nresponse saturation. This is reflected where the Bayesian 95th percentiles run up close to the prior upper\nbound of 1000 Ωm (Extended Data Figure 2). Thus, MT responses might underestimate the true\nresistivity of the freshwater bodies, and the actual salinity of the fresh aquifers can be lower than 1 psu,\nwhich is the threshold commonly used for drinking water. 3 psu contour in the salinity model is the upper\nlimit of the estimated salinity of fresh groundwater.\nWe added more details to that discussion on lines 443-453 to clarify our point.\nComment 9: L98: specify which isotopes or at least state if sediment or water is meant here\nReply 9: We clarified in the text that the isotopic data is for sediments.\nFor example, Goodbred et al. (2014) used 87Sr/86Sr, radiocarbon ages, and other chemical signatures to\ncharacterize the origin and depositional timing of sediments.\nWe only used the information from the cited papers to show the spatial variations of the Holocene-\nPleistocene sediment boundary.\nComment 10: L98/99: In Fig. 3 150-250m are given (instead of: 90m), or is another figure meant here?\nReply 10: on lines 98/99, we referred to hand tubewells data for lithology. 90 meters is the deepest depth\nthat the hand tubewells extend down to. The information from tubewells allows mapping Holocene-\nPleistocene sediment boundary, the LGMP, and Holocene sediment thickness. In Fig. 3, we used contours\nof 10-50-90 meters to depict Holocene sediment thickness and show LGMP distribution as orange colors.\nThe given depth ranges from 150 to 250 m are associated with the freshest part or lowest salinity at our\nMT sites (the square symbol).\nPoint-by-point response to referees’ comments\nIn Figure 3, we added the contour symbol to the legend to clarify Holocene sediment thickness and\nrelated the 150-250 m depth ranges with salinity information to Fig. 2b.\nComment 11: L101-103: connection to onshore aquifers?\nReply 11: Resistor R1 and R2 are onshore aquifers located within the low-lying coastal GBD. They are\nnot offshore aquifers. We revised in the main text to clarify that R1 and R2 are potentially connected to\ngroundwater recharge areas further north (Pleistocene terrace) and the West Bengal, respectively.\nAlthough R1 and R2 can connect to upgradient recharge systems, the extremely low topography in coastal\nGBD and low hydraulic gradients requires long timescales for natural horizontal recharge. The LGMP\nprotects Pleistocene groundwater from vertical mixing with shallow Holocene groundwater. The modern\nvertical recharge can happen where LGMP is terminated, or at the boundary between paleovalley and the\ninterfluve. Therefore, we proposed that R1 and R2 are primarily Pleistocene fresh groundwater.\nPlease see the revised “Fossil fresh aquifers” section in the main text for more information.\nComment 12: L12: Fig. 4 a and b cover depths up to 150 m, R1 and R2 fresh water extend (much) deeper\nand only the topmost parts overlap Fig 4a and b, this should be mentioned more clearly when deriving\nanalogies and processes\nReply 12: We agree with the reviewer. The Fig. 4a&b only show depth down to 150 m to explain the\nformation of paleovalley and LGMP distribution. Thus, we did not extend further down to the bottom of\nR1 and R2.\nWe added this clarification to the end of the Fig. 4 caption.\n“Note that Fig. 4a&b only shows depths down to 150 m to explain the formation of paleovalley\nand LGMP distribution, whereas Fig. 4c extends deeper to capture the full dimensions of two\nfresh aquifers, R1 and R2.”\nComment 13: L112: why are the deductions made from uncorrected Radiocarbon data? (especially as\nthis (published) data is used as a strong argument here for fossil water in the investigated system)\nReply 13: In Khan et al. (2019), there are two deep groundwater ages: As-affected areas and low-As\nareas. For the As-affected areas, the younger age of 2.0 +/- 0.6 kyrs might indicate the downward flow of\nHolocene As-contaminated groundwater rather than the true age of deep groundwater (>150 m). The As-\naffected areas do not have the LGMP to prevent vertical mixing and are associated with the Ganges\npaleovalley areas. For the low-As areas, the older age of 12.0 +/- 4.0 kyrs reflects no or less As\ncontamination. The deep groundwater in these areas is protected by LGMP from vertical mixing or\ndownward flow and is associated with the interfluve areas. Therefore, they showed groundwater ages of\nwater samples similar to aquifers in our study areas. Although the location of their study is not collocated\nwith our survey area, the two transects in their study go across the Ganges paleovalley located to the\nnorthwest of our MT survey.\nIn the supplementary materials of Khan et al. (2019), they have a section to explain why they did not use\nthe 14C age correction. The correction applied for the uncorrected ages would reduce the age by about 1.0\nPoint-by-point response to referees’ comments\n14C kyr, but they would have high uncertainty. Table S1 in Supplementary Materials of Khan et al. (2019)\nhave the corrected 14C age of groundwater for the southern transect.\n“14C Age correction: Converting radiocarbon ages to calendar ages is notoriously difficult for DIC in\ngroundwater because several poorly constrained corrections can subtract hundreds to a few thousand\nyears [Geyh, 2000]. One method is to infer an initial radiocarbon age from the youngest sample in a data\nset that does not contain any detectable bomb-produced 3H and therefore no bomb-produced 14C either\n[Hoque et al., 2012, Mihajlov et al., 2016]. Using existing dataset, this would amount to a decrease in the\nuncorrected ages by about 1.0 14C kyr (Supporting Figure S2). Additional corrections reflecting the\nrelease of carbon along the flow path have been proposed [Hoque et al., 2012, Mihajlov et al., 2016]\nconsidered on the basis of the stable isotopic composition (13C/12C) of DIC but they are highly uncertain.\nThe average δ13C of -4.5‰ (σ = 2.2‰, n= 33) in deep groundwater DIC within the study area shows\nconsiderably less depletion than deep groundwater in other parts of the basin, with an average δ13C of -\n16‰, (σ = 6‰, n= 39) (Figure 2c). We do not have sufficient data for the detailed geochemical modeling\nrequired for correction in such a case [Aravena et al., 1995] and no further corrections were attempted.”\nKhan et al. (2019) is not the only evidence that we use to prove that these aquifers are Pleistocene in age.\nStratigraphy and LGMP provide geological constraints. Multiple references cited in this work infer deep\naquifer systems in coastal BLD with similar depth range of our study have Pleistocene in age. We also\nadded the 14C activities from Majumder et al. (2011) in the discussion. Khan et al. (2019) also included\ngroundwater ages from other papers to make comprehensive analysis. Therefore, we used this published\ndata to provide geochemical constraints on groundwater ages in our study.\nIn the main text, I revised the section using the uncorrected age from Khan et al. (2019). I added more\nexplanations where the Pleistocene groundwater age is applicable. See the revised “Fossil fresh aquifers”\nsection for more details.\nComment 14: L106-118 preserved Pleistocene groundwater vs. resemblance to groundwater system\noffshore New Jersey: It should be clarified, if only the morphology or the (emplacement) processes are\nexpected to be similar\nReply 14: The point when comparing the morphology between these two aquifers is to show that R1 is\npotentially connected to upgradient recharge (albeit which would be extremely slow recharge due to the\nlow topography gradient) and not isolated pockets. The emplacement processes of the two aquifers might\nbe similar; however, the sediment supplies and geological settings of the two deltas are different. Adding\nthis complexity to the comparison will lose track of interpreting R1. Like our response to direct comment\nin the main text, we decided to neglect the comparison between the preserved Pleistocene groundwater in\ncoastal BLD with groundwater systems offshore New Jersey to keep track with the length of the paper\nand not losing focus on the main aquifer R1.\nComment 15: L127: refer to Fig. 2\nReply 15: We added Fig. 2 to the end of the sentence.\nComment 16: L152: Sea-Level Change\nReply 16: We fixed the typo.\nPoint-by-point response to referees’ comments\nComment 17: L302: refer to Fig. 1c\nReply 17: We referred to Fig. 1c on line 302.\nComment 18: L456: add link\nReply 18: We will make all data publicly prior to publication.\nFigures and Tables:\nComment 19: Are Figures 4a and 4c and Table 1 mentioned in the manuscript text?\nReply 19: We referred figures 4a and 4c to appropriate places in the main text. Previously, we used them\ntogether as Fig. 4.\nFor table 1, I did not mention it in the main text. Table 1 could go together with Fig. 2. Table 1 lists the\nbulk resistivity and estimated salinity of features in Fig. 2a and 2b. Although the color scales in Fig. 2a\nand 2b indicate clearly the values, table 1 provides the upper and lower limit of each feature. To save\nspace in the published version, we have eliminated this table.\nComment 20: L 470: please delete this line\nReply 20: We deleted this line.\nComment 21: Fig. 1: legends in a and b are too small, an in b the red entries are difficult to see\nReply 21: They are elevation with color-coded scale. The resolution of the figure is not great in the word\ndocument. We enlarged the text of the legends of Fig. 1a&b in the revised manuscript. We enlarged the\nlegends of Fig. 1a and 1b. We will upload the high-resolution figures with the revised manuscript.\nComment 22: L475: not clear how subaqueous delta depth an SoNG can be distinguished\nReply 22: The Swatch of No Ground (SoNG) is color-coded as blue, sea and lake, which has elevation\nless than 0 m. The depth of SoNG is based on bathymetry. If you think that the information is irrelevant,\nwe can erase it.\nComment 23: L481 not clear where shrimp farms are (dark blue to purple appears a bit masked)\nReply 23: The color of shrimp farms (dark blue to purple) is clear due to its high contrast with the pink\ncolor of fallow agricultural fields and green color of vegetation. Channels and river mouths have blue\ncolors, but their pattern can be distinguished from shrimp farms. Note that shrimp farms are primarily\nnext to channels.\nPoint-by-point response to referees’ comments\nComment 24: Fig. 2: The axis gives distances (+-60 km), where P11 appears to be the center. Why did\nyou choose this scale (and the center at P11) instead of 0 to 120 km? And why are Sundarbarns explicitly\nmentioned here?\nReply 24: P11 is very close to the center of the profile, but it is not at the center. This scale was used\nwhile we did inversion models for resistivity. The default setting of the code chose the center of the\ninversion profile as the origin. Readjusting the distance scale makes no difference in absolute values. We\nadded the Sundarbans here to help readers know the locations of the groundwater bodies and the saline\nwater gap. Readers can use the map in Fig. 1c and 3 to visualize the location of these structures. We added\nspecific four stations and North-South symbols with the same purpose.\nComment 25: Fig. 3: connection between Fig. 3 and Fig. 4 somewhat difficult\nReply 25: Figure 3 shows all available data that we gathered to explain the groundwater distribution in the\nmap view. Figure 4 shows a conceptual model to explain the formation of deep fresh groundwater R1 and\nR2 due to sea-level fluctuations with time in a cross-section view. The two figures have different\npurposes. You can use the Ganges paleovalley to connect the two figures in spatiotemporal scale.\nComment 26: Fig. 4: maybe insert R11 and R2 into Fig. 4c\nReply 26: We inserted R1 and R2 into Fig. 4c.\nResponse to Reviewer #3:\nComment 1:\nThe carbon-14 data used to support the conceptual model lacks geographic specificity, which leads to\napparent overinterpretation of the data. The data comes from one study (ref #5, Kahn et al. 2019). Most of\nthe carbon-14 ages in Kahn are Holocene and the subset that are Pleistocene are not identified and\nlocated relative to the current study area. Some of the points in Kahn are more than 100km from the\nclosest boundary of the study area. Extrapolating over such a distance is unreasonable in most contexts\nand certainly in a major delta complex such as the Ganges. Further, the authors argue for the water in R1\nand R2 to have been emplaced during the last glacial maximum -- only 7 of the carbon-14 dates in Kahn\nfall in the range of the LGM, which makes understanding the location of Kahn's data relative to the\ncurrent study even more crucial.\nConfidence in Pleistocene ages for the vast majority of Kahn's samples nearest to the subject study area is\nkey. If they're Holocene, the entire conceptual model falls apart. To invoke Pleistocene emplacement of R1\nand/or R2, the authors must unequivocally demonstrate the presence of Pleistocene water in close\nproximity to their transect.\nReply 1: In Khan et al. (2019), there are two deep groundwater ages: As-affected areas and low-As areas.\nFor the As-affected areas, the younger age of 2.0 +/- 0.6 kyrs might indicate the downward flow of\nHolocene As-contaminated groundwater rather than the true age of deep groundwater (>150 m). The As-\nPoint-by-point response to referees’ comments\naffected areas do not have the LGMP to prevent vertical mixing and are associated with the Ganges\npaleovalley areas. For the low-As areas, the older age of 12.0 +/- 4.0 kyrs reflects no or less As\ncontamination. The deep groundwater in these areas is protected by LGMP from vertical mixing or\ndownward flow and is associated with the interfluve areas. Therefore, they showed groundwater ages of\nwater samples similar to aquifers in our study areas. Although the location of their study is not collocated\nwith our survey area, the two transects in their study go across the Ganges paleovalley located to the\nnorthwest of our MT survey. Based on the similarity of lithological setting, we used Khan’s study to\nconstrain the imaged fresh aquifers R1 and R2.\nIn the supplementary materials of Khan et al. (2019), they have a section to explain why they did not use\nthe 14C age correction. The correction applied for the uncorrected ages would reduce the age by about 1.0\n14C kyr, but they would have high uncertainty. Table S1 in Supplementary Materials of Khan et al. (2019)\nhave the corrected 14C age of groundwater for the southern transect. The SWAA40557 sample has an\nuncorrected 14C age of 29.647 kyrs at 249 m depth and the corrected 14C age of 28.776 kyrs. The\nKHU40053 sample has an uncorrected 14C age of 15.118 kyrs at 162 m depth and the corrected 14C age of\n14.247 kyrs.\nWe agree with the reviewer that having deep groundwater age collocated with our MT survey will provide\nbetter constraints, but currently we do not have groundwater ages of R1 and R2. Thus, we used other\nevidence to constrain the imaged groundwater R1 and R2 (Fig. 3). The radiocarbon age from Khan et al.\n(2019) is one of the constraints. Khan et al. (2019) also included groundwater ages from other papers to\nmake comprehensive analysis. Along with the radiocarbon age, stratigraphy and LGMP provide\ngeological constraints. Multiple references cited in this work infer deep aquifer systems in coastal BLD\nwith similar depth range of our study have Pleistocene in age.\nWe are aware of the potential vertical mixing between shallow Holocene and deeper Pleistocene\ngroundwater due to the absence of LGMP. In these circumstances, deep groundwater may have Holocene\nages. For instance, Khan et al. (2019) showed that the ages of groundwater samples located within the\npaleovalley-areas with no LGMP-the groundwater age is Holocene. Also, the LGMP is laterally\ndiscontinuous across GBD that will lead to the variations of fresh groundwater distribution. Thus, our\nconceptual model still works on a regional scale.\nIn the main text, we revised the section using the uncorrected age from Khan et al. (2019). We added\nmore explanations where the Pleistocene groundwater age is applicable. See the revised “Fresh fossil\naquifers” section for more details.\nComment 2:\nIn the development of their conceptual model, the authors argue that the fresh water identified in R1 and\nR2 is relict, connate water emplaced during the LGM, as opposed to the terminal portion of an active\ngroundwater system. At several times they contradict this conceptualization by referencing upgradient\nconnectivity and making comparisons to the Ganges system to an offshore system in New Jersey (specifics\nnoted in comments in PDF file). The contradictory statements in the text undermine their conceptual\nmodel. Further, the northern extent of the R1 freshwater region extends an unknown distance outside the\nstudy area. There is no evidence presented for the reader not to envision this lobe of freshwater\ncontinuing upgradient as part of a continuous, active groundwater flow system.\nThe deficiencies in this study principally stem from its heavy reliance on data from existing published\nPoint-by-point response to referees’ comments\nwork in areas that do not exactly overlap with theirs to support their interpretations of their interesting\nMT transect results. Without site specific data to support their interpretations, the authors often find\nthemselves at odds with the data and with their proposed conceptual model. The original work (the MT\ntransect) appears sound and compelling; its interpretation much less so.\nReply 2: We argued that R1 and R2 are primarily fossil freshwater emplaced during the sea-level\nlowstand and potentially connected to upgradient recharge. Based on morphology of aquifers, LGMP\npaleosol map, trace of Ganges paleovalley, R1 and R2 could connect with upgradient recharge sources\nlocated over tens to hundreds of kilometers away from location of R1 and R2. Although they might be\nconnected to recharge sources, the recharge rates are extremely low due to low topography and hydraulic\ngradients. The recharge rates require thousands of years or more. Therefore, these aquifers can be\nconsidered as non-renewable water resources. We added references and revised the “Fossil fresh aquifers”\nsection in the main text to clarify our argument. We neglected the New Jersey comparison to focus on the\ninterpretation of R1 (see the responses to these specific comments).\nWe agree that R1 could extend to the northern part out of the MT profile.\nWe revised the interpretation based on all the comments. Please check the revised edition for more details.\nAgain, the radiocarbon age is not the only data we used to interpret our resistivity and salinity models.\nWe also used LGMP outcrop and subcrop as well as Holocene sediment thickness, the trace of\npaleovalley to explain our geophysical results. We also added the 14C activities from Majumder et al.\n(2011) in the discussion. The compelling parts of our work are imaging the full depth extent or thickness\nof groundwater systems rather than the top of the aquifers, the presence of LGMP controlling the\nfreshwater distribution, and the past sea-level cycles contributing to deep groundwater emplacement. The\nproposed conceptual model helps visualize these key results and dynamic processes involved in deep\ngroundwater emplacement.\nDirect Responses to Comments in the Manuscript:\nComment 1: Line 29-30: reference needed\nReply 1: We added van Geen et al. (2008) and Rahman et al. (2020) papers.\nComment 2: Line 50: “S” is capitalized in caption of fig 1\nReply 2: We changed small “s” to capital “S”.\nComment 3: Line 53: This colloquial sounding sediment size fraction seems out of place with the specific\nsize-fraction terms used in the rest of the sentence and prior sentence.\nReply 3: “coarse sand” and “clayey silt” are used to describe variety of grain size: coarse sand (0.5 - 1\nmm) to silt (0.05 - 0.002 mm) and clay (<0.002 mm). “Clayey silt” represents major grain size in the\nSundarbans mangrove forest, which is silt-dominated mixture with sand.\nComment 4: Line 55: Not on map. The reader has no context for where these are. Are they important to\neven mention?\nPoint-by-point response to referees’ comments\nReply 4: The two terraces are on Fig. 1b. I put “(Fig. 1b)” in the main text to show them. We think it\nworths to mention them because they have LGMP outcrops compared to LGMP subcrop, which is\nrelevant to Pleistocene fossil groundwater. They could be the location where groundwater recharge occurs\nfor Pleistocene aquifers, although the timescales could be thousand years or more.\nComment 5: Line 57: consider hyphenating: fresh-to-saltwater\nReply 5: We added hyphens.\nComment 6: Line 58-59: The statement about compaction appears unsupported. There are refs before it\nthat seem to document high sedimentation rates and refs at the end of the sentence that seem not to be\nrelated to the topic of compaction. Please clarify.\nIs there potential for mass loading (independently of compaction) to contribute to subsidence?\nReply 6: We included the subsidence to provide background information relevant to the study area. The\nsubsidence due to compaction contributes to the depth of Holocene-Pleistocene sediment boundary.\nThis boundary separates shallow Holocene groundwater, which has multiple problems mentioned above,\nto Pleistocene groundwater, which is fresh and potentially As-free.\nAlso, the subsidence creates the accommodation space for sediments and affects the evolution of river\nsystems in the study area.\nReferences at the end:\n28. Allison et al. (2003) indicate that the regional subsidence of Bengal Basin affects the evolution of\nHolocene sediments in the delta.\n29. Grall et al. (2018) show the subsidence increase in seaward direction.\n30. Steckler et al. (2022) provide quantitative evidence of compaction and subsidence. Subsidence and\ncompaction vary with depth, space, and time.\n31. Michels et al. (1998) support delta front seaward progradation. Label F in Fig. 4b demonstrates the\ndelta progradation.\nYes, mass loading plays a role in subsidence. In our study area, young unconsolidated sediments have\nhigh porosity. Therefore, adding mass on top of older sediments will cause compaction to occur.\nComment 7: Line 60: interesting but not related to the rest of the sentence and could be dropped with no\nloss if you need to save space or expand elsewhere.\nReply 7: We agree. So, we erased the phrase.\nComment 8: Line 63: Reference needed\nReply 8: We added Umitsu (1993) as reference for Pleistocene sea-level change in coastal BLD.\nPoint-by-point response to referees’ comments\nComment 9: Line 64: Misleading. This statement appears to be expanding on the statement in the topic\nsentence, which specifically mentions the GBD. Neither of the two references supporting this statement\nare from the GBD. Rephrase to clarify this is a generic world-wide statement, or cite supporting\ndocuments specific to GBD.\nReply 9: We changed references from global scale to local GBD by using Umitsu (1993),\nGoodbred and\nKuehl (2000), and Goodbred et al. (2003).\nComment 10: Line 101-103: This is speculative and seems plausible, however nothing in the remainder\nof this section definitely supports the statement.\nIt seems the water in both R1 and R2 could be part of modern flow systems that originate upgradient.\nSpecifically, R1 maintains its full thickness all the way to the left (north side) of figure 2. The upper limit\nand connectivity is unknown. On line 120, the authors acknowledge the seeming disconnection of R2 is\ndue to the geometry of their transect and the water is related to flowpaths originating to the NW of R2\nand off the N-S orientation of their MT survey.\nReply 10: The LGMP protecting the upper boundary of two aquifers from vertical mixing is a primary\nfeature to separate R1 and R2 from Holocene (modern) flow systems.\nBoth R1 and R2 could connect to upgradient recharge systems, but the extremely small gradient of low-\nlying coastal GBD requires timescales of thousands of years or more for recharge. Wilson and Goodbred\n(2015) (ref. 48) show that the flatter lower delta plain has gradient of 10-5, which is 1 meter change in\nslope over 100 km. Thickness of R1 extending towards the northern side at depth means Pleistocene\nfreshwater can extend farther north of the profile. However, the timescales required for recharge in our\nsurvey areas is not associated with modern flow systems.\nCurrently, we do not have isotope data at R1 and R2 to constrain the age and flow systems of the two\naquifers. Thus, we gather other information to explain the distribution of R1 and R2.\nThe Holocene isopach is one of the geological indicators for the age and extension of groundwater. We\nindicated in Fig. 3 that the Holocene sediment thickness near the vicinity of R1 is 90 m. Therefore, we\nstated that the aquifers here are likely associated with Pleistocene or older strata. The LGMP created a\nnatural barrier to prevent vertical mixing between modern (Holocene) groundwater with paleowater.\nComment 11: Line 106: Ref needed.\nReply 11: We added Sarkar et al. 2009 and Pickering et al. 2014.\nComment 12: Line 109-111: ref needed.\nReply 12: We moved the Hoque et al. (2014) and McArthur et al. (2008) to the end of the sentence.\nComment 13: Line 109: P-M reversed in acronym\nReply 13: We changed the order.\nComment 14: Line 111-112: unsupported speculation\nPoint-by-point response to referees’ comments\nReply 14: The LGMP is the fine-grained transgressive sediments mentioned above. The paleo-fresh\ngroundwater could not vertically mix with Holocene groundwater due to LGMP. I added McArthur et al.\n(2008).\nComment 15: Line 112-113: Kahn et al. clearly show the presence of some Pleistocene groundwater GW\nat 200-300m in the Ganges delta, however most of their dates are Holocene, including a large\nconcentration of dates around 10ka. It is not straightforward to identify the Pleistocene samples in Kahn\net al. that are within or reasonably proximal to the present study.\nBecause the authors of this study rely on the age dates in Kahn et al. to support their conceptual model\n(they are key to the assertion that the R1 and R2 are relict Pleistocene water), the onus is on the authors\nto equivocally demonstrate the Pleistocene samples in Kahn et al. are within the current study area (or\nmeaningfully proximal). It should not be up to the reader to piece this together. If it cannot be\ndemonstrated that there is Pleistocene water in or very near their study area, the conceptual model needs\nto be revised.\nReply 15: We do not have isotope data collocated with the MT survey profile. Thus, we used Khan et al.\n(2019) as a reference to provide constraints on groundwater ages. Since their samples cover both the\ninterfluvial and paleochannel areas, we can use them to correlate with LGMP and infer groundwater ages\nclose to our survey. Khan et al. also include other studies to compare and investigate the deep flow system\non a regional scale.\nKhan et al. (2019) is not the only evidence that we use to prove that these aquifers are Pleistocene in age.\nStratigraphy and LGMP provide geological constraints. Multiple references cited in this work infer deep\naquifer systems in coastal BLD with similar depth range of our study have Pleistocene in age. Khan et al.\n(2019) also included groundwater ages from other papers to make comprehensive analysis. Numerical\nmodeling from Michael and Voss (2008) indicates low regional horizontal recharge.\nIn Khan’s paper, there are two deep groundwater ages for As-affected areas and low-As areas.\nFor the As-affected areas, the younger age (2.0 +/- 0.6 kyrs) might indicate the downward flow of\nHolocene As-contaminated groundwater rather than the true age of deep groundwater (>150 m).\nFor the low-As areas, the older age (12.0 +/- 4.0 kyrs) reflects no or less contamination and can be a true\nage. The groundwater ages of these areas are more relevant to fresh aquifers in our studies since vertical\nmixing or downward flow is prohibited by LGMP. Mihajlov et al. (2016) show that deep aquifer (>150\nm) in central Bengal Basin was last substantially recharged around ~10 kyr ago. The last substantial\nrecharge at ~10 kyr was likely to occur in coastal Bangladesh and might be a reason for a large\nconcentration of dates around 10 ka. Therefore, the true deep groundwater ages must be older, and the 10\nka only shows the last time of recharge.\nMost low-As samples in Khan et al. have uncorrected 14C ages spanning within 10-20 kyrs and two\nsamples up to 30 kyrs. In their supplementary materials, they stated that the corrected 14C age may\ndecrease by 1.0 kyrs compared to the uncorrected ages. More correction methods can be applied, but they\nalso have high uncertainty. Table S1 in Supplementary Materials of Khan et al. (2019) have the corrected\n14C age of groundwater for the southern transect. The SWAA40557 sample has an uncorrected 14C age of\n29.647 kyrs at 249 m depth and the corrected 14C age of 28.776 kyrs. The KHU40053 sample has an\nuncorrected 14C age of 15.118 kyrs at 162 m depth and the corrected 14C age of 14.247 kyrs.\nPoint-by-point response to referees’ comments\nComment 16: Line 115: I'm not sure this is a useful comparison. The cited study indicates the\nPleistocene water is found in fine-grained units separated by saline water in higher-permeabilty units.\nThe model proposed by the authors herein on lines 101-103 seems to envision a pocket of relict\nPleistocene water that stands in isolation to the modern flow system. They imply this relict water is in a\npermeable unit (otherwise it could not be a developed as a drinking water resource). However, the notion\nof it being connected to modern upgradient systems arises later, including here, by leaning on this NJ\nstudy and later on lines 119-121. However, even after recognizing the potential continuity with\nupgradient systems, a static period of formation during the LGM is described on line 122.\nConnection with the upgradient GW system also is indicated by the morphology of R1 at the northern end\nof the transect -- there is no indication of a terminus and indeed looks very much like the\nnearshore/offshore terminus of the on-shore deep freshwater GW system described in the NJ article and\nin many other coastal GW systems.\nConsistency and clarity is needed in the authors' conceptualization of the system, and its alignment with\nexisting data in the Ganges delta.\nReply 16: First, we show the similarity in morphology of R1 in this study and New Jersey offshore\naquifer to argue that R1 might have upgradient recharge at very low rate due to low topography. Since we\ndo not have a full image of entire aquifer related to R1-reviewer points out that R1 maintains its thickness\nto the north-using this similarity might help clarify the point. Isolated pockets of freshwater will have\ndifferent shapes compared to the connected ones.\nIn New Jersey system, the freshwater stored in Miocene fine-grained sediments. We are aware of the\ncomplexity and differences between two delta systems in terms of stratigraphy, sediment supplies, and\ngroundwater ages. That’s why we only consider the shape of groundwater bodies.\nTo keep track of the length of the paper and not losing focus on the main aquifer R1, we neglect the\ncomparison with the New Jersey system.\nComment 17: Line 122: This is a significant overinterpretation of the age dates in Kahn et al that are\nbeing used to support a Pleistocene water source. Taking 20-26ka as the duration of the LGM, only 7 of\nthe wells in Kahn's dataset correspond to LGM and the location of those 7 wells is not provided to the\nreader of this article. As noted earlier, the onus is on the authors here to demonstrate those LGM samples\nare spatially relevant to their arguments.\nReply 17: We agree that we have minor constraints on R2 in terms of groundwater ages and LGMP\nextent. We changed the suggested timing of formation prior to the Holocene due to the dimensions of R2\nand stratigraphy. It is likely to have Pleistocene in age, but no further constraints on that.\nComment 18: Line 127-129: The traverse and paleochannel identified by Hoque is considerably north of\nthe end of MT transect herein. Their transect is about 23.3° latitude and the proposed paleochannel\nlocation shown on figure 3 herein is at about 22.3° latitude.\nPerhaps reference 43 (Sincavage et al. 2017) provides information on a more southerly channel but I was\nunable to access the article.\nReply 18: We traced the Ganges paleochannel further down to Sundarbans area in Fig.3 using its\nthalwegs. The transect of the Hoque et al. (2014) cover the interfluve and paleochannel areas. We had no\nPoint-by-point response to referees’ comments\nproblem with accessing Sincavage et al. (2017), so we can share this paper with the reviewer if necessary.\nThe study area of Sincavage et al. (2017) is in Sylhet, not near the southerly channel.\nComment 19: Line 133-137: Speculative\nReply 19: Yes, we clarify that we speculated this idea.\nComment 20: Line 148-151: Assigning an early Pleistocene data based on imaging alone is highly\nspeculative. Why could it have not occurred during a subsequent interglacial or even during the earlier\nPliocene? No data is provided to constrain. Even if the interpretation of R1 and R2 is correct, a LGM\ndate for them does not provide any sort of constraint on how long before their emplacement the saline\nwater in C2 must have occurred.\nReply 20: We revised the hypothesis of the emplacement of conductor C2. However, we do not have\ndirect constraints from boreholes, geophysical or geochemical data for C2. Therefore, we used offshore\nSangu-1 well and seismic data to constrain C2.\nWe agree that we do not have any absolute age constraint for the marine deposition. Thus, we change to\nrelative timing of events based on spatial relationship between C2 and R1 & R2.\nComment 21: Line 155-157: neither statement is satisfactorily substantiated\nReply 21: we added a combination of evidence to clarify the statement. More details are in the “Fossil\nfresh aquifers” and “High salinity zones” sections.\nComment 22: Line 289: I read through this but am not a geophysicist and the content was outside my\narea of expertise to critique.\nReply 22:\nComment 23: Line 386-387: Perhaps I misunderstand the context of the term “geothermal gradient” as\nused here. As I have experience it’s use, it refers to the increase in temperature with depth due to\nconductive hat movement out of the earth. In this context of this sentence, it is implied that the geothermal\ngradient is somehow related to air temperature.\nReply 23: The geothermal gradient here accounts for temperature increase with depth. Groundwater\ntemperature is potentially higher at depth due to heat transfer from surroundings. Therefore, the\ntemperature of groundwater is higher than the surface water temperature, which is between 25 and 28\ndegrees Celsius. Adding geothermal gradient factor to the conversion prevents overestimating fluid\nsalinity at deeper depth.\nComment 24: Line 389-390: This would seem to pose a problem for the presence of low-permeability\nunits that have been proposed and invoked as significant impediments to flux of Pleistocene water out of\nzones R1 and R2.\nReply 24: The proportion of sand is higher in Pleistocene sediments than in Holocene unit in coastal BLD\nshould not pose a problem for pumping groundwater out of the Pleistocene aquifers. Multiple deep wells\nhave been drilled to support local communities near Khulna city and elsewhere.\nPoint-by-point response to referees’ comments\nComment 25: Line 390-391: unusual phrasing\nReply 25: We changed “dominates” to “is higher than”\nThis phrase describes the conductivity effect of saltwater is the first-order control of the bulk\nconductivity. The conductivity effect of clay is lower than that of saltwater. Therefore, Archie’s law can\nbe applied to estimate fluid resistivity/conductivity without violating its assumption.\nComment 26: Line 390: “deeper” would perhaps be a bit more precise\nReply 26: We changed as suggestion.\nComment 27: Line 394: provided\nReply 27: We fixed it.\nComment 28: Line 492: define term for those unfamiliar\nReply 28: We added text next to PSS78 to clarify.\nComment 29: Line 495: Please make clear if this comment refers to both figures 2a and 2b. If only 2b,\nplease clarify which two scales it refers to as I interpret it as having just one scale ranging from yellow to\nblue.\nReply 29: We added text next to PSS78 to clarify."}, {"round": 2, "reviewer_comments": {"reviewer_2": "I thank the authors for thoroughly revising the manuscript. It has significantly improved, and many aspects have been clarified. Also, comments have been addressed in detail. Yet, though in discussions absolute age derivations have been largely changed to relative timing where direct constraints were missing, regarding uncorrected 14C ages I still would suggest to incorporate some (or a concise version) of the explanation into the manuscript or supplement that the authors provided in their reply to Comment 13 (reviewer 2) and “Comment 15 of the Direct Responses to Comments in the manuscript”, especially regarding the Khan et al. study. Also, highlighting the uniqueness of their study as suggested in the first revision was moderately addressed, which might still be slightly improved. As a minor point I suggest to check for spelling/typos (especially the revised passages).", "reviewer_3": "This is a review of the revised manuscript titled \"Deep fossil freshwater in coastal Bangladesh\" by Huy Le and others. The revised manuscript largely addressed comments I had and those of the additional reviewers. I appreciate the authors' improvements to the manuscript, in particular the \"Fossil Fresh Aquifers\" section. This section provides better background for someone not familiar with the impressive body of work on the GBD and more cohesive support for the authors' interpretations (despite the few remaining concerns about the emplacement mechanism and age of R1 discussed below).\nI reread and reviewed the entire manuscript for technical content but stopped commenting on editorial issues after the first 50 lines or so. There were several instances of either awkward phrasing or missing words. An editorial review prior to resubmission is recommended. The following items were noted for consideration by the authors.\nL22\nconsider adding \"the\" before \"location\"\nL34-35\nconsider \"...deep groundwater at a regional scale...\"\nL35\nconsider changing \"sustainability\" to \"susceptibility\"\nL69-70\nThe topic sentence on L64 limits important sea-level fluctuations to the Pleistocene\nL108-193\nI can accept R2 being relict Pleistocene water as it appears to be cut off from the modern flow system by the upgradient paleochannel and is more saline than R1. I have more difficulty with the certainty expressed in the age of R1 and its origin ascribed to a body of water \"formed\" during the Pleistocene. I think the issue is largely semantics.\nFirstly, the origin of the water. The language used in the text implies connate water, although that term is no longer used in the manuscript, e.g:\n\"freshwater bodies formed during the Last Glacial Maximum’s (LGM) sea-level lowstand and subsequently sealed by a finer-grained, low-permeability layer\" [L108-110]\nPoint-by-point response to referees’ comments – Second Revision\n\"remnants of preserved Pleistocene groundwater\" [L180]\nIntentionally or not, the language implies the water in R1 originated in that location and has remained there for 12,000+ years. Perhaps it has, however, that has not been demonstrated nor is it not supported by subsequent statements that indicate R1 is the distal end of a larger upgradient flow system (L123-125) with slow transport times (L123-125).\nSecondly, the timing of recharge of the groundwater in R1. The location of the 14C ages supporting the assertion of Pleistocene water is provided only as \"northwest of our profile\". Other than the fact those samples are somewhere in the western portion of the 100,000 km2 GBD (L50), the reader has no information whether these data are relevant -- were they collected 1 km, 10 km, or 100 km away? There is a significantly stronger case for transferability if they're nearby. Consider showing the locations on Figure 3 if they are within the map extent. The transport time of \"thousands of years or more\" (L124-125] is consistent with Holocene as well as Pleistocene recharge. The 14C data presented provide evidence for shallow-deep difference in age in another part of the GBD system having similar characteristics to the region of R1, which is useful, but not conclusive given the lack of information about their location. The evidence presented for Pleistocene water in R1 is not sufficiently compelling to make a definitive assertion as to its age, particularly in light of the wide range and complex distribution of groundwater ages presented by Kahn et al. (2019), from which the example 14C data were taken.\nConsider using the uncertainties surrounding the emplacement mechanism and timing of R1 in your favor instead of trying to force certainty with imperfect data. For example, if R1 is the distal end of a long, slow flow system, it is potentially a more usable and renewable resource than connate Pleistocene water. If the age of R1 is mid-Holocene (for example) rather than Pleistocene, it could indicate a more useful/renewable source. Since the groundwater system hosting the freshwater in R1 is relatively coarse-grained (capped by the clayey paleosol), could you pump water from R1 without substantial drawdown?\nCould you potentially induce a bit more recharge from modern recharge zones or set up managed aquifer recharge systems? The MT transect is a really useful first step to asking the next set of important questions relating to development and management of the freshwater resource in the GBD"}, "author_response": "Point-by-point response to referees’ comments – Second Revision\nDate: June 28, 2025\nManuscript Number:\nTitle of Article: Hidden deep freshwater reserves beneath salinity-stressed coastal Bangladesh\nName of the Corresponding Author: Huy Le\nEmail Address of the Corresponding Author: hdl2115@columbia.edu\nDear All,\nWe would like to thank the two reviewers for their second round of feedback and the editor for the time\nextension to edit our manuscript. We have incorporated the reviewer’s comments to revise the discussion\nof deep freshwater systems in coastal Bangladesh. Below we provide the reviewer comments in black text\nand our replies in blue.\nOnce again, thank you all for your consideration to our manuscript.\nBest regards,\nHuy Le and coauthors\nMAJOR REVISION\n Modify the interpretation of fossil freshwater to partial recharge freshwater reserve for R1\n Clarify the 14C age profile in Figure 3 and re-interpret groundwater age from Khan et al. (2019) to\nexplain partial recharge of R1.\n Expand the significance of this study in using MT method to find deep freshwater resources and\nthe widespread presence of paleo-valleys during the sea-level lowstand.\n Modify the title, abstract, subtitles, and discussion due to the change in R1 interpretation.\nREVIEWER COMMENTS\nReviewer #2 (Remarks to the Author):\nI thank the authors for their revisions. My comments have been largely addressed. An aspect that is not clear (and which relates to my former review suggestions regarding 14C) can still be found in lines 150 to 152, where the authors state that „Although uncorrected 14C ages are reported instead of corrected ages, which decrease the uncorrected values by 1.0 14C kyr, due to the high uncertainty associated with initial 14C content corrections, ..“ This reads as if all uncorrected ages can be corrected by 1.0 kyr, which does not appear to be the case. I assume the authors intended to state something else ( – and certainly not relate to all samples). And if the „1.0 kyr“ is meant to relate to specific samples/locations, some more information or context would be helpful (and it would be interesting to know why corrected 14C can be given in this case and not elsewhere in the study area, in case this is meant here).I would hence would suggest to rephrase this sentence, and also to cut this (5 lines) long sentence into separate sentences.\nReply: We thank reviewer 2 for addressing the remaining issues of our revised manuscript.\nThe uncorrected 14C ages from Khan et al. (2019):\nLine 139-151: We modified our interpretation from Khan et al. (2019) groundwater data. Given\nthe uncertainty of the uncorrected ages, we interpret the freshwater body R1 contains both\nHolocene and Pleistocene groundwater since the 14C age data within the interfluve cluster\naround 10-12 kyrs with two low-As wells exhibiting Pleistocene groundwater ages. Therefore,\nwe interpret R1 as a freshwater reserve containing paleowater with slow regional recharge\ninstead of primary fossil freshwater. We expanded the interpretation in this second revision\ncompared to the previous version but keep the main text to be succinct. Thus, readers can check\nthe original paper for the 14C age correction. We added the deep groundwater age samples from\nKhan et al. (2019) to Figure 3 in the Main Text for clarifying their location relative to our MT\nsurvey. I hope that the revised explanation addresses the remaining concern of reviewer 2.\nThe uniqueness of this study:\nIn the last section “Sea-level change and the formation of deep freshwater resources”, we\nemphasize the role of sea-level fluctuations and sedimentary processes in deep groundwater\nformation.\nLine 203-207: we added the significance of incised paleovalley controlling the freshwater\ndistribution. Incised paleovalley is a widespread phenomenon during sea-level lowstands in\nPleistocene. Paleo-river could erode soft, unconsolidated sedimentary layers and remove the\nimpermeable layer to confine deep groundwater system. We added three references associated\nwith paleovalley (refs 51, 60, 61). We also clarified that paleovalley can act as conduits for deep\ngroundwater salinization.\n51. Larsen, F., Tran, L. V., Van Hoang, H., Tran, L. T., Christiansen, A. V., & Pham, N. Q. Groundwater\nsalinity influenced by Holocene seawater trapped in incised valleys in the Red River delta plain. Nature\nGeoscience, 10(5), 376-381. (2017).\n60. Meyer, R., Engesgaard, P., & Sonnenborg, T. O. Origin and dynamics of saltwater intrusion in a\nregional aquifer: Combining 3‐D saltwater modeling with geophysical and geochemical data. Water\nResources Research, 55(3), 1792-1813. (2019).\n61. Ta, T. K. O. et al. Latest Pleistocene to Holocene stratigraphic record and evolution of the Paleo-\nMekong incised valley, Vietnam. Marine Geology, 433, 106406. (2021).\nLine 213-219: We emphasized the significance of using MT method to capture the full dimensions of\ndeep groundwater for volumetric estimation, compared to other EM techniques.\nWe double-checked the spelling and typos in the second revised manuscript.\nPoint-by-point response to referees’ comments – Second Revision\nReviewer #3 (Remarks to the Author):\nThis is a review of the revised manuscript titled \"Deep fossil freshwater in coastal Bangladesh\" by Huy\nLe and others. The revised manuscript largely addressed comments I had and those of the additional\nreviewers. I appreciate the authors' improvements to the manuscript, in particular the \"Fossil Fresh\nAquifers\" section. This section provides better background for someone not familiar with the impressive\nbody of work on the GBD and more cohesive support for the authors' interpretations (despite the few\nremaining concerns about the emplacement mechanism and age of R1 discussed below).\nI reread and reviewed the entire manuscript for technical content but stopped commenting on editorial\nissues after the first 50 lines or so. There were several instances of either awkward phrasing or missing\nwords. An editorial review prior to resubmission is recommended. The following items were noted for\nconsideration by the authors.\nL22\nconsider adding \"the\" before \"location\"\nL34-35\nconsider \"...deep groundwater at a regional scale...\"\nL35\nconsider changing \"sustainability\" to \"susceptibility\"\nL69-70\nThe topic sentence on L64 limits important sea-level fluctuations to the Pleistocene\nL108-193\nI can accept R2 being relict Pleistocene water as it appears to be cut off from the modern flow system by\nthe upgradient paleochannel and is more saline than R1. I have more difficulty with the certainty\nexpressed in the age of R1 and its origin ascribed to a body of water \"formed\" during the Pleistocene.\nI think the issue is largely semantics.\nFirstly, the origin of the water. The language used in the text implies connate water, although that term is\nno longer used in the manuscript, e.g:\n\"freshwater bodies formed during the Last Glacial Maximum’s (LGM) sea-level lowstand and\nsubsequently sealed by a finer-grained, low-permeability layer\" [L108-110]\nPoint-by-point response to referees’ comments – Second Revision\n\"remnants of preserved Pleistocene groundwater\" [L180]\nIntentionally or not, the language implies the water in R1 originated in that location and has remained\nthere for 12,000+ years. Perhaps it has, however, that has not been demonstrated nor is it not supported by\nsubsequent statements that indicate R1 is the distal end of a larger upgradient flow system (L123-125)\nwith slow transport times (L123-125).\nSecondly, the timing of recharge of the groundwater in R1. The location of the 14C ages supporting the\nassertion of Pleistocene water is provided only as \"northwest of our profile\". Other than the fact those\nsamples are somewhere in the western portion of the 100,000 km2 GBD (L50), the reader has no\ninformation whether these data are relevant -- were they collected 1 km, 10 km, or 100 km away? There is\na significantly stronger case for transferability if they're nearby. Consider showing the locations on Figure\n3 if they are within the map extent. The transport time of \"thousands of years or more\" (L124-125] is\nconsistent with Holocene as well as Pleistocene recharge. The 14C data presented provide evidence for\nshallow-deep difference in age in another part of the GBD system having similar characteristics to the\nregion of R1, which is useful, but not conclusive given the lack of information about their location. The\nevidence presented for Pleistocene water in R1 is not sufficiently compelling to make a definitive\nassertion as to its age, particularly in light of the wide range and complex distribution of groundwater\nages presented by Kahn et al. (2019), from which the example 14C data were taken.\nConsider using the uncertainties surrounding the emplacement mechanism and timing of R1 in your favor\ninstead of trying to force certainty with imperfect data. For example, if R1 is the distal end of a long, slow\nflow system, it is potentially a more usable and renewable resource than connate Pleistocene water. If the\nage of R1 is mid-Holocene (for example) rather than Pleistocene, it could indicate a more\nuseful/renewable source. Since the groundwater system hosting the freshwater in R1 is relatively coarse-\ngrained (capped by the clayey paleosol), could you pump water from R1 without substantial drawdown?\nCould you potentially induce a bit more recharge from modern recharge zones or set up managed aquifer\nrecharge systems? The MT transect is a really useful first step to asking the next set of important\nquestions relating to development and management of the freshwater resource in the GBD\nReply: We thank reviewer 3 for the thorough feedback on our revised manuscript. For this new revised\nversion, we modified and added details to address the remaining concerns.\nL22\nconsider adding \"the\" before \"location\"\nL34-35\nconsider \"...deep groundwater at a regional scale...\"\nL35\nconsider changing \"sustainability\" to \"susceptibility\"\nL69-70\nThe topic sentence on L64 limits important sea-level fluctuations to the Pleistocene\nWe changed the terms based on these recommendations.\nPoint-by-point response to referees’ comments – Second Revision\nLine 108-193:\nAgain, we thank reviewer 3 for the suggestions to improve our interpretation. See “Freshwater reserves”\nsection for more details. Here, we summarize the main change based on your feedback:\nWe added the groundwater age samples from Khan et al. (2019) to Figure 3 in our manuscript to show the\nlocation of their data relative to our MT survey.\nWe agree that the groundwater ages show both Holocene and Pleistocene ages. Therefore, we modified\nour interpretation for R1 from primary fossil groundwater to a freshwater reserve containing both\nHolocene and older groundwater due to regional recharge.\nWe reinterpreted Khan et al. (2019) data in lines 139-151. Given the uncertainties about the timing of R1\nand other potential emplacement mechanism, we propose that R1 is likely to form during the LGM sea-\nlevel lowstand with slow regional recharge during the Holocene.\nWe also checked words and contexts that imply fossil freshwater in the text to remove the ambiguity of\ntiming since MT only provides mainly spatial constraints, and we used other studies to put temporal\nconstraints on our results.\nWe also proposed two hypotheses of R2 formation. It could be a relict paleo-groundwater due to R2’s\nisolation from R1 and northern recharge. An alternative explanation is that R2, like R1, is a freshwater\nreserve with slow regional recharge from West Bengal. However, we had limited constraints on R2.\nWe also checked typos and weird phasing in our main text."}, {"round": 3, "reviewer_comments": {"reviewer_2": "I thank the authors for their revisions. My comments have been largely addressed. An aspect that is not clear (and which relates to my former review suggestions regarding 14C) can still be found in lines 150 to 152, where the authors state that „Although uncorrected 14C ages are reported instead of corrected ages, which decrease the uncorrected values by 1.0 14C kyr, due to the high uncertainty associated with initial 14C content corrections, ..“ This reads as if all uncorrected ages can be corrected by 1.0 kyr, which does not appear to be the case. I assume the authors intended to state something else ( – and certainly not relate to all samples). And if the „1.0 kyr“ is meant to relate to specific samples/locations, some more information or context would be helpful (and it would be interesting to know why corrected 14C can be given in this case and not elsewhere in the study area, in case this is meant here).I would hence would suggest to rephrase this sentence, and also to cut this (5 lines) long sentence into separate sentences.\nFurthermore, I would suggest to clarify line 163 where „.. the LGMP surface is not mapped“, whether it does not extend as far, or if it could not be mapped (for technical reasons), as in the current form the relation to „..reducing ist isolation..“ (L162) does not appear straightforward.", "reviewer_3": "The authors have provided an excellent revision to their manuscript. The arguments and conclusions are considerably stronger and make for a compelling interpretation of their excellent MT transect. I offer only a handful of very minor suggestions for clarification. Best wishes for successful publication!\nL81 A reference to Figure 2b would be helpful after \"...southward dip in its upper boundary.\" I found myself looking at figure 2a and trying to understand how the authors arrived at an aquifer depth of 600 m, which is clearly shown in Figure 2b. Alternatively, if it is desired to stay with the resistivity measurements, consider a different depth reference that coincides with one of the contour intervals on Figure 2a.\nL113-117 It seems as if two different mechanisms are provided to explain the capping of the freshwater sand units. The first sentence invokes rapid deposition during transgression, while the second sentence invokes a paleosol. Please clarify.\nL554 No location provided for data"}, "author_response": "Point-by-point response to referees’ comments – Final round of revision\nDate: October 2, 2025\nManuscript Number: NCOMMS-24-67705B\nTitle of Article: Buried deep freshwater reserves beneath salinity-stressed coastal Bangladesh\nName of the Corresponding Author: Huy Le\nEmail Address of the Corresponding Author: hdl2115@columbia.edu\nDear All,\nWe thank two reviewers for their approval of the previous version of our manuscript and suggestions to\ncomplete the final draft. We have incorporated these changes to the main text. Also, we provide the\nreviewer comments below in black text, and our replies in blue.\nBest regards,\nHuy Le and coauthors\nREVIEWER COMMENTS\nReviewer #2 (Remarks to the Author):\nI thank the authors for their revisions. My comments have been largely addressed. An aspect that is not\nclear (and which relates to my former review suggestions regarding 14C) can still be found in lines 150 to\n152, where the authors state that „Although uncorrected 14C ages are reported instead of corrected ages,\nwhich decrease the uncorrected values by 1.0 14C kyr, due to the high uncertainty associated with initial\n14C content corrections, ..“ This reads as if all uncorrected ages can be corrected by 1.0 kyr, which does\nnot appear to be the case. I assume the authors intended to state something else ( – and certainly not relate\nto all samples). And if the „1.0 kyr“ is meant to relate to specific samples/locations, some more\ninformation or context would be helpful (and it would be interesting to know why corrected 14C can be\ngiven in this case and not elsewhere in the study area, in case this is meant here).I would hence would\nsuggest to rephrase this sentence, and also to cut this (5 lines) long sentence into separate sentences.\nReply: We modified the sentence to clarify the point from lines 150 to 152.\nSince the original work from Khan et al. (2019) only reported the uncorrected ages for all samples in their\nmain paper, they reported the corrections in their supplementary materials and explained why the\ncorrection values were not used. However, the high uncertainty in initial 14C content corrections does not\nchange the primary findings in their paper that low-As groundwater are associated with deep aquifers\nwith Pleistocene in age. Therefore, we clarified in the previous version that we used the uncorrected\nvalues to interpret the groundwater ages and added details to explain why using these values is still\napplicable. We summarized their work below:\nPoint-by-point response to referees’ comments – Final round of revision\n“Although uncorrected 14C ages are reported instead of corrected ages, which decrease the\nuncorrected values by 1.0 14C kyr, due to the high uncertainty associated with initial 14C content\ncorrections13, the groundwater 14C data exhibiting both Holocene and Pleistocene ages…”\nIn the modified version, we cut off the confusing part “which decrease the uncorrected values by\n1.0 14C kyr” and separated the long sentence into two sentences. We focused on the main\nobservations and interpretation here and not delved into details about the corrections or\ngroundwater sample analysis in Khan et al. (2019).\nFixed:\nLine 148-152:\n“We reported the published uncorrected 14C ages instead of corrected ages because of the high\nuncertainty in correcting for the initial 14C content13, however, these 14C ages delineate a clear spatial pattern\nin the distribution of deep groundwater ages. The older groundwater in interfluve areas suggests that deep\nconfined systems protected by the LGMP, like the R1 reservoir, are likely to preserve paleowater at least\nfrom latest Pleistocene with limited slow regional recharge during the Holocene.”\nFurthermore, I would suggest to clarify line 163 where “... the LGMP surface is not mapped“, whether it\ndoes not extend as far, or if it could not be mapped (for technical reasons), as in the current form the\nrelation to “..reducing ist isolation..“ (L162) does not appear straightforward.\nReply: Thank you for your suggestion. We clarified the southernmost extension of the LGMP.\nLine 163-164:\n“Unlike R1, the LGMP surface is not mapped as far south as R2, reducing its isolation from saline\nHolocene groundwater and regional mixing.”\nFixed:\nLine 161-164:\n“The LGMP surface may end or become deeper than the maximum depth of tube wells at the southern\nlimit of data44. In the Sundarban Mangrove Forest, the LGMP cannot be mapped due to the lack of wells.\nIf the LGMP did not develop as far south as R2, it would reduce R2’s isolation from saline Holocene\ngroundwater and regional mixing.”\nReviewer #3 (Remarks to the Author):\nThe authors have provided an excellent revision to their manuscript. The arguments and conclusions are\nconsiderably stronger and make for a compelling interpretation of their excellent MT transect. I offer\nPoint-by-point response to referees’ comments – Final round of revision\na handful of very minor suggestions for clarification. Best wishes for successful publication!\nL81 A reference to Figure 2b would be helpful after \"...southward dip in its upper boundary.\" I found\nmyself looking at figure 2a and trying to understand how the authors arrived at an aquifer depth of 600 m,\nwhich is clearly shown in Figure 2b. Alternatively, if it is desired to stay with the resistivity\nmeasurements, consider a different depth reference that coincides with one of the contour intervals on\nFigure 2a.\nReply: Thank you for your clarification.\nWe agree that Figure 2b shows the depth extension of R1 up to 600 m. Initially, we used 3 Ωm contour\nto describe the resistivity anomalies, R1 and R2. After doing salinity conversion, we followed the\n3 psu contour in salinity (Fig. 2b) to interpret depth extension of R1 in the main text. These\nexplains why we used the “600-m” for depth of R1. To keep the description of dimensions of R1\nand R2 consistent with resistivity contours in Fig. 2a, we modified the “600-m depth” to “800-m\ndepth” on line 78.\nL113-117 It seems as if two different mechanisms are provided to explain the capping of the freshwater\nsand units. The first sentence invokes rapid deposition during transgression, while the second sentence\ninvokes a paleosol. Please clarify.\nReply: We clarify the mechanism protecting deeper freshwater from vertical mixing with shallow saline\ngroundwater. The LGMP formed during the glacial period is the primary factor protecting the freshwater\nlayer.\nLine 113-117:\n“A subsequent rapid sea-level rise flooded the coastal GBD45,48, trapping freshwater beneath fine-grained\ntransgressive sediments31,44,48. The Last Glacial Maximum Paleosol (LGMP)44, a stiff clay layer formed\nby weathering of the Pleistocene interfluves exposed during glacial periods (Ext. Fig. 7), acts as a barrier\ninhibiting vertical mixing between dense saltwater and lighter freshwater below44,48.”\nFixed:\nLine 111-115:\nWeathering of the Pleistocene interfluves exposed during glacial periods (Ext. Fig. 7) created a stiff clay\nlayer known as the Last Glacial Maximum Paleosol (LGMP)44. The subsequent rapid sea-level rise\nflooded the coastal GBD45,48, trapping freshwater beneath fine-grained transgressive sediments and the\nLGMP31,44,48. The LGMP acts as a barrier inhibiting vertical mixing between dense saltwater and lighter\nfreshwater below44,48.\nL554 No location provided for data\nReply: Thank you for pointing that out. We published the data on Zenodo.org website at\nhttps://zenodo.org/records/17353663. We clarify the data location in Data Availability Section."}]

Decision letter