Geophysical Delineation of Freshwater–Saline Water Interfaces in Coastal Area of Southwest Bangladesh

Insufficiency of potable water is acute in southwest (SW) coastal areas of Bangladesh. The local population ignores the depth to saltwater/freshwater interface causing many unsuccessful waters wells drilling. In this context, a combined use of borehole logs, geophysical well logs, vertical electrical soundings (VES), electrical resistivity tomography (ERT) and electrical conductivity (EC) of sampled waters was performed to identify saltwater/freshwater interface depths in this coastal part. The study shows that the depth to freshwater/saltwater interface varies from place to place occurring commonly between 190 to 285 m, and locally as shallow as 146 m. The shallow freshwater/saltwater interface depth is greatly influenced by the upconing of fresh water from the deep aquifer (DA) near the major rivers and coast compared to the landward part and is mixed with more saline waters above. Vertically infiltrated saltwater is the main cause of brackish water in the upper shallow aquifer (USA), which is hydraulically connected with the lower shallow aquifer (LSA), and not directly recharged from the Bay of Bengal in the south. The study will contribute to find out the depth of the potential freshwater aquifer and assess the aquifer vulnerability in the coastal area of SW Bangladesh.


Introduction
The coastal areas of Bangladesh covering an area of 47,201 km 2 accommodate approximately 46 million people. Knowing the depth at which fresh water occurs is of crucial importance to the people living in the saltwater-affected coastal area of SW Bangladesh. Different degrees of salinity (4000-25,000 µS/cm) affect more than 50% of the coastal areas of Bangladesh [1]. A total of 73% of the population in Bangladesh live in rural areas, and a majority of the rural people use hand-pump tubewell water for their drinking-water supply [2]. The diameter of the tubewell is normally 3.81 to 6.35 cm. Groundwater is abstracted from a tubewell by using a hand pump, as well as a centrifugal pump. The development of shallow or deep tubewells is sometimes limited because of the inaccessibility of aquifers [3,4]. About 61% of the population in the coastal area have serious health issues due to potable water shortage [5]. Mostly women and young girls suffer severely, as they have to collect water for domestic purposes. Hypertension, heart, and kidney disorders are common among them, as they drink insufficient water [6,7]. In some places, women and girls have to walk to a water source about five km away to collect fresh water for domestic purposes. There are some areas along the coast of Bangladesh, where water in both shallow hand-pump tubewells (SHTWs) and deep tubewells (DTWs) is not suitable for consumption due to excessive salinity. Tubewell failure is also a very common phenomenon in the coastal area, due to the unknown freshwater/saltwater interface depth in the aquifer. In many settlements, people store rainwater in manmade or natural ditches that serve as the only source of domestic water [3,8]. Sea-level change also affects the groundwater salinity and/or migration of salt water into groundwater aquifers in Bengal Basin [9]. That is another reason why delineation of the freshwater aquifer boundary in coastal areas is of particular importance is to establish a baseline. Extreme climatic events, sea-level fluctuations, and increased groundwater abstraction from coastal aquifers disrupt the natural groundwater equilibrium condition, resulting in shifting of the freshwater/saltwater interface inland, toward the abstraction well, and seawater upconing and/or moving locally upward [10][11][12][13][14].
The coastal aquifer's water is inferior in terms of quality and often beyond the potable limits. Most of the tubewells constructed in this area are saline [15,16]. Seawater intrusion in the freshwater zone and impact of saltwater upon public health and agriculture have been investigated by some researchers [17][18][19][20][21][22], but the exact depth of salt/freshwater boundary has yet to be resolved. This interface depth can be detected by downhole geophysical well logs, for example, electrical resistivity logging, and water sampling or by indirect surface methods, such as VES or ERT [23]. Both direct and indirect methods have been used successfully to determine aquifer characteristics and monitor saltwater intrusion in groundwater research because of the impact of salinity on EC [24][25][26][27][28][29][30]. There is almost no work done yet to find out the actual depth of the saltwater-freshwater interface in the coastal areas of Bangladesh. In this study, geophysical well logging, VES, ERT, and water chemistry were integrated to locate this interface in this area. This will certainly be helpful in planning and implementation of safe and cost-effective drinking-water facilities in the SW coast of Bangladesh.

