Surface and Crustal Study Based on Digital Elevation Modeling and 2-D Gravity Forward Modeling in Thandiani to Boi Areas of Hazara Region, Pakistan

Gravity data indicates that there is a regular relation between crustal structure, crustal density (composition), and surface ascension. In order to delineate surface and subsurface geological structure features, and to calculate the thickness variation of the crust and sedimentary/metasedimentary wedges, integrated approach of Geographic Information System (GIS) i.e. digital elevation models (DEMs) and two-dimensional forward modeling of gravity data were utilized, which provide the best results for the primary objectives. Tectonically, the study area lies in the Lesser Himalayas as well as to an extent in the sub-Himalaya, more concretely in the western limb of Hazara Kashmir Syntaxis. Topographic data was accumulated in XYZ coordinates utilizing point heights method, and DEMs generation, manipulation, interpretation, and visualization process were directed to surfer-15 and ArcGIS software. Determinately the visualization of surface geological structure in the form of DEMs were proposed. The gravity stations in single contour mode have been quantified by using Scintrex CG-5 gravity meter. The collected gravity data was processed by standardizing corrections, two-dimensional forward modeling along with gravity profile were utilized and bouguer anomaly map and gravity model was computed utilizing bouguer density of 2.4 g/cm 3 , where the subsurface structures are demarcated by the bouguer anomaly and gravity model. In summary this research has allowed the validation of surface and subsurface geological structure visualization. Digital elevation models provide a defensive prediction of the geological structure of the regional surface. The gravity model demarcated a series of stratigraphic units with density boundaries within the basement. The gravity model also suggests that the thickness of sedimentary/metasedimentary wedge in Thandiani area is 11.48 km and in Boi area, the thickness elongates to about 14.43 km. The total thickness of crust in Thandiani and Boi area is 49.53 km and 52.43 km respectively.


Introduction
The gravity method has numerous deals, which investigates anomalies caused by changes in the physical properties of subsurface rocks and requires fundamentally comparable interpretation techniques [1]. Especially gravity surveys have been widely used for crustal study, subsurface structure characterization and visualization, and proved to be a useful tool [2]. In addition, digital elevation modeling (DEMs) can be used to recognize surface geological structural trends, considering the crustal structure is closely related to the surface elevation through isostasy phenomena. The digital elevation models (DEMs) are one of the fundamental spatial datasets in geographic information systems (GIS). As sharp fluctuations in elevation is a sign of changes in the geological structure. Commonly, modeling surface structure can partly uncover the depth structure of the exploration target [3][4][5][6]. Crustal thickness diversifies significantly across continental margins, knowledge of the location of the crust-mantle interface (Moho), and the structure of the crust is imperative for understanding and evaluation continental rifting processes, geological evolution, and modern tectonic conditions. Gravity data usually investigate a regular relationship between crustal structure, crustal density (composition), and surface elevation. Geodetic computation has been made in the Himalayas for over 125 years. It has been great importance in the ascertainment of the deep structure of the earth's crust [7,8]. Gravity anomalies have been widely used to study the northwest Himalayan crustal configuration [8][9][10][11][12][13]. In the northern Pakistan, gravity measurement has been done principally by various researches [10,12]. Malinconico has conducted gravity and magnetic survey across the Main Mantle Thrust (MMT) and Main Karakoram Thrust (MKT); He was the first researcher who interpreted two-dimensionally the terrain-corrected gravity data. He suggested that the Indian plate is being thrust under the Eurasian plate and that the primary structure and this region dip towards the north. The northern structure being older, dip steeper than those of the south, which are younger [14,15]. Farah analyzed the Pre-Cambrian basement of the Punjab plain and Potwar plateau by gravity modeling in the south. It is suggested that the shield elements with increasing sediment cover can be followed by the northward movement underneath the southern thrust of Himalayas [16,17].
The Thandiani to Boi area of Hazara region lies in the western limb of Hazara Kashmir Syntaxis (HKS) which is highly folded and faulted and developed due to the stresses affected by the collision of Indian and Eurasian plates [18]. The present study area bounded by latitude 34°18′600″ N and longitude 73°26′11.99′′ E (Figure 1). This study is predicated on the integrated methods of digital elevation modeling (DEMs) and two-dimensional forward modeling of gravity data, proposing to identify the surface and subsurface structure trends to calculate the thickness of the crustal and sedimentary/metasedimentary wedge of the study area. Predicated on two-dimensional forward modeling method gravity anomaly was computed using GM-SYS profile modeling from Bouguer anomalies map of Thandiani to Boi area of Hazara region, Khyber Pakhtunkhwa (KPK) Pakistan. Projected gravity data have been practiced to generate a gravity model to illustrate the relationship of known geology with the gravity anomaly. The model parameters (e.g., dimensions, density, and depth) are adjusted until the best fit between the calculated gravity and observed anomalies were achieved.

