Geophysical Assessment of Seawater Intrusion into Coastal Aquifers of Bela Plain, Pakistan

Seawater intrusion is a major challenge in many coastal areas all around the world, mainly caused by over-exploitation of freshwater resources, climate change, and sea-level rise. Consequently, seawater intrusion reaches several kilometers inland, thus making the freshwater resources polluted and unsuitable for human use. Conventionally, the fresh-saline water interface is delineated by the number of laboratory tests obtained from boreholes. However, such tests suffer from efficiency in terms of data coverage, time, and cost. Hence, this work introduces Dar-Zarrouk (D-Z) parameters, namely transverse resistance (Tr), longitudinal conductance (Sc), and longitudinal resistivity (ρL) computed from non-invasive vertical electrical sounding (VES). Two-dimensional (2D) imaging of D-Z parameters provides a clear distinction of fresh-saline aquifers. Such techniques remove ambiguities in the resistivity interpretation caused by overlapping of fresh and saline aquifers during the process of suppression and equivalence. This study was carried out by 45 VES along five profiles in the coastal area of Bela Plain, Pakistan. D-Z parameters delineate fresh, brackish, and saline aquifers with a wide range of values such as freshwater with Tr > 2000 Ωm2, Sc < 3 mho, and ρL > 20 Ωm; saline water with Tr < 1000 Ωm2, Sc > 25 mho, and ρL < 5 Ωm; and brackish water with Tr between 1000–2000 Ωm2, Sc from 3 to 25 mho, and ρL between 5–20 Ωm. The D-Z results were validated by the physicochemical analysis using 13 water samples and local hydrogeological setting. The obtained results propose that D-Z parameters can be used as a powerful tool to demarcate the fresh-saline aquifer interface with more confidence than other traditional techniques. This geophysical approach can reduce the expensive number of borehole tests, and hence contributes to the future planning and development of freshwater resources in the coastal areas.


Introduction
In Pakistan, groundwater is an essential resource for sustaining the population and economy [1]. Recently, groundwater utilization in the coastal areas has increased to meet water demand for industrial, agricultural, domestic, and tourist uses [2]. Consequently, improper management, excessive development, and over-exploitation of freshwater have deteriorated both the quantity and quality of groundwater reserves [3]. The main concern is seawater intrusion, which has resulted in groundwater contamination and well reduction [4,5]. Seawater intrusion is a natural occurrence caused by the density difference between seawater and freshwater. Seawater intrusion occurs when, due to the over-extraction of groundwater resources, outflow exceeds the aquifer recharge which is dynamically connected to the sea [6,7]. Generally, problems arise, when the seawater enters the aquifer system due to excessive pumping beyond the natural recharge rate, or when the potentiometric surface lowered by the extreme pumping at local individual wells causes the upcoming of the natural boundary between saline water and freshwater [8,9]. In addition, coastal aquifers are sensitive to the sea level rise due to global warming, and the melting snowcaps add to the level of seawater encroachment and ease the accessibility of freshwater reserves [10]. Therefore, demarcation of the fresh-saline groundwater boundary is very important for extraction and management of freshwater reserves in the coastal areas.
Delineation of the saltwater intrusion has become a global issue within the last century. In the last two decades, many studies around the world were conducted to assess the fresh/saline aquifers [1,2,[10][11][12][13][14]. Classically, seawater intrusion is monitored using the physicochemical analysis of groundwater samples and piezometry through the measurements of water-level fluctuations [15,16]. For successful monitoring of saltwater encroachment, observation wells are required to be dispersed over the entire investigated area. However, the economic reason mainly causes the sparse density of test wells, and thus the distributed wells do not characterize the total area sufficiently [17]. With the intent to ease the lack of information and improve the monitoring of the aquifer system, an efficient technique is necessary. The surface geophysical methods such as vertical electrical sounding (VES) are known as one of the most proficient techniques to assess groundwater resources quantitatively and qualitatively [18][19][20][21]. Electrical resistivity is extensively used to delineate seawater intrusion into coastal aquifers. A clear difference in resistivity values of freshwater and saltwater is useful to assess the fresh-saline water interface in the coastal areas [22,23]. Low resistivity values show saltwater, whereas high values suggest freshwater. Subsurface resistivity depends on many factors such as water content, saturation degree, porosity, rock type, soil type, salinity, and mineralization of groundwater, etc. [12]. Conversely, some geophysical similarities cause uncertainty in the resolution. Sometimes, the interpretation shows an intermixing of resistivity values for fresh and saline groundwater [24]. Dar-Zarrouk (D-Z) parameters can remove the overlapping of fresh/saline aquifers, and provide a more certain, practical, and reliable solution to the interpretation [11,25]. Applications of the D-Z parameters, namely transverse resistance, longitudinal conductance, and longitudinal resistivity computed from VES measurements mark a clear boundary between saline and fresh groundwater. Many studies have used D-Z parameters for delineation of the fresh-saline aquifer interface [1,11,12,17,25]. However, these studies provide only a 1D (one-dimensional) delineation of fresh-saline aquifer. Mostly, ERT is used to delineate the 2D demarcation of fresh-saline groundwater, however, ERT cannot cover a large area. In order to get a clear insight into the fresh-saline aquifer, the 2D delineation of the aquifer system is essential. This work provides a new innovative approach for the 2D demarcation of fresh-saline aquifer using the D-Z parameters. Based on local geological settings and groundwater conditions, D-Z parameters were acquired for each distinct layer of fresh or saline water with various arrangements of resistivity and thickness. These parameters delineated saline and fresh aquifers with a wide range of specific values depending on local hydrogeological conditions and wells information. Freshwater was delineated by high values of transverse resistance and longitudinal resistivity but low values of longitudinal conductance, whereas saltwater was assessed with high values of longitudinal conductance but low values of the two other parameters.
During the last two decades, the coastal aquifers of Bela Plain have become polluted due to the over-extraction of groundwater resources mainly for agriculture demand. Groundwater, with seawater intrusion of more than 50 km from the Sonmiani Bay of Arabian Sea towards the inlands of the study area, has become unsuitable for drinking and agriculture purposes. The number of boreholes near the coast is discarded due to seawater intrusion. The main aim of the present investigation is to emphasize the application of D-Z parameters derived from VES data for 2D mapping of seawater intrusion into coastal areas of Bela Plain, Pakistan, and to display how this non-invasive geophysical technique can be used to reduce the number of expensive well samples of hydrogeochemical information. In addition, the monitoring boreholes data of physicochemical parameters were integrated with a geophysical interpretation to improve and validate the results of the D-Z parameters for the assessment of saline and fresh aquifers.

Study Area and Hydrology
Bela Plain is located in Lasbela District covering an area of 3100 km 2 . The study area is situated south of the Baluchistan Province between latitiude 25-26.55 • N and longitude 66-67 • E. Bela Plain is surrounded by hills of sedimentary rocks from three sides in the east, west, and north and Sonmiani Bay (Arabian Sea) in the south (Figure 1). The study area is located in the arid region having an average annual precipitation of 80-130 mm, high evaporation, and low recharge to groundwater. Agriculture is the main source of earning for people in the study area, which mainly relies on tubewells, streamflow, and rainfall. The investigated area has intermittent streams mostly, and the only perennial one is the Porali River in the southeast. The alluvial cover of unconsolidated deposits was formed during the quaternary age, and the errosion processes shaped the topographic relief of the study area.
The stream runoff and rainfall are the main sources of water in foothills and piedmont areas of Bela Plain. In northern areas, groundwater is mainly recharged by streamflows of Porali River through permeable beds. Whereas, rainafll is mainly the source of groundwater recharge for the rest of the areas in Bela Plain. In northern Mor and Par Ranges, Porali River at the Sinchi Bent guaging station, is the major source of water with an irrigation land of 4039 km 2 . The surface water hydrology, WASID (water and soil investigation division) is responsible for meausring the discharge rate of rivers and streams in the study area [26]. The mean annual flow is 120 million m 3 , which is mostly used in agriculture [27]. About 14-18 km downstream from the Sinchi Bent guaging station, Kud River joins the Porali River. At Mai Gondrani, the mean annual flow is 102 million m 3 [27]. Windar Nai with a drainage area of 1733 km 2 is a non-perennial stream, which direcly flows into the Arabian Sea in the south. Khantra Nai is another main non-perennial stream in the study area.
Gravel, sand, clay, and silt are the main lithologies in the investigated area as revealed by the drilling tests. However, gravel forms the most extensive and permeable aquifer, which lies in the unconsolidated deposits under the alluvial fans and piedmonts places. Generally, the study area has the unconfined aquifer system. However, the semi-permeable beds of silt and clay are dominant in the south and southwest parts where they make the confined aquifer sytem with saltwater. The water table varies between 0-15 m in the investigated area. The shape of the water table depth is of concave type in the southwest and middle of the investigated area along Porali River and Khantra Nai. The groundwater recharge rate is mainly controlled by the seasonal conditions of discharge from the surfacewater into the groundwater. The general flow direction of surface water in Bela Plain is from the northeast to the southwest towards the Arabian Sea. This investigation was conducted by the Water And Power Development Authority, Pakistan (WAPDA) under a project in order to assess seawater intrusion into the coastal area of Bela Plain for the development and management of freshwater resources. All the 45 VES were performed in the same time period/season (April/May, 2017) as the groundwater samples were collected from the 13 wells. The location of VES and wells inluding the hydrogeological setting of the investigated area is given in Figure 1a. Figure 1b shows a schematic diagram of seawater intrusion into the coatsal aquifer.    [27]); (b) a schematic diagram of seawater intrusion into the coatsal aquifer (modified after [18]).