Hydrology and Climate
The study area is situated in the lower part of the Ganges River tidal floodplain in the southwestern part of the Bengal basin bordered by the Bay of Bengal in the south ( Figure 1). It comprises of two coastal districts: Barguna and Patuakhali. The rivers flowing in the area include Payra in the north, Bishkhali and Baleshwar in the west, Buriswar in the middle, and Tetulia in the east. Other small rivers and channels include the Andharmanik, Galachipa, Kajol, Bhairob, Kopotakkho, Palla, and Rabnabad. Furthermore, the Ganges River flows to the northwest-northeast of the study area. Mostly dendritic pattern distributaries of the major rivers, and irregularly distributed ponds and beels (a lake like wetland) are present throughout the study area. In total 200,984 small and 494 big ponds are situated in this area [31]. Some of the ponds are being used as the sources of domestic water supply, especially during the dry season. All the rivers and channels are influenced by the low tide (1 m) and high tide (3.1 m) with respect to mean lower low water (MLLW) each day in the SW part [32]. Tropical monsoon is dominant in Bangladesh, with a huge amount of rainfall occurring from June to October and moderate-to-high temperature and humidity. The annual average temperature is 25 • C. Relative humidity is around 60% in the dry season and 98% during the monsoon. The average annual rainfall for the period 2000-2014 was 2570 mm [22]. Water 2021, 13, x FOR PEER REVIEW 3 of 24

Geology and Hydrogeology
The study area, a part of the Bengal basin, has relatively flat topography with very gentle slope. The surface geology is mainly formed by the tidal deltaic and marshy peat and clay deposits. The Recent sediments primarily deposited by the Ganges form the area with a complex interfingering of very fine to coarse sediments from several regressions and transgressions in the past.

Geology and Hydrogeology
The study area, a part of the Bengal basin, has relatively flat topography with very gentle slope. The surface geology is mainly formed by the tidal deltaic and marshy peat and clay deposits. The Recent sediments primarily deposited by the Ganges form the area with a complex interfingering of very fine to coarse sediments from several regressions and transgressions in the past.

Geophysical Prospection
The subsurface information for depth, geometry, and characteristics of a source aquifer can be determined by using geophysical surveys [40,41]. Studies related to saltwater investigations frequently use geophysical resistivity surveys because of the direct sensitivity to the water salinity [42][43][44][45].

Geophysical Well Logs
The natural gamma ray (GR) log measures the natural emission of total gamma radiation from a geological formation in a borehole [46]. The emission depends on the mineral composition of the rock or sediment. Sands (quartz) tend to be less radioactive than clays and fine-grained sediments. GR log is most frequently used for lithologic and stratigraphic correlation. Results of the natural GR logs are reported in counts per second (CPS), commonly used by mineral industry, or API units if calibrated [46,47].
Two electrodes typically spaced 16 or 64 inches apart during electrical resistivity logging in the borehole, called short normal (SN) and long normal (LN) resistivity log, respectively, are used in combination with two electrodes located at the surface. They measure the apparent resistivity in ohmmeters (Ω.m) and are used to detect water quality (mainly the salinity) and lithology (mainly the clay content). The SN log is more strongly influenced by the invaded zone. It is useful for detecting thin beds (>16 inch) and sometimes to measure formation porosity. The LN log measures an intermediate resistivity, which may permit calculation of both invaded zone resistivity ( ) and true rock resistivity ( ). GR and LN logs were combined to deduce the aquifer properties up to the depth of 380 m below the ground level (bgl). In this study, 14 natural GR and LN resistivity logs were used. The combined use of short and long normal may be very helpful in identifying intruded zones. However, SN is more influenced by the drilling fluid in the invaded zone, which may misinterpret the results. When SN shows higher resistivity than LN, this indicates a high resistivity mud has been used (Rmf > Rw). We have carefully checked both the SN and LN. We observed that SN has lower resistivity compared to LN in our studied

Geophysical Prospection
The subsurface information for depth, geometry, and characteristics of a source aquifer can be determined by using geophysical surveys [40,41]. Studies related to saltwater investigations frequently use geophysical resistivity surveys because of the direct sensitivity to the water salinity [42][43][44][45].

