Gravity Analysis for Subsurface Characterization and Depth Estimation of Muda River Basin, Kedah, Peninsular Malaysia

Abstract

Gravity survey is one of the passive geophysical techniques commonly used to delineate geological formations, especially in determining basement rock and the overlying deposit. Geologically, the study area is made up of thick quaternary alluvium deposited on top of the older basement rock. The Muda River basin constitutes, approximately, of more than 300 m of thick quaternary alluvium overlying the unknown basement rock type. Previous studies, including drilling and geo-electrical resistivity surveys, were conducted in the area but none of them managed to conclusively determine the basement rock type and depth precisely. Hence, a regional gravity survey was conducted to determine the thickness of the quaternary sediments prior to assessing the sustainability of the Muda River basin. Gravity readings were made at 347 gravity stations spaced at 3–5 km intervals using Scintrex CG-3 covering an area and a perimeter of 9000 km² and 730 km, respectively. The gravity data were then conventionally reduced for drift, free air, latitude, Bouguer, and terrain corrections. These data were then consequently analyzed to generate Bouguer, regional and total horizontal derivative (THD) anomaly maps for qualitative and quantitative interpretations. The Bouguer gravity anomaly map shows low gravity values in the north-eastern part of the study area interpreted as representing the Main Range granitic body, while relatively higher gravity values observed in the south-western part are interpreted as representing sedimentary rocks of Semanggol and Mahang formations. Patterns observed in the THD anomaly and Euler deconvolution maps closely resembled the presence of structural features such as fault lineaments dominantly trending along NW-SE and NE-SW like the trends of topographic lineaments in the study area. Based on power spectral analysis of the gravity data, the average depth of shallow body, representing alluvium, and deep body, representing underlying rock formations, are 0.5 km and 1.2 km, respectively. The thickness of Quaternary sediment and the depth of sedimentary formation can be more precisely estimated by other geophysical techniques such as the seismic reflection survey.

1. Introduction

Kedah is one of the northern states in Malaysia that is very well-known for its main and largest river, known as the Muda River. The Muda River flows from the East, channeling through Baling, Sik, and Sungai Petani districts down to lower land until the Strait of Malacca in the west. The water resources in the Muda River basin have been established over the year and consumed for irrigation and as potable water, mainly for Kedah [1], and the supplies have been extended to other states like Penang [2] with almost 80% dependency and at the Southern area of Perlis. Freshwater in this basin is very crucial for agricultural activities in Kedah and the industrial sector in Penang, as well as for domestic usage [3].

A decade ago, Malaysia had sufficient water supplies, however, the country has recently suffered from a shortage of clean water for domestic usage. This is due to increasing consumption and high-water demand exceeding the capacity of the river basin, especially in the escalation of development sectors and the growth trend of agricultural industries. Ghani et al. [4] reported that the Muda area experienced a moderately increasing rate of water usage from three major activities from 1984 to 2001, these are land clearing (1.2%), agriculture (1.35%) and urbanization (4%). In addition to the expansion of population density, the issue of water pollution also threatens the supply of clean water to certain areas and is affecting the ecology, especially in Kedah. In 2017, numerous complaints were reported regarding the failure of clean water and water shortage supply to rural areas such as Baling and Sik, Kedah. To overcome the issues, the Malaysian Government embarked on groundwater exploration, exploitation, and sustainability plans to optimally utilize the groundwater reserve.

What is still not clear, and has yet to be quantified, is the depth estimation of the Muda River basin. Even though there are many methods that can be applied to investigate the subsurface parameters particularly used as input in groundwater basin modeling, geophysical techniques are considered as among the best. However, an appropriate geophysical method to be carried out for an oriented investigation depends solely on the purpose of the study. Therefore, the aim of this study is to utilize the gravity method to investigate the subsurface structure and depths to sources of anomalies through extended gravity analyses, primarily by radially power spectrum and 3D Euler deconvolution. It is used to determine the subsurface properties according to anomalies that indicate the rocks’ density variations [5]. This potential method was employed because of its quick-moving survey, extensive coverage area and it is less expensive compared to other geophysical methods [6].