Geological Setting
The regional geology of the predicted region consists of sedimentary, metamorphic, and igneous rocks, which differed from Precambrian to recent (Table 1). The geological succession with a concise explanation of several stratigraphic units of the study area is given in Figure 2, which is approved by the Stratigraphic committee of Pakistan [18,19]. The basement rocks of the study area are not in the typical sequence and are separated by different faults. The thrust system of the East and West limbs of the Hazara Kashmir Syntaxis (HKS) was developed due to the collision of the Indian plate and the Eurasian plate. The Kashmir Boundary Thrust (KBT) is the youngest structural feature of the projected area, which cuts the Main Boundary Thrust (MBT), Panjal Thrust (PT), and Jhelum Fault (JF) near the apex [19]. Most of the study area is incorporated by sedimentary rocks; however, metamorphic and igneous rocks are also existed [20]. For two-dimensional gravity forward modeling, the geology of the region is divided into six units, which are Murree formation of Miocene age (sandstone, mudstone, and siltstone), Kuldana formation of Eocene age (shale and marl), carbonates from Cambrian to Eocene age (shale, limestone, and sandstone), Tanol formation of Precambrian age (pelitic and psammitic metasedimentary rocks with local subordinate graphitic schist, talc schist, and marbles), Hazara formation of Precambrian age (slates and phyllites), and Salt Range formation of Precambrian age, these all information is given in Figure 1 and Table 1. Several rock samples were collected from the study area, and their bulk density was measured in the laboratory (Table 2).

Regional Tectonic Setting
Geologically Pakistan's Himalayan organic belt has several unique geological features. The central section implements a simple profile of transcontinental-continental collision boundaries ( Figure 3). Along its 2500 km long structure begins with a remarkable regular series, so there is an exceptional possibility to study variations along the length of the structure, the tectonic belt, the volcanic plutonic arc, the island arc, and the length of the thrust belt ( Figure 3). For all these purposes and more, the Himalayas contribute a good prospect of geological phenomena [21]. The Himalayan range was created by the collision between the Indian and Eurasian plates [19]. The northern boundary of the Eurasian and Indian plates is an ophiolitic belt known as the Indus-Tsangpo Suture Zone (ITSZ) [22]. The ITSZ, which extending westward, bifurcated into two major thrust faults and surrounded a sequence of rocks. The MMT and Main Karakoram Thrust (MKT). Kohistan Island Arc is present between these thrust faults. The northern suture zone that is MKT separates the intrusive high-grade metamorphic rocks of the Eurasian Plate from the Kohistan terrain, whereas the southern suture MMT separates the Kohistan Island Arc from the Meta sediments of the northern margin of the Indian Plate. The Island arc terrain is composed of a sequence of high-density rocks, while the rocks of Indian and Eurasian plates are of types of lower density [23].  [14], the blue points indicates gravity stations, and red point represent the study region.  [11,13]. The inset shows the location of the study region.

Materials and Methods
The components of the methodology section can be divided into three parts, namely, digital elevation modeling (DEMs), bouguer anomaly mapping, and two dimensional gravity data forward modeling ( Figure 4). Conventional preparation is imperative for any sort of geophysical and geographical investigations. Before starting this comprehensive investigation, at our deskwork, we studied previous research to gather information about the geological and geophysical conditions of the study area. The present study appropriates various geological maps, topographic maps, remote sensing data, gravity data, and additional sources of information that could be beneficial for surface and subsurface structure modeling. The present study utilizing integrated software's such as GM-SYS, Surfer15, and ArcGIS. The primary consideration of gravity investigation was the establishment and modeling of base stations and auxiliary stations in the research area. The D-D / profiles were usually selected based on the main geological structures in the Thandiani to Bio area, the Hazara region of Pakistan.