Geophysical Well Logs
The natural gamma ray (GR) log measures the natural emission of total gamma radiation from a geological formation in a borehole [46]. The emission depends on the mineral composition of the rock or sediment. Sands (quartz) tend to be less radioactive than clays and fine-grained sediments. GR log is most frequently used for lithologic and stratigraphic correlation. Results of the natural GR logs are reported in counts per second (CPS), commonly used by mineral industry, or API units if calibrated [46,47].
Two electrodes typically spaced 16 or 64 inches apart during electrical resistivity logging in the borehole, called short normal (SN) and long normal (LN) resistivity log, respectively, are used in combination with two electrodes located at the surface. They measure the apparent resistivity in ohmmeters (Ω.m) and are used to detect water quality (mainly the salinity) and lithology (mainly the clay content). The SN log is more strongly influenced by the invaded zone. It is useful for detecting thin beds (>16 inch) and sometimes to measure formation porosity. The LN log measures an intermediate resistivity, which may permit calculation of both invaded zone resistivity (ρ i ) and true rock resistivity (ρ t ). GR and LN logs were combined to deduce the aquifer properties up to the depth of 380 m below the ground level (bgl). In this study, 14 natural GR and LN resistivity logs were used. The combined use of short and long normal may be very helpful in identifying intruded zones. However, SN is more influenced by the drilling fluid in the invaded zone, which may misinterpret the results. When SN shows higher resistivity than LN, this indicates a high resistivity mud has been used (Rmf > Rw). We have carefully checked both the SN and LN. We observed that SN has lower resistivity compared to LN in our studied resistivity logs.
Because it is more influenced by the resistivity farther away from the borehole (and thus less disturbed by borehole effects), LN usually provides a good approximation of formation resistivity (comprising matrix and pore water). We, therefore, usually emphasized LN only for this study. The SN was not used in this study, as it is strongly influenced by the invaded zone. These data were collected from the archive of DHPE, Dhaka, Bangladesh.
The qualitative interpretation of GR log was based on the identification of the sand (minimum cps) and shale (maximum cps) lines, in relatively constant zones of the log, from which a scale from 0 to 100% shale can be made. The shale line indicates a low permeability unit (clay or shale), and a deflection towards lower values indicates a decrease in the shale content and a higher permeability unit, such as sand and coarse grain sediments. The deflection in the LN log occurs towards left when ρw < ρmf and right when ρw > ρmf, where ρw and ρmf are formation water resistivity and mud-filtrate resistivity, respectively. At the same time, high resistivity indicates fresh water and the absence of clay/shale, whereas low resistivity identifies the presence of salt water/brackish water and/or shale/clay. Furthermore, the natural GR deflections were correlated with borehole log data to identify the different hydrostratigraphic boundaries. The deflections were also correlated with electrical resistivity data to identify fresh water and saline zones. For example, a saltwater sandy zone with higher total dissolved solids (TDSs) would be indicated by decrease in both natural gamma response and resistivity values, whereas freshwater-bearing aquifers with low TDSs would be delineated by increase in resistivity values and low gamma response. All the GR and LN resistivity logs were interpreted in the depth range from 50 to 380 m bgl (see Figure 3). The resistivity values for different lithological units were obtained from the LN resistivity log of the same borehole. The resistivity of pore water in a lithological unit was determined from the resistivity value of this unit by making use of Table 1. For example (Figure 3), when LN resistivity curve encompasses the sand and clay/clayey zones, the corresponding resistivity values are correlated with the formation resistivity (sand) and formation resistivity (clay) shown in Table 1, respectively. Following this correlation, we identified a series of water-resistivity values, such as B, MB, MF, and so on, from top to bottom of the log (see Figure 3 for an example). Saltwater/freshwater interface depth was then defined at the bottom of lithological unit containing saltwater and above the lithological unit containing fresh groundwater. The cross-sections were then constructed in different directions, as shown in Figure 1, to identify dominant litho-hydrologic units.

Resistivity Ground Surveys
Normally, a direct current (DC) is sent into the earth through a pair of current electrodes, and the resulting potential difference is then measured with the help of a pair of potential electrodes at the surface. The magnitude of the potential difference is closely related to the distribution of current in the subsurface. For a Schlumberger array configuration, where the current electrodes separation is much longer than the potential electrodes separation, Equation (1) is used to calculate apparent resistivity [48].
where a is the distance between the current electrode and the centre of the array, b the separation between potential electrodes, and R the resistance = ∆V I . When the electrode separation increases, the apparent resistivity is influenced by the resistivity of deeper layers. The Schlumberger array was used to conduct all the VES.
In the case of the Wenner configuration, the array spacing expands about the array midpoint, while maintaining an equivalent spacing between two adjacent electrodes, and Equation (2) is used to calculate apparent resistivity.
Water 2021, 13, 2527 where a is the electrode spacing.
Water 2021, 13, x FOR PEER REVIEW 7 of 24 , Figure 3. An example delineation of saltwater-freshwater interface in well A1.