2. General Geology

The Muda River catchment area, which occupies more than 45% of Kedah state (Figure 1), is considered one of the major river basins in Malaysia. It is in the southern part of Kedah and extends from north Perak in the south to north Penang and south Thailand in the west and north. The upper tributaries drain into mountainous areas with maximum elevations of about 1860 m above sea level. The lower reach of the Muda River catchment, comprising of about 45% of the total area, generally has a lower elevation ranging from a minimum 2.5 m at the river mouth to about 70 m in hilly areas towards the east [1]. Through this geology and geophysics link, the results revealed new information such as depth estimation, structural features, and lateral distribution of the Muda River basin.

Geologically, the Muda River basin consists of four major sedimentary formations (Figure 2) namely Mahang (Lower Ordovician-Upper Devonian), Sungai Petani (Lower Silurian), Semanggol (Lower Permian-Upper Triassic) and Saiong (Jurassic-Cretaceous). These are summarized in a stratigraphic chart as shown in Figure 3. However, almost 40% of the study area is occupied by Main Range granite in the center and Bintang Range granite in the southeastern parts. Mahang formation is mainly distributed in the southwest of the study area, consisting of arenite and silicate facies, argillitic facies and Limestone facies [7]. The formation can be correlated with Sungai Petani and Baling formations, which are separated by a granite body [8]. Semanggol formation comprises of conglomerate, rhythmic and chert facies [9]. However, in contrast, Foo [10] revealed that Southern Semanggol formation in north Perak consists of only two dominant lithofacies without chert facies, similar findings are also supported by Sajid et al. [11]. In northeast Kedah, the Saiong Bed was unconformably overlying the Semanggol formation [12] and distributed to southeast Thailand. A few Kodiang limestone outcrops in the northwestern part of the Muda River basin have a similar geological age as the Semanggol formation [13,14,15,16,17].

The most significant major fault in mainland Kedah is the NW-SE trending Bok Bak Fault zone. The Bok Bak Fault can be tracked from southwest Kelantan, passing through central Perak and Bukit Perak in Kedah and it terminates close to Pokok Sena in Kedah [18]. The Bok Bak Fault shows a distinct sinistral strike-slip displacement [19]. Burton [8] predicted the net displacement was 55 km. However, Raj [20] and Abdullah et al. [21] predicted the displacement to be 20 km and 25 km, respectively. Burley and Othman [22] conducted gravity studies on the opposite sides of the Bok Bak Fault and predicted the displacement to be 30 km based on the similarities in the steepness of granite boundaries.

Figure 3. Stratigraphy chart of the Muda River basin. Modified after [23].

Figure 3. Stratigraphy chart of the Muda River basin. Modified after [23].

3. Materials and Methods

A land gravity survey was carried out in the study area with a total of 374 gravity reading stations arranged at a spacing interval of 4 to 5 km covering an area of 9000 km² (Figure 4). The coverage, however, was hindered in the inaccessible high topographic and private areas at the northeastern part. The gravity survey was conducted by Burley and Othman [22] using the La Coste & Romberg (L & R) model-G gravity meter in 1988, originally to cover the whole of Peninsular Malaysia. All gravity data were calibrated to the International Gravity Standardization Network 1971 (IGSN 71) and then reduced for drift, free air, latitude, Bouguer, and terrain to produce corrected gravity data before further analyses in generating complete Bouguer anomaly (CBA) gravity, regional CBA, residual CBA, and total horizontal derivative (THD) anomaly maps by using Oasis Montaj software for qualitative interpretation. The 3D Euler Deconvolution and power spectral analyses were also carried out for quantitative interpretation, especially in estimating depths and dips of the geological formations and faults. The summarized workflow is shown in Figure 5.

Rock density of 2.40 mg m⁻³ is used for Bouguer corrections while the accuracy of the average elevation was about 1.5 m. The elevation measurements were made using two altimeters and the ‘single base’ system [22]. The Wallace and Tiernan altimeter was used for base stations and gravity stations, respectively, and, for data validation, the altimeters were separated by less than 100 m vertically and 15 km laterally [22].