VES Survey
Seawater intrusion into coastal aquifers is delineated by many techniques. It is never easy to directly interpret the measurements of apparent electrical conductivity quantitatively, since it depends on several properties, such as clay containing high concentrations of cations and dissolved electrolytes, temperature, porosity, etc. [28]. In addition, a high electrical conductivity shows saltwater, whereas its low values suggest freshwater [11]. This difference between the electrical conductivity of freshwater and saltwater is highly significant, and therefore, the VES method is more suitable to assess the fresh-saline groundwater interface in the coastal regions [29]. Geophysical methods such as vertical electrical sounding (VES) are widely used in environmental and hydrogeological studies given their ability to incorporate with other methods, such as hydrochemical and geological methods, through a comprehensive interpretation. A useful correlation is found between hydraulic and electrical properties of the saturated geological formations, which make the geoelectrical methods even more effective in such studies [30][31][32][33]. An integrated approach of the VES and geochemical method in combination with hydrogeological data has been used to evaluate the characteristics of groundwater and delineate the expansion of seawater intrusion into the coastal aquifer.
In this research, a VES survey was performed to compute the subsurface resistivity. Electrical resistivity depends on the characteristics of subsurface materials, such as rock texture, nature of mineralization, and electrical conductivity [34]. It varies from formation to formation. Coarse grains show high values of resistivity, whereas an increase in the clay content or salinity gives low values. In the VES survey, an electric current is injected into the subsurface through two electrodes known as current electrodes, and two other electrodes namely potential electrodes are used to measure the electrical potential. The measured voltage is converted into the resistance value using the resistivity meter. The anomalies are observed by plotting the resistance values on a log-log paper in the field. With each resistivity measurement, spacing between the current electrodes is increased to acquire the large depth of investigation. The signal to noise ratio is enhanced by increasing spacing between the inner potential electrodes. In this survey, 45 VES measurements ( Figure 1a) were acquired by a resistivity meter (Syscal R1 Plusmodel, Iris) using the Schlumberger array. The VES data points were divided into five profiles, i.e., A, B, C, D, and E as shown in Figure 1a. The spacing of current electrodes (AB) and potential electrodes (MN) was adjusted at all depths according to a relationship, i.e., 2(AB/2) > 5 MN [35]. The minimum and maximum half-current electrodes spacing (AB/2) was used as 1 and 200 m, respectively, whereas the potential electrodes spacing (MN) varied between 0.25-20 m. The resistivity measured directly from the VES data is apparent resistivity. This resistivity needs to be inverted, as it does not represent the true values of subsurface layers [36]. The apparent resistivity data were inverted by the IPI2WIN software version 3.1.2c, developed by Moscow State University [37]. The software program compared the theoretical curve with the inverted (true) resistivity and the calibrated field curve (Figure 2a). The shape of the modeled VES curves shows variations in different layers with different resistivities [35]. The VES models divide the subsurface formation into distinct layers with a specific value of thickness and resistivity. The correct interpretation of VES models is performed using the local information acquired from boreholes data and hydrogeological setting.