Resistivity Ground Surveys
Normally, a direct current (DC) is sent into the earth through a pair of current electrodes, and the resulting potential difference is then measured with the help of a pair of potential electrodes at the surface. The magnitude of the potential difference is closely related to the distribution of current in the subsurface. For a Schlumberger array configuration, where the current electrodes separation is much longer than the potential electrodes separation, Equation (1) is used to calculate apparent resistivity [48].
where a is the distance between the current electrode and the centre of the array, b the separation between potential electrodes, and R the resistance = ∆ .
When the electrode separation increases, the apparent resistivity is influenced by the resistivity of deeper layers. The Schlumberger array was used to conduct all the VES.
In the case of the Wenner configuration, the array spacing expands about the array midpoint, while maintaining an equivalent spacing between two adjacent electrodes, and Equation (2) is used to calculate apparent resistivity.
where is the electrode spacing. In this study, VES 1 to VES 19 were conducted by a McOHM (Model-2115) OYO corporation, Japan, and taken from the MoST report, 2016 [49]. The rest of the VES and ERT were conducted by using a SYSCAL-R2 resistivity meter. This device is a high-power fully automated resistivity meter for DC electrical survey. The measuring process was adjusted to get the highest possible accuracy in the field conditions. The length of the current electrode spacing (AB/2) varies from 170 to 495 m, depending on the access to the survey area. ERT was carried out over a distance of 960 m, with 40 m unit electrode spacing. Soundings at 27 stations and one ERT were carried out to understand the groundwater quality and the geological structure of the area.
Obvious outliers were first removed from the measured resistivity data. The filtered data were then inverted in terms of 1-D multilayer earth models. The starting model was obtained from the available geological and hydrogeological information about the study area, followed by refinement of the initial model, using an iterative approach. Finally, the IPI2WIN-1D computer program [50] was used to obtain the final resistivity model. The inverse problem is solved by using the minimum number of layers required to fit the data (blocky inversion) and the regularized fitting minimizing algorithm, using Tikhonov's approach to solve ill-posed problems [51]. The results from the blocky inversion were considered to generate sharp contrast during inversion. A priori information about different layers depths and resistivities was used for regularizing the inversion and get more realistic solution [52].
• ERT Electrical imaging is a popular geophysical method in hydrogeology, given its sensitivity to crucial parameters, such as salinity and clay content [53]. ERT measurements were collected by repeating resistivity measurements with different electrode spacing at different locations. In contrast to VES, it can thus be used to characterize subsurface lateral resistivity variations, as well. For this study, a Wenner electrode configuration was used, and the apparent resistivity was calculated by Equation (2). The survey was carried out in the north-south direction perpendicular to the coast of the Bay of Bengal, starting approximately 50 m from the seashore (Figure 1). The Res2DInv software was used to invert the measured apparent resistivity data to create a model of the subsurface resistivity based on the iterative smoothness-constrained least-squares method [54,55].
The DOI can be defined by the depth below which the electrical resistivity of the subsurface no longer affects the measured resistivity data [56]. The sensitivity of ERT data decreases with increasing depth. The decreasing trend of resolution and sensitivity with depth is well documented [57]. Therefore, it should be addressed in calculation to avoid any misinterpretation or overestimation of subsurface resistivity structures. This is particularly more important in saline grounds. The electrical current flow tends to concentrate in highly conductive zones which will significantly decrease the DOI [58]. The DOI was calculated following the Equation (3) by Oldenburg and Li [56]. In the DOI method, two model inversions (ρ 1 and ρ 2 ) were carried out with different reference models (ρ re f 1 and ρ re f 2 ). The values of 0.1 and 10 times the average apparent resistivity values were taken for ρ re f 1 and ρ re f 2 , respectively. DOI = log 10 ρ 1 − log 10 ρ 2 log 10 ρ re f 1 − log 10 ρ re f 2 The discrepancy of the inverted model is insignificant when the value of DOI is approaching zero [56], and then the inverted model can be considered as a reliable model. In contrast, the DOI value around 1 indicates significant discrepancy of the inverted model. The DOI values between 0.1 and 0.2 are often used to delimit the reliable part of the inverted model [56][57][58][59][60]. In this study, we used the subjective value 0.1 to define the DOI of the inverted model.