3.1. Gravity Anomaly Separation

The observed gravity anomaly, $G_{obs}$, consists of high and low frequency components representing shallow and deep geological structures, respectively. These frequency components are known as residual (high), $G_{res}$, anomaly and regional (low), $G_{reg}$, anomaly:

$$G_{obs}=G_{res}+G_{reg}$$

(1)

Some mathematical methods are required to separate the map data into two components [24]. The separation of these components must be performed in the spatial frequency domain via fast Fourier transform (FFT) known as a low-pass filter and a high-pass filter based on the method given by Harrison and Dickinson [25]. The regional gravity anomaly map is a product of a low-pass filter (long wavelength) that shows deep-seated structure while the residual gravity anomaly map is produced by the high-pass filter (short wavelength) that shows the shallow structure.

3.2. Total Horizontal Derivative Analysis

Total horizontal derivative (THD) is a well-known filtering method for potential methods to look at a different perspective of gravity anomaly, particularly to highlight a high gradient area that probably represents structural trends [23]. It measures the change rate of the gravity field in two horizontal directions (x and y), as shown in the following equation:

$$THD \left(x,y\right)=\sqrt{{(\frac{\delta g}{\delta x})}^2+{(\frac{\delta g}{\delta y})}^2}$$

(2)

where, $\frac{\delta g}{\delta x}$, $\frac{\delta g}{\delta y}$ are x and y directions, respectively.

The THD map is produced from the residual gravity anomaly map for further shallow lineament analysis that can be correlated with surface topography lineament analysis. This THD map will help in providing more lineament coverage, especially in the flat region, like the southwestern part of the study area.

3.3. 3D Euler Deconvolution

The aim of 3D Euler deconvolution (3DED) is to determine the geometry, field source location and depth estimation of gravimetric anomalies [26,27,28,29]. The 3DED has successfully been used in discovering fault lineaments distribution including the dip direction and extension below the subsurface [30]. This method was created by Thompson [31] and Reid et al. [32] who developed the equivalent method which operated on magnetic data. Later, the method has been introduced to be used in gravity data such as in Fairhead et al. [33] and Huang, et al. [34]. The 3DED equation given by [32] is:

$$N\left(\beta -g\right)=\left(x-x_0\right)\frac{\delta T}{\delta x}+\left(y-y_0\right)\frac{\delta T}{\delta y}+\left(z-z_0\right)\frac{\delta T}{\delta z}$$

(3)

where $g$ is the gravity anomaly detected at $\left(x,y,z\right)$ and $x_0$, $y_0$ and $z_0$ is the position of the gravity source. $N$ is a constant defined as the structural index (SI) depending on the type and shape of the source body such as geological contacts, faults, cylinder, or sphere (

Table 1. Structural index for gravity.

Source	Gravity S.I.
Sphere	2.0
Horizontal	1.0
Fault (small step)	0.0

). The 3DED is commonly applied in magnetic [35] and gravity interpretation since it requires a general input of gravity source geometry and requires no information of the acceleration vector. However, several SI values are required to be tested to find the optimum input to give a good clustering solution with a particular size of window. Reid, Allsop, Granser, Millett and Somerton [32] stated that window size must be small enough to differentiate anomalies resulting from different sources yet large enough to incorporate deep sources and spatially extensive sources. Chenrai et al. [36] also stated that a large window size could detect very deep basement sources.

3.4. Power Spectrum Analysis

Spector and Grant [37] had introduced a spectral factorization method to describes the frequency content of a signal based on a finite set of data [38], thus estimating the depth basements based on gravity sources that were controlled by wavenumber components. Spectral analysis of a gravity profile can be derived in the discrete Fourier Transform mathematically as the equation below:

$$\Delta g \left(x_i\right)={\displaystyle \sum }_n^N=1 [A_n cos\frac{2\pi n{x}_i}{L}+B_n cos\frac{2\pi n{x}_i}{L}]$$

