Abstract: Deterioration of water quality in the coastal zones of Bela plain of Pakistan due to saltwater infiltration into the freshwater aquifer has become a major concern. The present study was conducted to estimate the depth of the fresh water-salt water interface using appropriate analytical techniques to data obtained from test boreholes. The theoretical values were not valid for the study area especially where the sedimentary deposits are not homogeneous or isotropic, commonly assumed for such solutions. Zones of fresh and salt water are stratified and a maximum chloride value occurred near to the surface indicating salt water intrusion in the upper layers. The lower layers due to their low permeability appear to be resisting the inland encroachment of sea water. The extent of intrusion was related to the rate of ground water outflow from the aquifer to the sea. The numerical and analytical results were compared and showed good agreement.
INTRODUCTION
Sea intrusion into coastal aquifers has become a major concern (Batayneh, 2006) due to the changes in quantity and quality of fresh groundwater discharged to coastal salt water systems. Therefore, understanding of saline water intrusion is essential for the management of coastal water resources (Ginzburg and Levanon, 1976). The location of fresh water-sea water interface in a coastal aquifer is highly significant for use and management of available water resources. The flow of freshwater in the aquifer towards sea is necessary to avoid salt water from encroaching farther into the aquifer (Cooper et al., 1964). When a fresh water aquifer is connected to a surface body of salt water, such as the sea or a salt lake, the denser salt water will tend to penetrate some distance into the aquifer until a state of balance with fresh water is established which is called as equilibrium state of sea water intrusion.
Previously, salt water intrusion into coastal aquifers has been the focus of many researchers. A number of researchers have been engaged in the area of seawater intrusion. Bear and Dagan (1964) and Shamir and Dagan (1971), among others, published solutions based on a sharp interface assumption. Mualem and Bear (1974) considered the effect of discontinuous stratifications on the shape of the interface.
| Fig. 1: | Location map of Bela plain |
Henry developed the first analytical solution that includes the effect of dispersion. Pinder and Cooper Jr. (1970) used the method of characteristics to solve the convective-dispersive transport equation with a constant dispersion coefficient. Lee and Cheng (1974) used stream functions to obtain a steady-state solution of the same equation. In working with confined coastal aquifers, most researchers considered the existence of an exit face that connects the aquifer directly to the sea.
The study of this problem arose from the increasing demand of subsurface water supplies, especially in the arid areas, due to domestic and irrigation projects adjacent to the coastal regions around the world. To manage the water resources amicably in coastal areas, it is pertinent to predict the changes in groundwater quality with special emphasis on water salinity within the aquifer due to any future abstraction plans. Proper future planning and regular monitoring of water quality changes in coastal areas is a pre-requisite to implement feasible correction measures (Mercado, 1985).
Location of study area: Beta Plain of Lasbela district covers an area of about 3100 km2. It is located in the southern part of Baluchistan between the latitude 25°15'-26°25' North and longitude 66°0'-66°55' East (Fig. 1). The area is surrounded by mountains in the east, west and north Miani Harbor, a bay of Arabian Sea, is the southern boundary. It lies in an arid zone with mean annual rainfall between 75-125 mm, potential evaporation is very high thus resulting in very low recharge to the ground water.
Agriculture is the principle economic source in the area which depends on rainfall, stream flow and open wells. The streams in the area are mostly intermittent and Porali River is the only perennial one, located on the southern side. It is bounded by sedimentary rocks on the east, west and north while in the south, it extends to the Arabian sea. The floor of the plain constitutes the unconsolidated deposits of quaternary age. The topography of the area has been formed by erosion.
Bela Plain consists of five main physiographic units i.e. Mountains, piedmont area, flood plain, Porali plain and rolling sand plain (Fig. 2). The main mountains are: Mor Range (north-eastern side), Par Range (south-eastern side), Halar Range (north-western side) and Haro Range (south-western side). These mountains are composed mainly of dense, indurated, sedimentary and igneous rocks.
Hydrology: In the piedmont and foothill areas of Bela Plain, water is derived directly from precipitation and the stream flow originates from the mountains (Fig. 3). In the area underneath and adjacent to Porali River and other streams, the ground water is recharged by stream runoff through permeable beds. Whereas, in the remaining part of the area, ground water is recharged mainly from percolation of precipitation.
Porali River is the main stream of the region. Its main sources of water are in Mor and ParRanges in the north, with drainage area of 4040 km2 at Sinchi Bent gauging station. The discharge measurements of streams and rivers are done by Surface Water Hydrology, WASID (WAPDA, 1973). The average annual flow was 120 million m3 from 1962-1971 and is mostly used for irrigation. The Kud River flows into the Porali about 13-16 km downstream from the Sinchi Bent gauging station. The average annual flow was 102 million m3 recorded at Mai Gondrani from 1962-1971. Windar Nai is a major non-perennial stream with a drainage area of 1735 km2 and flows directly into the sea. Another prominent non-perennial stream is Khantra Nai.
Test drilling revealed that the most permeable and extensive aquifer in the unconsolidated deposits is gravel that underlies the alluvial fans and piedmont areas. The groundwater in major part of Bela Plain is generally unconfined except in the south-western part where it is confined beneath with semi-permeable beds of silt and clay. The depth to water table varies from 11 m (near the sea) to 260 m (near the mountains).
In south western and almost all the middle portion of the study area, the water table depth zones are concave in shape along Khantra Nai and Porali River (Fig. 3). The conditions of recharge depend on the seasonal surface water discharge into the existing ground water. Almost all over the Bela Plain, the direction of ground water movement is in conformity with surface water flow (Fig. 4). In the south-eastern part of the plain, the flow lines are converging into a narrow strip which is indicative of very low hydraulic conductivity of this area.
Chemical quality of water: Water samples were collected from 25 shallow wells from the whole study area. The distribution of total water salinity (TDS mg L-1) into different zones is presented (Fig. 5). The TDS values of groundwater ranged from 280 to 9080 mg L-1. Out of 25 total samples, the TDS values of only 12 water samples were less than 1500 mg L-1 and is useable. The remaining 13 water samples were highly saline and unfit for irrigation (WAPDA, 1973). In general, the shallow water in eastern part of Bela Plain, beginning from Naka Khari to area north of Bela town, is fresh with low salt concentration.
| Fig. 2: | Physiographic features of Bela plain |
| Fig. 3: | Map of water channels and water table depth zones |
| Fig. 4: | Flow net of Bela plain |
| Fig. 5: | Map of chemical quality of ground water in test holes wells |
MATERIALS AND METHODS
Two approaches were followed for the study. First dealt with the salt water-fresh water interface in homogeneous and isotropic aquifers and the second dealt with the layered aquifers.
Interface in homogeneous and isotropic aquifers: Several analytical methods were reviewed to locate the fresh water and saline water interface in homogeneous and isotropic aquifers (Table 1). Based on the data available and the simplicity of the methods, three methods were selected for use in the study. In order to apply the analytical methods, the study area was divided into three zones based on the topography and aquifer thickness (Fig. 6). Due to the non-availability of the bed rock data, the aquifer thicknesses was assumed based upon the lithological data.
| Zone 1: | It is the western portion. Its left boundary is Haro Range and right boundary crosses Lasra Dhora while running parallel to Khantra Nai |
| Zone 3: | It is the eastern portion which is in fact the Windar Valley, with Windar Nai flowing in it, Pab Range and a small portion of Mor Range is bordering it on the right side |
| Zone 2: | It is lying in between Zone 1 and Zone 3. The outlet of all the three zones is in Miani Horbar as shown Fig. 6 |
Groundwater discharge was calculated in the three zones by applying Darcy’s Law:
Q = K•I•A
where, Q is the discharge (m3 sec-1), I is the hydraulic gradient of the ground water (%) and A is the area of the aquifer (m2):
A = K•H = (Length•Thickness)
where, H is the thickness of the aquifer in m, L is the length of the aquifer in m:
Q = K•I•L•H
| Table 1: | Classification of methods for determining the interface location |
| Fig. 6: | Map of saline water classification zones |
| Table 2: | Hydrologic characteristics of aquifer in the three zones |
| *Discharge per unit length of the aquifer | |
So, the discharge per unit length of the aquifer becomes:
q = Q/L = K•I•H
Groundwater discharge to the sea calculated by Darcy's Law, total thickness of the aquifer, hydraulic conductivity and the hydraulic gradient of each zone is given in Table 2. Hydraulic gradient was estimated with the help of water table elevation contours.
In the present study, the following analytical methods were used and the results obtained by these methods were represented in graphical form Herzberg (1901) equation.
The study assumed a static equilibrium between aquifer fresh water and sea water and a sharp interface around the middle of the mixed region. Considering the equilibrium between pressures from either side on the interface at point P gives:
Therefore:
where, D is depth of interface below the mean sea level, λf is fresh water specific weight, taken as 1 g cm-3, γS is salt water specific weight, taken as 1.025 g cm-3, h is water table altitude.
This leads to:
D = 40 h
This relation is valid only under the following conditions:
| • | Salt water must be stagnant, so that the head of the salty ground water below the interface is equal to that of sea |
| • | Flow in the fresh ground water is horizontal |
| • | Water table or piezometric surface lies above the sea level and slopes downward toward the sea |
Where the groundwater flow is nearly horizontal, the G-H relation gives satisfactory results. Only near the seashore, where vertical components of water flow become pronounced, causes significant errors in the position of the interface occur (Todd, 1980).
For applying G-H method, water table profile up to the sea was needed which is not available for the study area. In order to apply the above method, Dupuit equation was applied to estimate water table profiles in the three zones. Dupuit equation requires the head difference in ground water and the respective distance which was solved by the following equation:
where, h(x) is the required water table elevation at any point x, H and h are the upstream and downstream head of ground water respectively, L is the horizontal distance existing between these heads.
Kashef method (Kashef, 1983): To determine the configuration of the salt water-fresh water interface in the coastal unconfined aquifers, some selected rigorous mathematical solutions were adopted according to Kashef (1983). For the zones very close to the sea shore, two simple equations were derived on the basis of rigorous solutions according to Charmonman (1965) and Kashef and Csallany (1977).
In the case of horizontal fresh water outflow to the sea, the solution was deduced from the solution of Charmonman (1965) and for vertical outflow face, the solution was deduced from the author's own solution (Kashef and Csallany, 1977). He also provided a solution for interface location by modifying the G-H equation.
Procedure:
| • | First of all discharge per unit length is determined by applying Darcy’s Law |
| • | Then M is determined by applying the equation: |
where, α is the excess specific gravity:
where, γf is Fresh water specific weight, as 1 g cm-3, γS is Salt water specific weight, as 1.025 g cm-3, α is 0.025, L is the total length; the sea has intruded into the aquifer. With the value of N calculated, L can be found out. x is taken as horizontal distance inland from the sea shore. Then following equations are used for locating the interface.
| • | When the fresh water outflow to the sea is vertical, the solution was deduced from the solution of Kashef and Csallany (1977) as given by the following equation: |
where, D is depth to interface (outflow face being vertical):
| • | In the case of horizontal outflow, the solution was derived from the solution of Charmonman (1965) and was: |
DH is depth to interface (outflow face being horizontal).
| • | Solution derived, by using G-H approximation: |
Hence, it gives the modified G-H interface.
The values of DV and those of DH are always larger than D (Kashef, 1983).
Gupta and Gaikwad (1987) developed a model to predict the equilibrium location of upconed interface due to a horizontal well in an unconfined coastal aquifer of finite thickness. In the present study only that portion of the model was considered which accounted for a specified natural flow condition i.e., the interface profile and length of intrusion are determined.
In this model an unconfined coastal aquifer of uniform thickness was considered. The following assumptions were considered:
| • | Ground water flow is steady |
| • | Aquifer is homogeneous and isotropic |
| • | Tidal and barometric fluctuations of the water level are neglected |
| • | No infiltration to the aquifer |
| • | Interface is sharp |
The following procedure was adopted to predict the equilibrium location of upconned interface.
Procedure: The length of the domain where the interface intersects the base of the aquifer is determined from the equation:
where, q is the fresh water discharge to the sea per unit length of the aquifer, γf is fresh water specific weight, as 1 g cm-3, γS is salt water specific weight, as 1.025 g cm-3, α is 0.025.
This equation is then used to find the interface profile for the natural flow condition:
Where:
| β | = | 1, for horizontal outflow face |
| β | = | 0.55, vertical outflow face |
| D | = | depth to the interface |
| = | X/L | |
| M | = | L/H |
RESULTS AND DISCUSSION
G-H method (Herzberg, 1901): Data in Table 3 show some selected points in the three zones and the depth of salt water-fresh water interface calculated by G-H method at these points. These points were selected for comparison with the other methods used in the analyses.
In zones I and 3, the head difference of 6 m was considered. The difference between the two heads cannot be taken more than 6 m because G-H relation simply multiplies the water table elevation with 12.2 m, i.e., 6.1x12.2 m = 243.84 m, whereas, thickness of the aquifer in these zones is assumed to 215 m only. Similarly the head difference in zone 2 is taken as 10 m.
Water table profiles and the respective depths to the interface in the three zones are shown in Fig. 6. The largest sea water extent of 13.7 km was observed in zone 3 due to the existence of flat hydraulic gradient as compared to the other zones. Because, flat gradients resulted from high hydraulic conductivities allow greater sea water intrusion.
| Table 3: | Calculation of fresh water and sea water interface position by Ghyben-Herzberg method at some selected points |
| *Total length of intrusion | |
| Table 4: | Calculation of fresh water-sea water interface position by Kashef’s method at selected points |
| Table 5: | Calculation of freshwater-seawater interface position by modified ghyben-herzberg method at some selected points |
| Table 6: | Calculation of freshwater-seawater interface position by Gupta and Gaikwad method at some selected points |
The least sea water extent of 10.5 km was observed in zone 2. Sharp hydraulic gradients were observed here due to low hydraulic conductivity. This property of aquifer here is resistant to sea water intrusion. However, the sea water extent of 12.1 km was found in zone 1 which is between the two other zones due to the existence of intermediate hydraulic gradient.
Kashef method (Kashef, 1983): Figure 8 shows the interfaces for the three zones with the outflow face being horizontal and vertical, respectively according to Kashef Method (1983). In the graphical presentation, the two interfaces are nearly identical. But the calculated values were greater for DH than DV values, but the difference increased as the sea is approached (Table 4).
A comparison of three graphs showed that zone 2 has the deepest interface and the largest extent due to the depth and thickness of aquifer (275 m). It was reported that the extent of intrusion increases with an increase in the thickness and hydraulic conductivities of aquifer (Bear, 1979). The results obtained by Kashef method are given in Table 4. The sea water interface calculated with modified Ghyben-Herzberg formula by Kashef (1983) for some selected points are given in Table 5 (Gupta and Gaikwad, 1987).
Data in Fig. 9 show the interfaces for three zones with the outflow face being horizontal and vertical, respectively according to Gupta and Gaikwad (1987). The results were almost identical to those obtained by Kashef (1983). The results of some selected points is Table 6 indicated that the modified G-H formula by Kashef (1983) shows the interfaces at shallow depths as compared to those found by Kashef (1983) and Gupta and Gaikwad (1987) including also the G-H method (Herzberg, 1901). Although the calculated values by this method are far greater than those by the other method, the depths to interfaces in three zones are:
DG-H<DG-HMOD<DGUP<DKAS
where, DG-H is depth to interface by G-H method (Herzberg, 1901), DG-HMOD is depth to interface by Modified G-H method by Kashef (1983), DGUP is depth to interface by Gupta and Gaikwad (1987), DKAS is depth to interface by Kashef (1983) method.
A comparison of the interface locations as calculated by three methods is presented in Fig. 7-9 (Modified G-H, Kashef and Gupta). Data in Fig. 6 shows the areal extent (zones) of saline water intrusion as calculated by the analytical methods.
The G-H Method (Herzberg, 1901) showed the largest area for sea water intrusion, whereas, Kashef (1983) and Gupta and Gupta (1987) methods showed the identical areas on this scale i.e., 2.5 cm and 10 mL. But the study of results obtained by Kashef (1983) showed less extent of sea water intrusion than that determined by Gupta and Gupta (1987). The estimated extent of saline water intrusion applying G-H method for homogeneous conditions is presented in Table 7. It was assumed that the permeable layers of zone 1 and zone 2 are 150 m thick and that of zone 3 is 180 m thick.
Chemical analyses helped to set some criteria in order to check the variability and extent in aquifer water quality due to sea water intrusion. An increase in chloride (Cl) content indicated sea water intrusion. However, the variations in the concentration of other ions were frequently ignored due to their modification resulting from the interaction with aquifer minerals (Mercado, 1985).
| Fig. 7: | Profiles of Ghyben-Herzberg (G-H) interface in classified saline zones |
| Fig. 8: | Profiles of Kashef interface in classified saline zones |
| Fig. 9: | Profiles of Gupta-Gaikwad interface in classified saline zones |
| Table 7: | Calculation of freshwater-seawater interface position by Ghyben-Herzberg method at some selected points assuming two layered aquifer |
| Table 8: | Mean values of Anion (WAPDA) |
| The average anion values are taken for each well in ppm | |
Harris (1967) stated that Cl contents in water in excess of 50 mg L-1 indicate salt water contamination. Similar results were reported by Studies in stratified aquifer of Barrier Island by Chow (1964) who suggested Cl content in average sea water to be 18,930 mg L-1 as stated by Reilly and Goodman (1985).
WAPDA (1973) analyzed the chemical quality of all the wells of Bela Plain. They selected 5 wells which were found to be contaminated by the sea water intrusion on the basis of the order of abundance of the following anions (Table 8):
| • | Chlorides (Cl) |
| • | Sulphates (SO4) |
| • | Bicarbonates (HCO3) |
It was pointed out that while considering the sea water contamination, the Cl to HCO3 ratio is important to diagnose sea water intrusion. Because sea water is low in HCO3 so an increase in this ratio might be due to an increase in Cl contents resulting from sea water. But no criteria was set up for this ratio, these ratios are shown in Table 8.
According to Termblay et al. (1973), the ratio of Cl contents to TDS in the sea water remains almost constant when the sea water is mixed with fresh water. They considered Cl content of 20 mg L-1 as belonging to uncontaminated ground water, Cl content greater than 40 mg L-1 indicates salt water contamination and Cl content in excess of 100 mg L-1 was considered to be a part of the zone of diffusion.
According to Todd (1980), the Cl concentration in natural water is commonly less than 10 mg L-1 (milligram per liter) in humid regions but up to 1000 mg L-1 in more arid regions; about 19,300 mg L-1 in sea water and as much as 200,000 mg L-1 in brines.
Based upon the above discussion, the groundwater or freshwater having Cl content of 10,000 mg L-1 or more will be considered as being contaminated by sea after intrusion and the water with Cl content equal to or less than 1000 mg L-1 fresh water.
Considering the chemical quality of ground water in well No. 26 and 4 (Fig. 10), the transition (in horizontal as well as in vertical direction) from fresh to saline water is observed as sea is approached. In the case of well no 26, both the Cl and TDS values concentration decreased with depth. The values near the surface were 16873 and 35800 mg L-1 for Cl and TDS, respectively which were very close to sea water quality. Whereas, the values of Cl and TDS were 1180 and 2540 mg L-1, respectively at a depth of 185 m.
The distribution of Cl contents with depth in the groundwater in test well No. 26 is shown in Fig. 10. A zone of transition existed from fresh water (chloride content of 1000 ppm) to saline water (chloride content of 10,000 ppm).
A similar pattern of water quality deterioration except the value of Cl shows that well 4 lies almost in the fresh water zone. The Cl contents rose to 3445.9 mg L-1 at a depth of 55-60 m. As the location of well No. 4 is far away from the sea and near to Porali River, the decrease in Cl values might be due to the flushing of this zone by the fresh water from Porali River. However, Price (1996) reported that in arid areas, the potential evapotranspiration generally exceeds the rainfall for most of the year; the aquifer outcrop may be a zone of concentration rather than dissolution of salts. In this way high salinity level in ground water can be formed at or near the surface.
Well No. 13 showed a decrease in Cl contents with depth which might be attributed to the movement of high saline water to aquifer due to high permeability of aquifer layers in this zone. Well No. 7 located along the Titian Nai showed low concentration of Cl (Fig. 10) ranging from 308-1212 mg L-1 which lies almost in the fresh water zone. Well No. 16 located far inland shows low Cl contents ranging from 524.6-861.4 mg L-1 which lies in the fresh water zone. Also, Thoriako Dhora, originating from the Mor Range, might be responsible for improving the chemical quality of this well.
Well No. 14 is located above well No. 16 (Fig. 10). might have caused contamination in the latter well due to abstraction of water. But it seems that the area here is being excessively recharged by the numerous streams from the Mor Range which checks the contamination.
At the location of well No.18, the ground water is confined below thick shale layer. There was found a zig zag pattern of Cl concentration from a depth of 30-80 m below the sea level. This pattern is due to the permeability of the formations encountered in this well i.e., a permeable layer is present at 75 m depth which abruptly increased the Cl value to more than 12000 mg L-1.
| Fig. 10: | Chloride distribution in different wells |
However, based on the criteria for groundwater contaminated with sea water, it can safely be concluded that the water of this well is being affected by sea water intrusion through the permeable formations. Because, low chlorides contents were found in less permeable clay layers due to increased resistance of the clay layers.
Well No. 17 is located at a distance of 3.5 km from sea and is 1.5 km south of Windar Hal. The depth to water table is 14.5 m and high evapo-transpiration may be responsible for high salinities at shallow depth.
Now consider that the Cl content is attaining the maximum value of 15395 mg L-1 in the layer of coarse to medium sand with gravels at a depth of 67-70 m below mean sea level. The high Cl content (15395 mg L-1) might be indicating the presence of interface penetrating the pervious layer. Furthermore, operation of well No.1 might have extended the zone of dispersion in well No. 17 or the tidal fluctuations, that are more pronounced near the shore, also directly contribute to variations in chloride content (Termblay et al., 1973).
Well No. 8 is located about 6.4 km inland of well No. 17, 2.6 km north of Mobar Dhora and 3.8 km south of Windar Nai. The well water in not affected by sea water intrusion which might be attributed to its greater distance from the sea and flushing of permeable zones by the surface sources or fresh ground water (Table 9).
The convex shape of the water table contours along Windar Nai indicated that this stream is recharging the ground water (Fig. 4). This is clearly shown from the good quality of ground water at shallow depths in this area.
| Table 9: | Comparison between G-H interface and the available data |
The zonal interpretation of water quality data in the three important aspects of the pattern of saline water intrusion in the area is given below:
| • | Water quality deterioration towards the sea is common to all zones which is an apparent indication of saline water intrusion |
| • | The extent of sea water intrusion in zone 1 and zone 2 is more than twice as compared to zone 3 which shows greater outflow from the aquifer in zone 3 compared to the others |
| • | Pattern of distribution of Cl contents does not agree to the theoretical or analytical approach due to the heterogeneity and configuration of the strata |
CONCLUSION
Application of Ghyben-Herzberg method showed the largest extent of sea water intrusion as compared to other three methods. Because, it assumed zero outflows from the aquifer. According to this method, sea water intrusion was observed up to 6 km inland. The data of well No. 3 and 17 indicated that the depth to water contamination is occurring at about depths where Ghyben-Herzberg interface was located.
The high Cl contents found in the water of shallow well No. 26, 13 and 18 might have been transported from the chloride rich soils of the sea shore. In conclusion, in the non-homogeneous sediments of the study area, the depth to the fresh water-sea water interface responded to the differential flushing of stratified sediments of different permeability. Also, no regular pattern of chlorides is found in the ground water. In the presence of surface fresh water sources, high quality ground water was found in the subsurface permeable formations and relatively saline water in the impermeable formations. In those areas, where the recharge is less or absent, sea water seems to be intruded through the permeable layers, but the impermeable layers were least affected by the sea water.