Determination of Pore-Water Resistivity
We used the LN measured apparent resistivity, ρ LN , as an approximation of the formation resistivity, ρ LN ≈ ρ t . Archie's law [24] describes the bulk resistivity (ρ t ) as being related to pore-water resistivity (ρ w ) through the formation factor, F, as in Equation (4).
This relation works well in sandy formations. However, in shaly formations, the relationship between ρ t and ρ w is more complex, since the surface conductivity of clay particles must be accounted for. Patnode and Willie [61] prescribed an additional term to include the resistivity ρ m ( ·m) of the sediment/rock matrix.
In this study, for clayey deposits, we used a formation factor (F) equal to 2 and matrix resistivity of 50 Ωm, as deduced from detailed studies in the coastal plain of Belgium (e.g., Reference [62]). Table 1 shows the relation we used between groundwater resistivity and Water 2021, 13, 2527 9 of 23 formation resistivity for the different water-quality classes that were defined by De Moor and De Breuck [63]. We extended the table given by Walraevens et al. [64] for sands with F = 4 to include also the case for clayey deposits. Depth-wise lithological zones have been defined based on natural GR response and validated them with available nearest borehole. Shallow aquifers are affected by the presence of salt water. Both salinity and clay in the formation reduce the bulk resistivity and will lead to low resistivity values. The low resistivity values indicating saltwater intrusion might be equivocal. Nonetheless, it should be noted that the low resistivity value of clay may indicate freshwater condition, according to the Equation (5).
The resistivity values for different lithological units were obtained from the LN resistivity log of the same borehole. The resistivity of pore water was determined from the corresponding resistivity value of the same lithological units, according to Table 1. For example (Figure 3), when the LN resistivity curve encompasses the sand and clay/clayey zone, the corresponding resistivity values are correlated with the formation resistivity (sand) and formation resistivity (clay), as described in the Table 1, respectively. A similar technique has been applied for the interpretation of VES. Each VES having different zones of resistivity was checked carefully and correlated with the nearby borehole log. The depth at which the resistivity value of pore water obtained from VES becomes more than 6 Ω.m was taken as the beginning of freshwater depth. Usually, there is a steep change in resistivity at the interface between the VES layers, and 6 Ω.m is mostly a trustworthy value to differentiate them: layers saturated with fresh water will show higher resistivity, while layers with saline water will show lower resistivity [65]. The transition thus corresponds to a resistivity limit of 6 Ω.m and is referred to as the freshwater-saltwater interface depth, where the upper part is saline (<6 Ω.m) and the lower part is freshwater (>6 Ω.m) in the study area. The interface depth thus represents the depth at which fresh water can be found.
The detailed lithology of the sampled tubewells was studied. EC of water samples from the nested wells (a set of tubewells having different depths and drilled very near to each other) was also examined. EC of the sampled tubewells was assumed to be the same for the whole lithologic unit (well screen is placed in the lower part of the sand layer); the top of that lithologic unit was chosen as the starting depth of the fresh water. EC values from sampled tubewells were converted to resistivity values by considering the formation factor 4 (well screen in the sand layer). These resistivity values were then correlated with different water-quality types given in Table 1. The formations having pore-water resistivity more than 6 Ω.m in the abovementioned datasets were termed as freshwater zone, and formations with resistivity less than 6 Ω.m were termed as saline water zone.

Saltwater-Freshwater Interface Map
This is the first attempt to construct a saltwater and freshwater interface map of the study area by considering geophysical logs, geophysical survey, groundwater flow direction, water chemistry, and the borehole logs with chloride data. Firstly, 14 resistivity logs, 27 VES, 1 ERT, and water EC of the 36 selected tubewells were examined. Considering the resistivity values, local geology, geomorphology, and other influencing factors, such as population growth centers, we produced a map of the freshwater-saltwater interface depth. Firstly, all points with the interface depths were posted on a plain piece of paper. The contour lines were then drawn manually and scanned and finally digitized to present a version, using QGIS 3.16. The depth of the saltwater-freshwater interface was differentiated by a series of colors ranging from dark green to red, where the dark green color represents the shallowest and red color represents the largest depth to fresh water.

Geophysical Cross-Section
Example cross-sections are given in Figure 4A-C. In these sections, three distinctive layers of different salinity types are identified, following the hydrostratigraphy of the study area. In all cross-sections, the upper part is mostly brackish (B), with its actual depth levels varying from place to place. This zone extends on average up to 120 m, although it starts at around 100 m and ends at 200 m near the coast. A less saline (moderately brackish (MB)) zone up to the depth of around 245 m is present below the top zone. The thickness of this zone is higher in the P7/P8 section than that in the G2/G1 section. This zone is absent near the coast and reaches a maximum depth of 240 m in the Patuakhali town. A layer of weakly fresh (WF) water is present below 180 m depth in the Patuakhali town and at 200 m in the coast. The presence of a moderately brackish (MB) zone below 350 m depth indicates either saline pockets or a saline-water zone ( Figure 4B,C). There is not any evidence of direct seawater intrusion observed in any part of the study area. Some more exploratory DTWs below 350 m with geophysical logs and water sampling are needed, especially in the coastal area, to confirm this interpretation.
From these sections, it is also evident that seawater is not coming directly from the sea, but rather is leaching from the surface. Some recent studies [22,66,67] have shown that high salinity in the study area up to depth of ca. 200 m bgl is due to vertical infiltration of dissolved evaporite salts, rather than direct saltwater intrusion from the Bay of Bengal.   From these sections, it is also evident that seawater is not coming directly from the sea, but rather is leaching from the surface. Some recent studies [22,66,67] have shown that high salinity in the study area up to depth of ca. 200 m bgl is due to vertical infiltration of dissolved evaporite salts, rather than direct saltwater intrusion from the Bay of Bengal.