(4)

where n is the total amount of data in $x$ direction, $x_i$ is the data interval in direction $x$, $n$ is the harmonic number, $A$ and $B$ are the coefficient of cosine and sine, L is the length of the cross-section of the anomaly. Finally, the average depth of basement rock was estimated by calculating the slope gradient of the power spectrum curve and substituting the result of the slope in the following relationship:

$$z=-\frac{m}{4\pi }$$

(5)

where z is the depth and m is the obtained slope from the power spectrum curve [39].

Gravity anomalies produced by a shallow structure are more dominated by high wavenumber components than those resulting from a deeper source [35]. This effect may be quantified by computing the power spectrum of the anomaly since the long power spectrum has a linear gradient whose magnitude is dependent upon the depth of the source [37].

4. Results and Discussion

4.1. Qualitative Interpretation

Complete Bouguer gravity anomaly (CBA), regional CBA anomaly, residual CBA anomaly and THD anomaly maps are shown in Figure 6a–d, respectively. The CBA map with terrain correction reflects the major geological and structural patterns. The regional anomaly map represents the long-wavelength gravity source which was produced from deep-seated bodies and was obtained by applying a low pass frequency filter of 30,000 m wavelength cut-off. The residual CBA map was produced by applying a high pass filter of 30,000 m wavelength cut-off that only allowed the short-wavelength gravity source to be retained. The THD map was produced from further filtration of the residual CBA map that marks the maximum ridge of the anomaly gradient that delineates geologic contacts or faults, which distinguish different densities of rock [40,41].

4.1.1. Complete Bouguer Gravity Anomaly Map

Based on the map shown in Figure 6a, the pattern of the Bouguer gravity anomaly can be classified into three distinctive zones of low (LG), medium (MG), and high (HG) gravity anomaly values ranging from, −30 to 0 mGal, 1 to 10 mGal, and 11 to 30 mGal, respectively. The LG regions can be observed and situated mostly in the southeastern to eastern part that consists of the Mahang Formation and granitic relief; mark the area of high topographic land as shown in Figure 4. The MG region is comprised of Mesozoic sedimentary formations, which are Sungai Petani, Kubang Pasu and Semanggol Formations, while the (HG) zone has been identified in Quaternary sediments at the western side, near to the coastal area of the study area.

Generally, the gravity anomaly values decrease laterally from southwest to northeast indicating the decrease in density due to the thinning of Mesozoic sediments towards the eastern and southeastern parts. On the other hand, the shallowing of the granitic basement as the presence of two major uplifted granite ranges, namely Main Range and Bintang Range at the eastern and southeastern parts, plays a major role in the sediment thinning. The negative gravity anomaly region, (LG), indicates the presence of Main Range and Bintang Range above the surface that have an effect of the low-density root of mountains and performed isostatic equilibrium phenomenon, causing a further gravity deficit. Consequently, the high-density earth crust compensates for these granite ranges. The cause of the low Bouguer anomaly value of the granite might be due to its less mean density compared to the mean density of the surrounding rock [22]. Loke et al. [42] reported that the gravity anomaly shows a prominent gravity minimum centered over the Main Range which is due to a thickening of granite batholith in the area. This idea is strongly supported by and reported in Ooi et al. [43].

Variations in gravity anomalies were normally caused by variation in the density of the subsurface rocks [44] and they usually indicate faults or lithological contacts [45]. As aforementioned, the negative gravity anomaly covered in the granite and Mahang formation has continued southward, indicating the presence of extended granite basement overlain by younger formations and the continuity of the Main-Bintang granite ranges. Sungai Petani, Kubang Pasu and Semanggol formations were classified by irregular medium (MG)–high (HG) gravity anomaly zones, suggesting the deposition of relatively dense sedimentary formations. The highest gravity (HG) anomaly zone was detected in the Quaternary alluvium overlying the Pleistocene Simpang formation and indicating deep sedimentary basin thickness [46].

4.1.2. Regional CBA Anomaly Map