Digital Elevation Models (DEMs)
The digital elevation model is a digital illustration of three-dimensional information (X,Y,Z) of the continuous topography of the bare earth surface in a particular reference coordinate system [4]. The surface geological structure of the Thandiani to the Boi region in the Hazara region is considerably complicated. In this study, the analysis of bare terrain was conducted based on the DEMs, as dramatic fluctuations in elevation are signs of variations in the geological structure [5]. There are ordinarily two primary data structures, the grid structure and the TIN (triangular irregular network) structure, where DEMs data can be stored. In the grid data structure, only the elevation (z) of each node of the grid is recorded. All fluctuations of the terrain cannot be covered by the cell size of the mesh [24]. Therefore, we utilize The TIN structure to examine the surface structure changes in the study area. In the TIN data structure, the surface elevation X, Y, and Z of the terrain specific point are recorded, which describes a more realistic structure. There are numerous techniques to accumulate the DEM data, i.e., point height, stereoscopic images, active remote sensing methods (radar, lidar) [25]. This study practices the point-height method because its high-resolution images and their characteristics, such as real-time data, temporal data, and the accuracy of horizontal and vertical measurements are most popular. First, a two dimensional and three-dimensional digital elevation model based on elevation data was assembled for Thandiani and Bio region. Subsequently, two-dimensional and three-dimensional DEMs were produced for the regional zone. Terrain data were accumulated in XYZ coordinates from Google Earth, and DEM generation, DEM manipulation, DEM Interpretation, and DEM Visualization process were directed to Surfer15 and ArcGIS software (Figure 4). A DEM based on contours consists of a set of contours that are equal to the corresponding elevation value, "Mark". We use a 50-meter interval between profiles. The contours were stored in ordered coordinates and can be thought of as polylines or polyline arcs containing elevation values.

Data Acquisition
The relative gravity meters CG-5 Scintrex with a perceived resolution of 1µGal were used for gravity data acquisition ( Figure 5). Three GPS receivers were used to obtain each stations location coordinates (Table 3), including TOPCON's GRS-1 type and the 5700 L1 type of two TRIMBLERs. The GRS-1 utilizes dual frequencies, while the 5700 L1 uses one frequency to obtain location information [26]. This information helps us to understand the structural trends in our field of study. The dataset was executed in a single profile mode, and the main and auxiliary base stations were established. The main base station for gravity measurement was built in the campus of Azad Jammu University. It was connected to the National Base Station of the Church of Christ Lal Kurti Rawalpindi. The gravity of the Rawalpindi National Base Station (RNBS) was beforehand connected with the Teddington and Washington Pendulum Stations [27]. The auxiliary base stations were established in Boi. As a consequence, all measurement dimensions have been incorporated into the International Gravity Network (IGN) [28]. The gravity value of the base station is shown in Table 3. Since the study persistence was to delineate the subsurface geological structure more clearly, it was decided to conduct a single profile gravity survey of the geological structure in the design area. Meanwhile regional gravity measurements, the interval between stations are usually administered 1 km, along a single section.

Data Processing
The main purpose of data processing is to improve the observation signal-noise ratio to achieve a reliable explanation [29]. By applying the standard correction to each station's measurements, the original gravity data were reduced to produce the complete Bouguer anomaly given in equation (1). For obtained data, all changes in the Earth's gravitational field must be corrected, not due to differences in the density of the underlying rock ( Figure 6). In order to calculate and interpret the observed gravity data, it is essential to determine the density and topographic correction of the subsurface rock units [29]. Several authors examine different methods for weight reduction density determination for gravity data. Utilizing the F-H method proposed by Parasnis (Figure 6), the method is based on the minimum correlation between Bouguer gravity anomaly and Bouguer correction [30]. The reduction density of 2.4g/cm 3 is determined to be the gradient of the line (Figure 2) to correct the gravity data. Table 2 shows the absolute gravity value, free air anomaly, and the range of values for the complete Bouguer anomaly. The ultimate goal of the reduction in gravity data is to produce a complete Bouguer anomaly. Based on the above procedure, the Bouguer anomaly is given by the following equation: Where "∆g" shows Bouguer anomaly in m Gal , 'g obs " represent observed gravity value (m Gal ), "γ" shows normal gravity (m Gal ), "βh" represent free-air correction (m Gal ), "h" is the Elevation (m), "G" represents gravitational constant (6.67×10-11m 3 kg -1 s -2 ), "ρ" is the estimated density (g/cm 3 ), "2πGρh" show bouguer correction (m Gal ), and "ρT" represent terrain correction value (m Gal ).