VES
The 1-D models obtained from VES show that the investigated subsurface comprises five or six geo-electric layers, showing different salinity zones, with a maximum investigation depth of less than 285 m. The topmost zone up to ca. 15 m depth from the VES located near the rivers has relatively high resistivity values in T1-3, P1-3, A1-2, VES 13, and VES 14. In groundwater-quality class, this corresponds to a WF zone. This can be the effect of freshwater invasion from the surrounding surface freshwater bodies. There are some layers of very low resistivity values underneath the topmost soil up to larger depth. The markedly low resistivity values of these layers could be attributed to the effect of seawater intrusion, but the main water table is mostly above the sea level and has a flowing trend from the recharge area in the north to the Bay of Bengal in the south [22]. The water-quality class of these layers varies from B to S in some cases. The resistivity sounding curves at different locations show different behaviours; the inversion results are also different. This may be due to the different levels of saltwater effects and lithological variation in the study area. It is evident from the model that resistivity is very low (0.98 to 1.68 Ω.m), up to 100 m depth. The presence of clay and silt layers inhibits a rapid downward movement of saltwater. Below this zone, there is a zone of pore-water resistivity values around 6 Ω.m at the bottom of the model, indicating the possible presence of fresh water below ~200 m depth.
The root-mean-square (rms) errors between the calculated and observed resistivity values range from 2.4% to 12%. Example models obtained at the measured points VES T1 and VES 19 are given in Figure 5, where the model response presents a good fit with

VES
The 1-D models obtained from VES show that the investigated subsurface comprises five or six geo-electric layers, showing different salinity zones, with a maximum investigation depth of less than 285 m. The topmost zone up to ca. 15 m depth from the VES located near the rivers has relatively high resistivity values in T1-3, P1-3, A1-2, VES 13, and VES 14. In groundwater-quality class, this corresponds to a WF zone. This can be the effect of freshwater invasion from the surrounding surface freshwater bodies. There are some layers of very low resistivity values underneath the topmost soil up to larger depth. The markedly low resistivity values of these layers could be attributed to the effect of seawater intrusion, but the main water table is mostly above the sea level and has a flowing trend from the recharge area in the north to the Bay of Bengal in the south [22]. The water-quality class of these layers varies from B to S in some cases. The resistivity sounding curves at different locations show different behaviours; the inversion results are also different. This may be due to the different levels of saltwater effects and lithological variation in the study area. It is evident from the model that resistivity is very low (0.98 to 1.68 Ω.m), up to 100 m depth. The presence of clay and silt layers inhibits a rapid downward movement of saltwater. Below this zone, there is a zone of pore-water resistivity values around 6 Ω.m at the bottom of the model, indicating the possible presence of fresh water below~200 m depth.
The root-mean-square (rms) errors between the calculated and observed resistivity values range from 2.4% to 12%. Example models obtained at the measured points VES T1 and VES 19 are given in Figure 5, where the model response presents a good fit with measured data and there are notable variations from zone to zone in resistivity values obtained by 1-D inversion. The 1-D models of the investigated VES reveal six geo-electrical layers. The model layers show an alternation of clayey and sandy materials from top to bottom. The soundings show decreased resistivity up to ca. 200 m bgl, followed by a rising resistivity trend. The deepest layer (layer six), having relatively high resistivity values, could be attributed to the freshwater-bearing formation. The VES T1 is located just a few meters from a river and shows a relatively high resistivity value (ca. 6 Ω.m), up to 17 m in depth. On the other hand, VES 19 is located ca. 500 m meters away from a river and shows low values, up to 22 m, compared to VES T1. The relatively high resistivity value (18.4 Ω.m) is observed at the former site below 146 m, while, at the latter site, below 205 m in depth.  Figure 6A shows the inverted resistivity model after four iterations, with a 2.6% final RMS error. The model shows that the site is characterized by relatively low background  Figure 6A shows the inverted resistivity model after four iterations, with a 2.6% final RMS error. The model shows that the site is characterized by relatively low background resistivity (1-6 Ω.m), depicting a general trend of low to high resistivity values from top to bottom. The low resistivity values in the profile infer the presence of high salinity. Very low (<1.75 Ω.m) but irregularly distributed resistivity is observed from the ground surface to 120 m depth. Extremely low resistivity (~1 Ω.m) observed locally in two places in the upper part of model is due to infiltration of saltwater. The presence of a pond in the middle of the profile is indicated by a relatively high resistivity value (>2 Ω.m). The gradual increase of resistivity with depth is observed below 120 m depth. The saline-water-filled pores of fine-grain-size sediments are responsible for the low resistivity values. The extent of saltwater movement is mainly controlled by the vertical infiltration. There is no evidence of lateral migration of salt water directly from the Bay of Bengal in the south of the profile.