Dar-Zarrouk Parameters
Occasionally, the VES interpretation causes ambiguity in the resolution and resistivity values, which shows the overlapping of demarcation of saline and fresh groundwater [24]. D-Z parameters remove any such uncertainty in the resistivity interpretation, and provides a clear interface between the freshwater and saltwater. They present a more practical, definite, and reliable solution to the interpretation [11,25]. D-Z parameters mark a clear interface between saline and fresh groundwater with a wide range of their values without overlapping of values or aquifers. D-Z parameters, namely transverse resistance (T r ), longitudinal resistivity (ρ L ), and longitudinal conductance (S c ) were first introduced by Mailet in 1947. These parameters offer sufficient resolutions to subsurface resistivity for the delineation of seawater intrusion [11]. They efficiently mark a distinct boundary between saline and fresh groundwater with a specific range of their values [12]. These parameters are estimated through different arrangements of thickness and resistivity acquired from the VES models (Figure 2b), and thus they give a better understanding of the subsurface discrete layers for the evaluation of a boundary between saline and fresh groundwater [11,38,39]. Any geoelectrical layer in a VES model is defined by two fundamental properties, namely resistivity (ρ) and thickness (h) [40]. Various arrangements of resistivity (ρ) and thickness (h) for each lithologic unit of VES model describe the D-Z parameters [11,41]. Transverse resistance was estimated by the following arrangement of h and ρ [11,41]: where h and ρ are the aquifer thickness and resistivity of each layer in the VES model. T r is the transverse resistance measured in units of Ωm 2 . Longitudinal conductance was calculated by the following formula [11,41]: where S c is the longitudinal conductance of each saturated layer, and it is measured in Siemens (mho). The longitudinal resistivity is acquired using the following equation [11,41]: where ρ L is the longitudinal resistivity of each layer, and is estimated in Ωm. Equations (1)-(3) were used to estimate the D-Z parameters of each distinct layer of the VES model ( Figure 2b). D-Z parameters were estimated for all 45 VES models to get the 2D mapping of discrete saline and fresh groundwater in Bela Plain. Specific values ranges of T r , S c , and ρ L for saline, brackish, and fresh aquifers were obtained based on the hydrogeological setting and wells information of the investigated site. To obtain specific values ranges of the D-Z parameters, the lithological logs data of the 13 wells were integrated with the D-Z parameters of the selected 13 VES nearby the wells. An example of such integration between the selected VES 9 and the lithological well log 3 is shown in Figure 2a. However, these values ranges of D-Z parameters are not unique and may vary from area to area based on the local hydrogeological setting.

Geochemical Method
The physicochemical analysis was conducted for 13 groundwater samples obtained from the pumping wells with an average depth range of 60-110 m. The laboratory analysis of groundwater samples was conducted using the ion chromatography technique by the Pakistan Council of Research in Water Resource (PCRWR). The physicochemical study was performed by the recommended limits of the World Health Organization [43] for main anions such as chloride (Cl − ), sulphates (SO4 2− ), nitrate (NO3 − ), and bicarbonates (HCO3 − ), as well as main cations including magnesium (Mg 2+ ), sodium (Na + ), potassium (K + ), calcium (Ca 2+ ), and other geochemical parameters such as electrical conductivity (EC), pH, and total dissolved solids (TDS) [44]. The standard procedures method was used to determine the concentrations of calcium, bicarbonate, sodium, sulphate, potassium, chloride, and magnesium [44]. Bicarbonates and chlorides ions were obtained by the volumetric method. The UV-visible spectrophotometer was applied for the analysis of sulphates. The atomic absorption spectrometry was performed to analyze the main cations. The concentration of cations and anions was measured in milligrams per liter (mg/L). The concentration of the electrical conductivity in microsiemens per centimeter (μS/cm) was determined at 25 °C. The consistency of the groundwater quality parameters was confirmed using the ionic balance of water, and the acceptable reliability was found between +5% and −5%.

Interpretation of Resistivity Data
Resistivity measured in the field is not a true representative of subsurface geologic units, and therefore, apparent resistivity is inverted to true resistivity through an inversion procedure of the software [1]. Resistivity varies even with a minor change in lithology. The resistivity of sand is lower than the resistivity of gravel, similarly, the resistivity of silt/clay is lower than that of sand. Moreover, freshwater shows a high resistivity, while saline water indicates a low resistivity. Mostly, gravel/sand reveal freshwater, whereas clay delineates saltwater [12]. However, with a homogeneous and anisotropic subsurface, uncertainties can be expected in VES methods for the