Data Interpretation
The Bouguer anomaly is caused by the heterogeneity of the distribution of density in the subsurface. Once the gravity data is reduced to the Bouguer anomaly, the next step usually involves meshing the data to generate a map, followed by interpretation [31]. Therefore, the interpretation of gravity data is based on the Bouguer anomaly. Various geological structures are approximately linear and can be solved by two-dimensional interpretation. Different authors discussed the steps involved in interpreting gravity anomaly data [32]. In this study, two-dimensional forward modeling were used to exhibit the basements deep subsurface structure.

2D Forward Modeling
The GM-SYS profile of Geosoft Oasis Montaj 8.0.1 software was used for the interpretation of gravity data. The GM-SYS profile is a program for determining the gravity and magnetic response of geological cross-sections. Forward modeling includes building hypothetical geological models and assessing geophysical responses [32]. Gravity modeling is carried out along selected D-D / profile by utilizing the Paranis approach (Figure 6), which combines previously known structural and geological information in the study area [28]. This method was used to calculate exceptions. The observation of the gravity values along this segment were taken from the Bouguer anomaly map. According to the formation characteristics and density response, the subsurface layer is divided into different parts, namely, alluvial, Murree formation, carbonate, Hazara formation, salt range formation, and crystal crust. The density values of several rocks were perceived ( Table 2). The computational error range of our model is from 0.54 to 0.773, which is mathematically sensible. In the profile, the basement section shows different formations included in the model, their respective property densities during the modeling process. We select different colors to distinguish the formations following the information for the basement from the geological with assumed density. In the gravity modeling, the density values are adjusted until a reasonable fit is achieved between the observed calculated gravity anomalies.

Digital Elevation Models
Prior to the appearance of remote sensing and digital elevation models, field geologists were only capable of obtaining large-angle views of specific geological areas through geological maps [33]. In areas with reduced outcrops, field surveys often provide local observations that make it difficult to infer and incorporate more regional structures. With the availability of geographic maps and visualization tools such as ArcGIS Pro, Google Earth Pro, and Surfer 13, it is reasonable to check large areas [6]. Global remote sensing mosaics with sub-resolution, such as Google and Bing Web services, can broadly identify and depict geological structures. In this section, DEMs are commonly built using topographic data, predicated remote sensing techniques. The quality of DEMs is a measure of how accurate elevation is at each pixel (absolute accuracy) and how accurately the morphology is presented (relative accuracy). Various factors impersonate an important role in the quality of DEM-derived products i.e., terrain roughness, sampling density (elevation data collection methods), grid resolution or pixel size, interpolation algorithm, vertical resolution, and terrain analysis algorithm [34]. The efficiency of using a DEM for surface geological structure characterization is demonstrated here by the agreement between previous surface geological study and DEM derived data. The advantage of such a method is to give an interpretation at the scale of the slope, which is not possible when carrying out fieldwork at the regional outcrop scale. The result shows that analysis based on DEM allows identification of large scale structure and surface of failure and regional tectonic features. The above analysis shows how the use of DEM can help obtain the first interpretation of instability of surface geological and terrain structure. The surface structural trend of Thandiani and Boi area is presented in Figure 7 and Figure 8. The regional topographic surface map is presented in Figure 9.

Bouguer Anomaly Map
The Bouguer anomaly map in the projected area is compiled in the 1:500,000 range with a section interval of 2 mGal. The Bouguer anomaly map exhibits the presence of different anomaly zones in the west and north of the studied area. The gravity values of the Thandiani to the Boi region range from -124 M gal to -234 M gal (Figure 10. b). In general, high and low gravity anomalies were perceived in the northwest and southwest regions of the study area, respectively. These changes in the Bouguer anomaly map reveal the heterogeneity of the distribution of subsurface density. The anomaly is characterized by mutually alternating local maxima and minima.
The size and shape of the anomaly imply that it likely resulted from complex source structures. Due to the composition and structure of the earth's crust, gravity anomalies vary from place to place. The purpose of crustal deformation measurement is to determine the relative position or gravity changes of points on the earth's surface, to achieve geometric and physical information about crustal deformation; the concept of measurement is displacement, and the relative change of strain and gravity values [35]. The Bouguer gravity reduction density was reduced with 2.67 g/cm 3 . The contours on the gravity anomaly map usually trend northwest to the southwest, and gravity descends to the northeast. The trend of bearing force in the study area indicates that the thickness of the crust is increasing towards the northeast. According to the Bouguer anomaly contours, the area is divided into three regions, northwest, central, and southwest. In the southwest of the Thandiani region, the area is relatively stable, with fewer changes in gravity shown on the map. In this area, less deformation is observed in the rock sequence of the surface. In the central area between Thandiani and Boi, gravity values are relatively changing. In the Northwest Boi region, gravity anomalies increased, indicating changes in crustal thickness.