Discussion
The potential of an aquifer for freshwater extraction depends on the hydraulic conductivity, storativity, transmissivity, and the quality of water. The varying lithology in the studied borehole logs includes alternating clay, clayey sand, silt, silty clay, very fine sand, fine sand, and medium sand. It is typical of the SW coast of Bangladesh [22]. The coastal plain sand below 200 m depth is the major source of fresh groundwater for domestic purposes in the study area.
The integrated study of natural gamma and resistivity logs, VES, ERT, and EC of sampled waters in this study revealed three hydrostratigraphic sections. The top section is USA, corresponding to low-resistivity brackish water to moderately salt water in some places. Its thickness is around 100 m on average, and it rests on a moderately brackish LSA zone. A few VESs near the rivers indicate the USA top reflected by high resistivity values. These localized high resistivity values may be due to the effect of dry soil, heavy monsoon downpour, and infiltration from the surrounding rivers, ponds, and lakes. This zone extends down to maximum 15 m in depth. However, extracting groundwater from this zone will not be suitable for drinking and domestic use and should be discouraged because of the potential risk of pollution due to poor sanitation and agricultural activities [68].
The second zone, extending between 100 and 220 m, represents a moderately brackish water zone. The transition between the high-salinity upper zone and relatively lowersalinity middle zone is gradual. The zone below ca. 200 m in depth has higher resistivity values, indicating weakly fresh water. Most of the pumping water wells are located in this zone to provide water for domestic purposes. The depth to this weakly freshwater zone varies from place to place, generally between 190 and 285 m, with a minimum depth of 146 m ( Figure 7A). The relatively shallow interface depth is observed mainly along the rivers and near the coast. The groundwater flow system of the study area is somewhat complex, but we find some kind of fresh groundwater upconing in the areas of major rivers. The zoomed in section of the southern part ( Figure 7B) shows a similar pattern of upconing of fresh groundwater from the DA, and that is mixed with some more saline groundwater above. The map ( Figure 7B) also clearly shows that the interface depth decreases near the main rivers and coast, whereas the depth relatively increases towards the continental part of the study area. The water-salinity-type profiles ( Figure 8A,B) based on VES show a general trend of decreasing salinity with depth. The salinity along the major rivers and coast is low compared to the landward part. This phenomenon also indicates the upconing of fresh groundwater from the DA and mixing with saltwater from the shallow aquifers above. We would recommend to further investigate the hypothesis with different scenarios, such as GW abstraction, sea-level rises due to climate change, etc. The use of radiocarbon ( 14 C), along with stable isotopes, could be very helpful to reveal the age and recharge mechanisms of the water in the study area, and we plan to perform this in the future.
The presence of some saltwater-contaminated ponds corresponds to low resistivity values at shallow depth. It is worth noting that there is no sign of saltwater movement from the sea to the landward side of the study area. Rather, saltwater is infiltrating from the top layers downwards. Extensive pumping from the deeper aquifers, especially in the towns and population growth centers, causes vertical infiltration of salt water from the upper aquifers and is also responsible for the salinization of the DA. The salinity in the top layer is the result of the leaching of evaporite deposited on the depressed/low-lying areas and mixing with marine-induced flood water [22,66]. Some authors [69,70] also showed that lateral seawater intrusion is slow compared to the downward infiltration of saltwater as a natural phenomenon in the coastal area of SW Bangladesh. This could also be a cause of variation of the depth to saltwater-freshwater boundary in the study area. Groundwater flows from recharge areas having high hydraulic heads to low hydraulic heads in the coastal area. Low-density fresh water tends to remain above high-density salt water, and this probably can explain the presence of some freshwater zone due to accumulation of rainwater at the top of the upper shallow unconfined aquifers, where water is mostly brackish in some studied VESs. Many researchers [71][72][73][74][75][76] have documented repeated sea-level changes that influence and change the direction of groundwater flow in the coastal area of Bangladesh. The presence of saline water below 300 m in depth in some studied resistivity logs in the SW corner of the study area results from the Holocene marine transgression between 10,000 and 6000 BP, when the sea level increased fast, and the flow situation may have changed severely. At this time, seawater was trapped in the back-barrier dunes and lagoon environments.  Figure 7A for clarity of displaying trend of interface depth in the study area.  Figure 7A for clarity of displaying trend of interface depth in the study area.  The interface depth was plotted against the distance to the important rivers ( Figure 9). This graph illustrates that the interface depth drops with increasing distance to the river, indicating that the salinity depth distribution is influenced locally by the upward flow of fresh groundwater from the deep aquifer. A noticeable upward flow of fresh groundwater from the deep aquifer is marked vertically up to 150 m and laterally within 1500 m distance from the rivers. Beyond 1500 m from the river, the interface depth remains practically constant (Figure 9). The presence of fresh groundwater at a shallow depth in some studied VESs in the northern part of the study area near the river indicates freshening by river water and rainwater in the upper part of the shallow water table aquifer. Very similar results were also observed by Sarker et al. [22], using chloride profiles, where the depths of groundwater samples were 30 m to ca. 100 m, acquired perpendicular to the rivers in six different locations over the study area.
heads in the coastal area. Low-density fresh water tends to remain above high-density salt water, and this probably can explain the presence of some freshwater zone due to accumulation of rainwater at the top of the upper shallow unconfined aquifers, where water is mostly brackish in some studied VESs. Many researchers [71][72][73][74][75][76] have documented repeated sea-level changes that influence and change the direction of groundwater flow in the coastal area of Bangladesh. The presence of saline water below 300 m in depth in some studied resistivity logs in the SW corner of the study area results from the Holocene marine transgression between 10,000 and 6000 BP, when the sea level increased fast, and the flow situation may have changed severely. At this time, seawater was trapped in the backbarrier dunes and lagoon environments.
The interface depth was plotted against the distance to the important rivers ( Figure  9). This graph illustrates that the interface depth drops with increasing distance to the river, indicating that the salinity depth distribution is influenced locally by the upward flow of fresh groundwater from the deep aquifer. A noticeable upward flow of fresh groundwater from the deep aquifer is marked vertically up to 150 m and laterally within 1500 m distance from the rivers. Beyond 1500 m from the river, the interface depth remains practically constant (Figure 9). The presence of fresh groundwater at a shallow depth in some studied VESs in the northern part of the study area near the river indicates freshening by river water and rainwater in the upper part of the shallow water table aquifer. Very similar results were also observed by Sarker et al. [22], using chloride profiles, where the depths of groundwater samples were 30 m to ca. 100 m, acquired perpendicular to the rivers in six different locations over the study area.