Geochemical Method
The physicochemical analysis was conducted for 13 groundwater samples obtained from the pumping wells with an average depth range of 60-110 m. The laboratory analysis of groundwater samples was conducted using the ion chromatography technique by the Pakistan Council of Research in Water Resource (PCRWR). The physicochemical study was performed by the recommended limits of the World Health Organization [43] for main anions such as chloride (Cl − ), sulphates (SO 4 2− ), nitrate (NO 3 − ), and bicarbonates (HCO 3 − ), as well as main cations including magnesium (Mg 2+ ), sodium (Na + ), potassium (K + ), calcium (Ca 2+ ), and other geochemical parameters such as electrical conductivity (EC), pH, and total dissolved solids (TDS) [44]. The standard procedures method was used to determine the concentrations of calcium, bicarbonate, sodium, sulphate, potassium, chloride, and magnesium [44]. Bicarbonates and chlorides ions were obtained by the volumetric method. The UV-visible spectrophotometer was applied for the analysis of sulphates. The atomic absorption spectrometry was performed to analyze the main cations. The concentration of cations and anions was measured in milligrams per liter (mg/L). The concentration of the electrical conductivity in microsiemens per centimeter (µS/cm) was determined at 25 • C. The consistency of the groundwater quality parameters was confirmed using the ionic balance of water, and the acceptable reliability was found between +5% and −5%.

Interpretation of Resistivity Data
Resistivity measured in the field is not a true representative of subsurface geologic units, and therefore, apparent resistivity is inverted to true resistivity through an inversion procedure of the software [1]. Resistivity varies even with a minor change in lithology. The resistivity of sand is lower than the resistivity of gravel, similarly, the resistivity of silt/clay is lower than that of sand. Moreover, freshwater shows a high resistivity, while saline water indicates a low resistivity. Mostly, Water 2020, 12, 3408 8 of 23 gravel/sand reveal freshwater, whereas clay delineates saltwater [12]. However, with a homogeneous and anisotropic subsurface, uncertainties can be expected in VES methods for the characterization of subsurface lithologies [36]. Therefore, the VES models and local hydrogeological information obtained from wells are correlated to constrain the subsurface into discrete lithologic layers. In this study, the non-predictive bias evaluation of the subsurface was obtained by correlating 45 VES models with lithological logs of 13 boreholes (Table 1). A comparison between the interpreted VES 9 and lithological well log 3 is shown in Figure 2a. Similarly, the calibration (comparison) between the 45 VES models and 13 lithological well logs was performed to obtain a resistivity values range for each layer containing fresh-saline or brackish aquifer with specific lithology. This calibration delineated the subsurface into several lithological units, such as dry strata/topsoil cover (above water table) having resistivity >35 Ωm, clay (below water table) with saline water having resistivity <7 Ωm, clay-sand (below water table) with brackish water and resistivity between 4-24 Ωm, sand (below water table) having a freshwater and resistivity range of 18-42 Ωm, gravel-sand (below water table) with freshwater and resistivity from 34 to 60 Ωm, and gravel (below water table) having freshwater and resistivity >50 Ωm (Table 1). The calibration shows that the interpreted layers have very narrow resistivity ranges with an overlapping nature. The clay of saline water having resistivity <7 Ωm intermixes with the clay-sand of brackish water with resistivity between 4-24, and thus causes an overlapping of saline-brackish aquifer with resistivity between 4-7 Ωm. Similarly, the clay-sand of brackish water having a resistivity range of 4-24 Ωm interferes with the sand of freshwater with resistivity from 18 to 42 Ωm, and hence creates an overlapping of brackish-fresh aquifer with a resistivity range of 18-24 Ωm. In addition, Table 1 shows the intermixing of sand with gravel-sand for resistivity overlapping of 34-42 Ωm, and gravel with gravel-sand for resistivity overlapping of 50-60 Ωm. Therefore, such overlapping causes an ambiguity in resolution for the separation between saline and fresh groundwater. This implies that the VES technique cannot interpret the subsurface with a high resolution, especially when the lithological units have some similarities, such as clay-sand and clay or sand and clay-sand. However, D-Z parameters estimated from the VES models can remove any such ambiguity, and thus provide a more confident resolution to the interpretation of fresh-saline aquifers. These parameters have a wide range of values to mark a clear interface between freshwater and saline water.