2D Gravity Model along the D-D′ Profile
The two-dimensional gravity modeling (along profile D-D′ selected from Thandani to the Boi area) was carried out to a depth of 54 km employing software-based on the method of Talwani and Paranis approach [32,35]. As a result the gravity model were constructed which reveals that the thickness of the crust increases to the northeast in the study area. However, in the north-east, KBT is divided between the Abbottabad formation and the Muree formation of molasses. Following the gravity model, the above part shows a graph between the observed gravity values (pink) and the calculated gravity values (green) (Figure 11). The Y-axis exhibits gravity in M gals and the X-axis in meters. The calculated model explicates that a thin layer of salt is present in the subsurface that pinches out toward the northeast of the study area. Following the result of 2D gravity modeling, it is suggested that the normal crustal thickness in the southwest part (Thandiani area) has been assumed to be 49.53 km. In the northwest part (Boi area) is 52.43 km. The thickness of the sedimentary/metasedimentary wedge in the Thandiani area is 11.48 km. In the Boi area, the thickness extends to about 14.43 km ( Figure 11).
In the gravity model, the geological structure approximates a horizontal prism with a finite length and polygonal cross-section. The model takes into account the constraints provided by deep-source gravity anomalies. According to the literature, density values have been assigned to non-survey parts of the subsurface. The density distribution is gradually altered to achieve the best fit between observing gravity and calculating gravity [28]. Given the primary constraints, the best fit between the measured and observed anomalies can be accomplished by assigning a slightly lower density to the south-west and north-east subsurface. For calculating and observing mutations in gravity data, modeling is based on three steps, namely sediment, Moho effects, and the joint effects of Moho and sediment on the 38 km thick Indian shield crystal crust [36]. The density ratio assigned to the mantle and geological bodies was 2.95 g/cm 3 compared to the average density of crystal basements. In comparison, the density assigned to the mantle was 3.30 g/cm 3 .
Concerning the shape of sedimentary/metasedimentary wedges, the most conservative solution means simplified plane geometry. Finally, the low gravity in the Thandiani to Boi area compensates with the effects of regional density discontinuity (e.g., Moho and the top of the crystal basement). The density contrasts of the subsurface vary, such as (Alluvium, Muree formation, carbonates, Hazara formation, Salt Range formation, crystalline crust, and Moho) have -0.

Conclusion
This study discusses: (1) The interpretation of uncertainty of surface geological and terrain structure by utilizing digital elevation models (DEMs), (2) the evaluation and characterization of subsurface geological structure features, basement configuration, thickness variation of crust and sedimentary/metasedimentary wedges and to comprehend the relationship of known geology with the gravity anomaly. Firstly surface topographic structure was described employing Geographic Information System (GIS), i.e., DEMs. Terrain data were accumulated in XYZ coordinates from Google Earth, DEMs generation, DEMs manipulation, DEMs interpretation, and DEMs visualization process were directed into Surfer15 and ArcGIS software. Subsequently, two and three-dimensional color embossed shaded relief DEMs were constructed, which exhibit topographic structural features, allowing them to identify a surface geological structural trend associated with the subsurface structure. Secondly, the interpretation of gravity data was performed using two-dimensional density forward modeling employing the GM-SYS software. Bouguer anomaly map depicted the crustal variations across the study area, and different formations along their densities were identified. Bouguer anomaly map is decreasing gradually towards the northeast due to an increase in crustal thickness. Less deformation is observed in the subsurface rock sequence. In the Northwest of the study area, gravity anomalies increased, indicating changes in crustal thickness. The gradual thickening of crust towards the northeast of the study area is associated with the collision of Indian and Eurasian plates. A thin layer of salt is present in the sub-surface that pinches out toward the northeast of the study area. Based on the two-dimensional computed gravity model, the thickness of the sedimentary/metasedimentary wedge in the Thandiani area is 11.48 km. In the Boi area, the thickness extends to about 14.43 km. The total thickness of the crust in the Thandiani and Boi area is 49.53 km and 52.43 km, respectively.