The regional CBA anomaly map in Figure 6b has a gravity anomaly range from −12.5 to 18 mGal and can be classified as two main regions: LG (−12.5 to 0 mGal) and MG (1 to 18 mGal). The gravity reading increases towards the western part indicating the development of crustal density towards the sea. Moreover, this phenomenon has proved the statement regarding the thickening of Mesozoic sediment towards the western region. The dense contour in the central part of the map, which split the negative and positive anomalies, suggests the extension of granite and surrounding rock formation contact at depth. The zone of LG and MG has an equal part whereby the increasing gravity anomaly readings were increasing proportionally from east to west, or from high topographic to low topographic regions. However, in the vicinity of the Muda River basin, the LG zones dominated in approximately 80% of the area.

4.1.3. Residual CBA Anomaly Map

The residual anomaly map was the final corrected gravity response and was used for processing and interpretation phases [41]. Figure 6c shows the responses of high-frequency data that represents the shallow subsurface structures. The residual gravity values responded greatly in terms of sedimentary formations classification and lithological contacts compared to the CBA map (Figure 6a). For instance, the Semanggol formation has a consistent MG anomaly value (4 to 10 mGal). The gravity contour on the Bok Bak Fault; the main structure dividing Main Range granite with the Semanggol formation, Sungai Petani formation has become denser and more obvious. The four LG spots were detected to highlight a better representation of granitic bodies at high topographic regions especially at northern (Main Range), western (Jerai granite), eastern (Main Range), and southern (Bintang Range) of the study area.

4.1.4. THD Anomaly Map

The THD anomaly map was derived from the residual CBA map. Figure 6d shows a dense and conspicuous feature of steep gravity gradient in the central, northeast, and southwest parts of the study area. In general, two dominant lineaments trending NW-SE were observed at the central part of the study area. One of these lineaments was perfectly overlapped on Bok Bak Fault zone when THD and structural maps were superimposed to each other. Since one of the lineaments is perfectly matched to the Bok Bak Fault zone, the remaining lineaments close to it were also most likely associated with possible faults which have yet to be determined.

Another NE-SW strong gradient phenomenon is observed both northeast and southwest of the study area. Structural lineament analysis was then conducted by tracing the lineaments observed in both the THD anomaly map and structural maps. Trends of both lineaments are compared with topography map lineaments analysis, which represents general surface structures. The lineaments trends deduced from the topography map and THD residual anomaly map are shown in Figure 7a,b, respectively. The generated synoptic Rose diagrams on structural trends in Figure 8 inferred the delineation of the major lineaments in the NW-SE direction followed by the NE-SW direction.

The topographic lineament is, in general, quite similar to the result reported by Azman et al. [47] stressing the presence of a dominant lineaments trend in central Kedah along NE-SW and NW-SW. The THD residual anomaly structures trend resembles the analysis of Burley and Othman [22] who stated that the NW and NE trending sections are along 325° and 37°, respectively. In general, the lineaments of the THD residual anomaly map that represent shallow subsurface structures have mimicked the trend to the surface topography lineaments with less than 8° dipping angles.

4.2. Quantitative Interpretation

The Euler Deconvolution and power spectrum analyses are used for quantitative interpretation to estimate any subsurface structure. The Euler deconvolution is used for resolving the 3D components of structural elements while the power spectrum analysis is used for estimating the subsurface source depth and possible interfaces.

4.2.1. Euler Deconvolution

In this study, Euler’s solutions of the CBA map are generated using a structural index of 0.0 to delineate faults/lineaments contact with depth estimation based on Euler anomalies. A window size of 5 km × 5 km and 20 km × 20 km grid cells with a depth of tolerance of 10% were applied (Figure 9). Both Maps were overlayed on top of the Kedah state map with the known location of faults/lineaments adopted from the geological map of Kedah. A high concentration of Euler anomalies was observed along the Bok Bak Fault and northeastern part of the Muda River Basin. A moderate concentration of Euler anomalies can be spotted at the northern, western, southwestern parts of Kedah and south-southwestern parts of the Muda River Basin. The depth range of Euler anomaly sources are ranging from less than 2 km to more than 10 km.