Transverse Resistance
An important parameter of the D-Z parameters, namely transverse unit resistance (T r ), was estimated for each layer of the 45 VES models along profiles A-E through the arrangement (Equation (1)) of aquifer thickness and resistivity of the layer. Transverse unit resistance was used Water 2020, 12, 3408 9 of 23 to demarcate fresh, brackish, and saline water with a wide range of T r values based on the wells information and hydrogeological setting of the investigated area. High values of T r suggest that freshwater has sand or gravel as the main lithology, while its low values indicate saline water has clay as the main lithology. In this study, transverse unit resistance delineated saline water with T r < 1000 Ωm 2 containing clay, brackish water with T r between 1000-2000 Ωm 2 having clay-sand, and freshwater with T r > 2000 Ωm 2 having sand, gravel-sand, and gravel ( Table 2). The values range of T r for fresh, brackish, and saline water shows no intermixing, which suggests that the aquifers delineated by T r have no intermixing of saline-brackish or brackish-fresh aquifers. The distribution of transverse unit resistance values over 45 VES in Figure 3a suggests that T r increases from the sea (in the south) towards the inland (in the north). Based on the interpretation of T r , Figure 3b provides the demarcation of saline, brackish, and fresh aquifers at different depths (50, 100, 150, and 200 m). At a 50 m depth, T r varies between 50-4100 Ωm 2 with a T r range of 50-574 Ωm 2 for saline water, 1280-1705 Ωm 2 for brackish water, and 2520-4100 Ωm 2 for freshwater. At a 100 m depth, the T r values range is 10.9-36,400 Ωm 2 with T r between 10.9-574 Ωm 2 for saline aquifer, 1280-1705 Ωm 2 for brackish aquifer, and 3000-36,400 Ωm 2 for fresh aquifer. At 150 and 200 m depths, T r changes from 10.9 to 36,400 Ωm 2 having T r of 10.9-92.8 Ωm 2 for saline aquifer, 1280-1705 Ωm 2 for brackish aquifer, and 2730-36,400 Ωm 2 for fresh aquifer. Generally, T r (from 2520 to 36,400 Ωm 2 ) for freshwater increases with depth, whereas T r (from 574 to 10.9 Ωm 2 ) for saltwater decreases with depth, which suggests seawater intrusion with an increasing depth (Figure 3b). The 2D transverse unit resistance was obtained and interpreted for saline, brackish, and fresh groundwater along five profiles (A-E) (Figure 4a,b). Transverse unit resistance values vary along profile A to E with Tr between 11.3-28080, 8.8-6750, 11.5-9951, 11.3-9951, and 10.5-9951 Ωm 2 , respectively ( Figure 4a). Along profile A, the saline aquifer was delineated with Tr between 11.3-384 Ωm 2 for VES 11-13, brackish aquifer was revealed by Tr from 1386 to 1650 Ωm 2 for VES 8-10, and fresh aquifer was interpreted with a Tr range of 2700-28080 Ωm 2 for VES 1-9 (  The 2D transverse unit resistance was obtained and interpreted for saline, brackish, and fresh groundwater along five profiles (A-E) (Figure 4a,b). Transverse unit resistance values vary along profile A to E with T r between 11.3-28080, 8.8-6750, 11.5-9951, 11.3-9951, and 10.5-9951 Ωm 2 , respectively ( Figure 4a). Along profile A, the saline aquifer was delineated with T r between 11.3-384 Ωm 2 for VES 11-13, brackish aquifer was revealed by T r from 1386 to 1650 Ωm 2 for VES 8-10, and fresh aquifer was interpreted with a T r range of 2700-28080 Ωm 2 for VES 1-9 (Figure 4b). Seawater intrusion for profile A was observed about 45 km towards the inland along VES 8-13. Along profile B, freshwater was revealed by a T r range of 2788-6750 Ωm 2 for VES 14-16, brackish water was interpreted with T r of 1376 Ωm 2 for VES 16 and 17, and saline water was assessed with T r between 8.8-384 Ωm 2 for VES 11 and 17-22 (Figure 4b). Seawater intrusion of about 75 km was found along this profile for VES 11 and 16-22. Along profile C, saline water was evaluated with T r between 11.5-180 Ωm 2 for VES 20 and 26-29, brackish water was interpreted by T r of 1440 Ωm 2 for VES 26, and freshwater was delineated with T r between 3330-9951 Ωm 2 for VES 23-25 (Figure 4b). Seawater intrusion of 50 km was delineated along profile C with VES 20 and 26-29. Along profile D, freshwater was revealed by a T r range of 2275-9951 Ωm 2 for VES 23 and 30-33, brackish water was interpreted with T r between 1512-1705 Ωm 2 for VES 31-34, and saline water was delineated with T r between 11.3-390 Ωm 2 for VES 13, 34, and 35 (Figure 4b). This profile reveals seawater intrusion of about 27 km along VES 13 and 31-35. Along profile E, saline aquifer was revealed by T r between 10.5-560 Ωm 2 for VES 40-45, brackish aquifer was interpreted with T r 1650-1700 Ωm 2 for VES 38-42, and fresh aquifer was delineated with T r 2520-9951 Ωm 2 for VES 23 and 36-39 (Figure 4b). Seawater intrusion of about 65 km was revealed along profile E with VES 38-45. Figure 5a shows the integrated 2D maps of transverse unit resistance, and Figure 5b provides the demarcation of saline, brackish, and fresh groundwater for the integrated 2D T r maps. Generally, seawater intrusion of about 27-75 km is delineated from the Arabian Sea in the southwest towards the inland of Bela Plain in the northeast (Figure 5b). intrusion of about 65 km was revealed along profile E with VES 38-45. Figure 5a shows the integrated 2D maps of transverse unit resistance, and Figure 5b provides the demarcation of saline, brackish, and fresh groundwater for the integrated 2D Tr maps. Generally, seawater intrusion of about 27-75 km is delineated from the Arabian Sea in the southwest towards the inland of Bela Plain in the northeast (Figure 5b).

Longitudinal Conductance
Longitudinal unit conductance (S c ) was estimated using a combination (Equation (2)) of aquifer thickness and resistivity obtained from 45 VES models along profiles A-E. Based on the wells information and hydrogeological setting of the study area, a specific values range of longitudinal unit conductance was obtained to delineate fresh, brackish, and saline aquifers with a clear interface. Low values of S c reveal freshwater with gravel-sand, whereas high values delineate saline water mainly with clay. Saline aquifer was delineated with S c > 24 mho having clay, fresh aquifer was demarcated by S c < 3 mho with sand, gravel-sand, and gravel, and brackish aquifer was interpreted with S c between 3-24 mho containing clay-sand (Table 2). Figure 6a shows mapping of the longitudinal unit conductance, whereas the demarcation of saline, brackish, and fresh groundwater at different depths (i.e., 50, 100, 150, and 200 m) is given in Figure 6b mho for brackish aquifer, and 0.15-2.4 mho for fresh aquifer. At a 150 m depth, Sc varies between 0.23-1460 mho with a Sc range of 80-1460 mho revealing saline aquifer, 4.1-21.5 mho delineating brackish aquifer, and 0.23-2.2 mho interpreting fresh aquifer. At a 200 m depth, Sc ranges from 0.07 to 1460 mho with Sc between 42-1460 mho for saline aquifer, 4.1-20 mho for brackish aquifer, and 0.07-2.2 mho for fresh aquifer. Longitudinal unit conductance values for saline water increase with depth from 34.1 to 1460 mho and decrease for freshwater with depth from 2.4 to 0.07 mho, and thus suggest an increase in seawater intrusion with depth ( Figure 6b).  Figure 8a provides the integration of 2D S c maps. The fresh, brackish, and saline aquifers delineated by the integrated 2D S c maps over the entire study area are given in Figure 8b. The integration of 2D S c maps delineates seawater intrusion of about 35-90 km from the seaside in the southwest towards the investigated area in the northeast (Figure 8b).

Longitudinal Resistivity
Longitudinal resistivity (ρL) was estimated for 45 VES models along profile A-E from a relation (Equation (3)) obtained by the arrangement of aquifer thickness and longitudinal conductance. Longitudinal resistivity revealed saline, brackish, and fresh aquifers with a wide range of values based on local hydrogeological setting and wells information. Generally, a high ρL suggests freshwater mainly with gravel-sand, while a low ρL reveals saline water with clay as the dominant lithology. In this investigation, longitudinal resistivity delineated the saline aquifer with ρL < 5 Ωm having clay, the brackish aquifer with ρL from 5 to 20 Ωm containing clay-sand, and freshwater with ρL > 20 Ωm having sand, gravel-sand, and gravel ( Table 2). Mapping of the longitudinal resistivity over 45 VES is given in Figure 9a Figure 9 suggests that longitudinal resistivity increases (from 34 to 432 Ωm) with depth within the fresh aquifer and decreases (from 4.1 to 0.1 Ωm) with depth within the saline aquifer, which shows seawater intrusion from the seaside with the increasing depth (Figure 9b).

Longitudinal Resistivity
Longitudinal resistivity (ρ L ) was estimated for 45 VES models along profile A-E from a relation (Equation (3)) obtained by the arrangement of aquifer thickness and longitudinal conductance. Longitudinal resistivity revealed saline, brackish, and fresh aquifers with a wide range of values based on local hydrogeological setting and wells information. Generally, a high ρ L suggests freshwater mainly with gravel-sand, while a low ρ L reveals saline water with clay as the dominant lithology. In this investigation, longitudinal resistivity delineated the saline aquifer with ρ L < 5 Ωm having clay, the brackish aquifer with ρ L from 5 to 20 Ωm containing clay-sand, and freshwater with ρ L > 20 Ωm having sand, gravel-sand, and gravel ( Table 2). Mapping of the longitudinal resistivity over 45 VES is given in Figure 9a Figure 9 suggests that longitudinal resistivity increases (from 34 to 432 Ωm) with depth within the fresh aquifer and decreases (from 4.1 to 0.1 Ωm) with depth within the saline aquifer, which shows seawater intrusion from the seaside with the increasing depth (Figure 9b). 2D longitudinal resistivity maps were obtained and interpreted for saline, brackish, and fresh groundwater zones along profiles A-E (Figure 10a,b). Longitudinal resistivity varies along profile A to E with ρL between 0.1-432, 0.1-256, 0.1-321, 0.1-321, and 0.5-321 Ωm, respectively (Figure 10a). The fresh, brackish, and saline aquifers revealed by the longitudinal resistivity are given in Figure  10b. Along profile A, saline aquifer was assessed with ρL ranging from 0.1 to 3.2 Ωm for VES 11-13, brackish aquifer was revealed with ρL from 10 to 18 Ωm for VES 8-11, and fresh aquifer was interpreted by ρL between 35-432 Ωm for VES 1-10. Seawater intrusion along profile A was delineated about 45 km inside the study area along VES 8-13. Along profile B, fresh aquifer was delineated by ρL between 34-256 Ωm for VES 14-16, brackish aquifer was demarcated with ρL 8 Ωm for VES 17, and saline aquifer was assessed with ρL from 0.1 to 3.  2D longitudinal resistivity maps were obtained and interpreted for saline, brackish, and fresh groundwater zones along profiles A-E (Figure 10a,b). Longitudinal resistivity varies along profile A to E with ρ L between 0.1-432, 0.1-256, 0.1-321, 0.1-321, and 0.5-321 Ωm, respectively (Figure 10a). The fresh, brackish, and saline aquifers revealed by the longitudinal resistivity are given in Figure 10b. Along profile A, saline aquifer was assessed with ρ L ranging from 0.1 to 3.2 Ωm for VES 11-13, brackish aquifer was revealed with ρ L from 10 to 18 Ωm for VES 8-11, and fresh aquifer was interpreted by ρ L between 35-432 Ωm for VES 1-10. Seawater intrusion along profile A was delineated about 45 km inside the study area along VES 8-13. Along profile B, fresh aquifer was delineated by ρ L between 34-256 Ωm for VES 14-16, brackish aquifer was demarcated with ρ L 8 Ωm for VES 17, and saline aquifer was assessed with ρ L from 0.1 to 3.2 Ωm for VES 11 and 17-22. Seawater intrusion of about 75 km was revealed along profile B for VES 11 and 17-22. Along profile C, saline aquifer was evaluated by ρ L between 0.1-2 Ωm for VES 20 and 26-29, brackish aquifer was evaluated with ρ L 9 Ωm for VES 26, and fresh aquifer was delineated by ρ L 37-321 Ωm for VES 23-25. Seawater intrusion of 50 km was delineated along this profile with VES 20 and 26-29. Along profile D, freshwater was revealed with ρ L ranging from 37 to 321 Ωm for VES 23 and 30-34, brackish aquifer was assessed with ρ L between 11-18 Ωm for VES 31 and 33-35, and saline aquifer was delineated by ρ L between 0.1-3.  Figure 11a provides the integration of 2D longitudinal resistivity maps, whereas the separation of saline, brackish, and fresh groundwater zones in the integrated 2D ρ L maps is given in Figure 11b. The longitudinal resistivity delineated seawater intrusion of about 25-75 km from the Arabian Sea in the southwest towards the study area in the northeast (Figure 11b). zones in the integrated 2D ρL maps is given in Figure 11b. The longitudinal resistivity delineated seawater intrusion of about 25-75 km from the Arabian Sea in the southwest towards the study area in the northeast (Figure 11b). Saline, brackish, and fresh groundwater revealed by transverse unit resistance, longitudinal resistivity, and longitudinal unit conductance shows good matching. D-Z parameters delineate seawater intrusion of about 25-90 km from the sea towards the inland. Hence, D-Z parameters demarcate three different aquifers of fresh, brackish, and saline water without an overlapping character. The study of D-Z parameters suggest that the overlapping of saline, brackish, and fresh aquifers caused by intermixing of similar lithologies (clay-sand and clay or sand and clay-sand, etc.) in VES models can be successfully removed by the use of D-Z parameters in any aquifer system. Hence, D-Z parameters provide more confident solutions for the delineation of fresh-saline aquifers.  Saline, brackish, and fresh groundwater revealed by transverse unit resistance, longitudinal resistivity, and longitudinal unit conductance shows good matching. D-Z parameters delineate seawater intrusion of about 25-90 km from the sea towards the inland. Hence, D-Z parameters demarcate three different aquifers of fresh, brackish, and saline water without an overlapping character. The study of D-Z parameters suggest that the overlapping of saline, brackish, and fresh aquifers caused by intermixing of similar lithologies (clay-sand and clay or sand and clay-sand, etc.) in VES models can be successfully removed by the use of D-Z parameters in any aquifer system. Hence, D-Z parameters provide more confident solutions for the delineation of fresh-saline aquifers.

Physicochemical Analysis
D-Z parameters computed from the resistivity measurements delineated saline, brackish, and fresh groundwater zones in the studied area. The demarcation of these aquifers was confirmed by the physicochemical analysis. Such analysis was conducted for 13 groundwater samples collected from the wells located in different parts of the investigated site. The location of groundwater samples is shown in Figure 1a. The statistical analysis of physicochemical parameters was performed for main cations namely calcium (Ca 2+ ), magnesium (Mg 2+ ), sodium (Na + ), and potassium (K + ), as well as the main anions such as chloride (Cl − ), nitrate (NO3 − ), bicarbonates (HCO3 − ), and sulphates (SO4 2− ), and the parameters including total dissolved solids (TDS), pH, and electrical conductivity (EC) ( Table 3). The suggested guideline of the World Health Organization (WHO) was used to assess the groundwater quality [42,43]. Here, we assessed the saline, brackish, and fresh groundwater by the suggested range of WHO for main cations, anions, and pH, as well as TDS and EC ( Table 3). The parameters which exceeded the suggested limit of WHO were interpreted as brackish-saline aquifers and those which did not exceed the permissible range were evaluated as freshwater.
The interpretation of the physicochemical study reveals that four samples (1, 2, 5, and 8) lie in the fresh aquifer zone, five samples (3, 6, 9, 10, and 12) delineate the brackish aquifer, and the remaining four groundwater samples (4, 7, 11, and 13) reveal the saline aquifer ( Figure 12). Fresh groundwater samples (Figure 12) lie in the zone of the fresh aquifer delineated by D-Z parameters (Figures 3-11). Similarly, locations of five brackish-groundwater samples match with the brackish aquifer revealed by D-Z parameters (Figures 3-11). Moreover, saline groundwater samples are located in the zone of the saline aquifer as demarcated by D-Z parameters (Figures 3-11). The results suggest a strong correlation for saline, brackish, and fresh groundwater revealed by

Physicochemical Analysis
D-Z parameters computed from the resistivity measurements delineated saline, brackish, and fresh groundwater zones in the studied area. The demarcation of these aquifers was confirmed by the physicochemical analysis. Such analysis was conducted for 13 groundwater samples collected from the wells located in different parts of the investigated site. The location of groundwater samples is shown in Figure 1a. The statistical analysis of physicochemical parameters was performed for main cations namely calcium (Ca 2+ ), magnesium (Mg 2+ ), sodium (Na + ), and potassium (K + ), as well as the main anions such as chloride (Cl − ), nitrate (NO 3 − ), bicarbonates (HCO 3 − ), and sulphates (SO 4 2− ), and the parameters including total dissolved solids (TDS), pH, and electrical conductivity (EC) ( Table 3). The suggested guideline of the World Health Organization (WHO) was used to assess the groundwater quality [42,43]. Here, we assessed the saline, brackish, and fresh groundwater by the suggested range of WHO for main cations, anions, and pH, as well as TDS and EC ( Table 3). The parameters which exceeded the suggested limit of WHO were interpreted as brackish-saline aquifers and those which did not exceed the permissible range were evaluated as freshwater.
The interpretation of the physicochemical study reveals that four samples (1, 2, 5, and 8) lie in the fresh aquifer zone, five samples (3, 6, 9, 10, and 12) delineate the brackish aquifer, and the remaining four groundwater samples (4, 7, 11, and 13) reveal the saline aquifer ( Figure 12). Fresh groundwater samples ( Figure 12) lie in the zone of the fresh aquifer delineated by D-Z parameters (Figures 3-11). Similarly, locations of five brackish-groundwater samples match with the brackish aquifer revealed by D-Z parameters (Figures 3-11). Moreover, saline groundwater samples are located in the zone of the saline aquifer as demarcated by D-Z parameters (Figures 3-11). The results suggest a strong correlation for saline, brackish, and fresh groundwater revealed by physicochemical and D-Z parameters. Hence, saline, brackish, and fresh groundwater zones revealed by D-Z parameters of the VES method are validated by the physicochemical analysis, which suggests that this approach can be successfully used as an alternate to the expensive laboratory tests analysis in order to assess seawater intrusion into fresh aquifers for coverage of large areas. revealed by D-Z parameters of the VES method are validated by the physicochemical analysis, which suggests that this approach can be successfully used as an alternate to the expensive laboratory tests analysis in order to assess seawater intrusion into fresh aquifers for coverage of large areas.

Discussion
This research was conducted for the demarcation of saline, brackish, and fresh groundwater in a coastal area (Bela Plain) of Pakistan. Due to the over-exploitation of freshwater resources, especially for agriculture purposes, the study area is facing a problem of seawater intrusion, and thus, it is a big challenge for the proper management of freshwater resources in Bela Plain. Seawater intrusion is a universal problem, especially in the coastal areas. Conventionally, seawater intrusion

Discussion
This research was conducted for the demarcation of saline, brackish, and fresh groundwater in a coastal area (Bela Plain) of Pakistan. Due to the over-exploitation of freshwater resources, especially for agriculture purposes, the study area is facing a problem of seawater intrusion, and thus, it is a big challenge for the proper management of freshwater resources in Bela Plain. Seawater intrusion is a universal problem, especially in the coastal areas. Conventionally, seawater intrusion into coastal aquifers is delineated by laboratory tests using the groundwater samples. However, such tests are expensive, and still cannot assess large areas. Geophysical methods especially VES techniques have been widely used for an evaluation of contaminated aquifers. The use of such methods is inexpensive, non-invasive, and user friendly. However, in some cases, the subsurface lithologies show similarities such as clay-sand and clay or sand and clay-sand, which cause the intermixing of resistivity values. Such intermixing in resistivity values of subsurface layers produce an overlapping of lithologies, and hence cause an overlapping of fresh-saline aquifers associated with these lithologies. In this work, we have used Dar-Zarrouk parameters to solve this problem. D-Z parameters namely longitudinal resistivity (ρ L ), longitudinal conductance (S c ), and transverse unit resistance (T r ) estimated from resistivity models delineate the fresh-saline aquifer interface without causing any intermixing of saline, brackish, and fresh groundwater. Several studies have used such parameters for the delineation of fresh-saline aquifers, however, none of these studies provides 2D mapping of the fresh-saline interface.
In addition, the 2D ERT cannot cover large areas. Therefore, this work is the first ever investigation, which provides 2D delineation of fresh-saline aquifers using the D-Z parameters computed from VES.
The results of T r , S c , and ρ L in Figures 3-11 are highly agreeable to each other, especially the fresh-saline aquifers delineated by T r and ρ L show good matching. However, there are still some minor discrepancies. For instance, the S c results do not show perfect matching with those of the other two D-Z parameters. The brackish-fresh interface of S c is delineated further inland (north/northeast) as compared to that of T r and ρ L , while the saline-brackish interface of T r and ρ L is demarcated more towards the sea (south/southwest) as compared to that of S c . This inconsistency is due to the nature of these parameters. High values of T r and ρ L with a wide range (i.e., 2000-29000 Ωm 2 and 20-440 Ωm, respectively) delineate freshwater, while low values of these parameters with a small range (i.e., <2000 Ωm 2 and <20 Ωm, respectively) reveal brackish-saline water. Conversely, low values of S c with a low range (i.e., <3 mho) delineate freshwater, whereas its high values with a wide range (i.e., 3-1200 mho) demarcate the brackish-saline aquifer. The comparison of saline-brackish and brackish-fresh interface revealed by the three D-Z parameters and the resultant interface of these parameters including the topography of the investigated area is shown in Figure 13. The resultant fresh, brackish, and saline aquifers of the D-Z parameters ( Figure 13) are in good agreement with those of the hydrogeological map of the investigated area ( Figure 1a).

Conclusions
In this study, the calibration between the formation resistivity computed from 45 VES models and lithological logs constructed from 13 boreholes constrained subsurface units into six layers, such as clay having resistivity <7 Ωm with saline aquifer, clay-sand having a resistivity range of 4-24 Ωm with brackish aquifer, sand having resistivity between 18-42 Ωm with fresh aquifer, gravel-sand This innovative geophysical approach is applicable in any aquifer system under any hydrogeological setting. This investigation clearly marks the fresh-saline boundary (i.e., intermixing of fresh and saline aquifers), and proposes that the inland areas in the north and northeast of the nearby brackish-fresh boundary are the most sensitive areas, which can be easily intruded by seawater in the future. A 2D/3D ERT survey is suggested to be conducted in the sensitive areas near the brackish-freshwater interface for future management of freshwater resources in Bela Plain.

Conclusions
In this study, the calibration between the formation resistivity computed from 45 VES models and lithological logs constructed from 13 boreholes constrained subsurface units into six layers, such as clay having resistivity <7 Ωm with saline aquifer, clay-sand having a resistivity range of 4-24 Ωm with brackish aquifer, sand having resistivity between 18-42 Ωm with fresh aquifer, gravel-sand with resistivity between 34-60 Ωm containing fresh aquifer, and gravel with resistivity greater than 50 Ωm having fresh aquifer. This calibration shows the intermixing of resistivity values for subsurface layers, which cause an overlapping of fresh-brackish and brackish-saline groundwater zones. Based on the local wells information and hydrogeological conditions, D-Z parameters delineated saline, brackish, and fresh groundwater with a wide range of values. The saline aquifer was delineated by T r < 1000 Ωm 2 , ρ L < 5 Ωm, and S c > 24 mho, brackish groundwater was revealed for T r between 1000-2000 Ωm 2 , ρ L range of 5-20 Ωm, and S c between 3-24 mho, and fresh aquifer was demarcated by T r > 2000 Ωm 2 , ρ L > 20 Ωm, and S c < 3 mho. These parameters delineate seawater intrusion of about 25-90 km from the Sonmiani Bay towards Bela Plain. Saline, brackish, and fresh groundwater demarcated by 2D mapping of D-Z parameters was confirmed by the physicochemical analysis of 13 groundwater samples and the hydrogeological map of the investigated area. The results confirm that the qualitative analysis of D-Z parameters is consistent with the regional geology of the study area. Hence, D-Z parameters can delineate seawater intrusion into fresh aquifer through clear demarcation of 2D/3D mapping of fresh-saline aquifers, and thus reduce the number of expensive laboratory tests. This economical approach is effective in both homogeneous and heterogeneous aquifers.