Conclusions
Knowing the depth at which fresh water occurs is of crucial importance to the people living in the saltwater-affected coastal area of SW Bangladesh. Geophysical well logs (natural GR and LN electrical resistivity) are integrated with borehole logs, VES, and ERT from the SW coastal area of Bangladesh to infer the source and location of subsurface saltwater and freshwater zones. We were able to identify the existence of three major aquifers with different water qualities. The upper aquifer above 100 m bgl depth is mostly brackish and sometimes highly to moderately saline in nature and prone to surface contamination. High salinity in this aquifer is not due to lateral saltwater intrusion from the Bay of Bengal revealed by ERT and VES. The discontinuous aquitard below the USA allows saltwater movement to the hydraulically connected LSA at depths of 100-200 m; it bears mostly moderately brackish water. In this aquifer, the salinity level is a bit lower compared to the upper aquifer. A continuous clay/silty clay layer below the LSA inhibits saline-water migration to the DA encountered below 200 m bgl. Water quality in this aquifer is weakly fresh. The presence of salt water below 300 m in the SW corner of the study area revealed in some resistivity logs, VES, and EC of the sampled groundwater probably results from the influence of sea-level rise during the Holocene transgression. From the equivalent EC of water in the aquifer, estimated from the geophysical and hydrochemical data, we mapped the salinity in the aquifer qualitatively and determined the saltwater/freshwater interface depth. The maps show that the shallowest interface depth is observed along the major rivers and coast, due to upconing of fresh groundwater from the DA. A relatively higher depth of interface is observed in the landward part of the study area. The interface depth versus distance to the major rivers shows a positive correlation; that is, the interface depth is increasing with increasing distance to the river, indicating that the freshwater distribution is influenced by the upward flow of DA water. The salinity profiles indicate a general trend of decreasing salinity with depth, but, more importantly, near the major rivers compared to the more landward part. This also indicates the vertical flow of DA water, and that is mixed with more saline waters above. It is recommended that the future groundwater exploitation be done in the DA, 200 m bgl.
This study, therefore, indicates the importance of the combined use of geophysics and hydrochemistry in delineating accurate water-quality type in aquifers useful in well design, construction, development, and groundwater abstraction, especially in complex hydrogeological coastal area.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.