Bok Bak Fault has the most significant number of Euler anomalies along the fault zone, with a shallow depth ranging from 2 km to 4 km at the northwestern part and increasing towards the southeastern direction with a depth ranging from 6 km to more than 10 km. Another major Euler anomaly observed in the northeastern part of the Muda River Basin is interpreted as a fault boundary separating from the Semanggol Formation and Mahang Formations with granite. This fault boundary has a shallow depth ranging from 4000 m to 6000 m separating Semanggol and Mahang Formations and increases towards the southwestern direction with a depth ranging from 6000 m to 8000 m that separates Semanggol Formation and granite. Similar Semanggol-granite fault boundaries are also observed in the southern part of the Muda River Basin trending in the NE-SW direction near the Bok Bak Fault and the southwestern part of Kedah, trending in the NNE-SSW direction. Euler anomalies are also determined in the western part of Kedah, which might represent a fault contact boundary between the Quaternary sediments and Kubang Pasu Formation.

4.2.2. Radially Power Spectrum Analysis

This analysis was conducted at 88 location points of 10 km × 10 km interval all over Kedah focusing on the Muda River basin, as shown in Figure 10. Each location point has a midpoint of a 20 km horizontal line except for four points in the southern part that were used to enhance the data. To control and validate the data, each profile was set to have 13 data points with the center of the lines representing the depth of deep and shallow sources. Figure 11 shows selected spectral analysis curves of point 1 to 8 with the natural logarithm of the power spectrum against the wavenumber that describes the deterioration or decrease in the value of the power spectrum by observing in line with all other points. The depth due to gravity sources for each profile was measured from the slope of the curves based on the radial power spectrum. The low wavenumber power spectra curves pattern was interpreted as representing a deeper source and vice versa.

Figure 12a,b shows maps of shallow to deep gravity source of the study area, respectively. In general, the depth of the shallow gravity source increases from east to west with a maximum depth of 0.76 km from the surface located in the southwestern part of the Muda River basin. A similar increasing trend pattern can also be observed on the map of deep gravity source with a maximum depth of 1.2 km from the surface. Figure 13 shows an overview of the shallow and deep gravity sources’ trend with the topography. The basin area has a bigger radius and is flatter in the shallow source at −400 m as compared to the deeper source at −1190 m, resembling a typical basic structure which is depressing downward. The red curtain shows the downward extension of Main Range Granite that is associated with the Bok Bak Fault. Figure 14 shows four profiles across the topography with shallow and deep source maps. Based on these profiles, the maximum depth area, with an average of 1.1 km, is located almost in the middle of the profile where most of the formations, lithological contact and weak zones were located. The trend of the graphs also showed the thinning of sediments towards the eastern part of the study area (coastal), if the granite bodies were not considered. These indicate that subsurface structures play a significant role in establishing groundwater basin, basement depression and such sedimentary thickness may alarm further groundwater accumulation investigation in the area.

5. Conclusions

The lineament analysis from the THD residual anomaly map and topography map that represent shallow subsurface and surface structures have shown an identical trend of NW-SE and NE-SW. The THD residual anomaly can be used to enhance the lineament data for general structural analysis, especially on the flat or non-topographic area. Two depth sources generated by spectral analysis to estimate the major sedimentary layer thickness of the basin for the presence of groundwater and its subsurface structure contribute towards developing a groundwater system. The shallow sources layer has a thickness range of 300 m to 760 m with an average thickness of 460 m, while the deep sources layer has a thickness range of 170 m to 875 m with an average thickness of 565 m.

The study outcomes show that the optimal depths were found around the most prominent structural area at the southwestern part of the Muda River basin that contained three major lithologies and formation contacts (Sg. Petani formation, Semanggol formation and Main Range granite) separated by two major faults (Bok Bak Fault and Kampung Padang Tembak fault). Hence, this area contains considerable sedimentary thickness that may suggest further groundwater investigation. However, a more comprehensive analysis of other geophysical methods is required to get a more precise result.