ISSN 1794-6190 e-ISSN 2339-3459 https://doi.org/10.15446/esrj.v28n4.114856 EARTH SCIENCES RESEARCH JOURNAL Earth Sci. Res. J. Vol. 28, No. 4 (December, 2024): 397 - 408 G EO PH Y SI C S Study of the structural characteristics of the Earth’s crust in Cao Bang province and adjacent areas using gravity data Thanh Duong Van1*, Hung Pham Nam1, Trieu Cao Dinh2, Trong Cao Dinh1, Tuan Thai Anh1, Bach Mai Xuan1, Luc Nguyen Manh1 1. Institute of Geophysics, VAST, A8/18 Hoang Quoc Viet Road, Hanoi, Vietnam 2. Institute for Applied Geophysics, VUSTA, 210 Doi Can Road, Hanoi, Vietnam *Corresponding author: dvthanh@igp.vast.vn How to cite item: Thanh, D. V., Hung, P. N., Trieu, C. D., Trong, C. D., Tuan, T. A., Xuan, B. M., & Manh, L. N. (2024). Study of the structural characteristics of the Earth’s crust in Cao Bang province and adjacent areas using gravity data. Earth Sciences Research Journal, 28(4), pending pages. https://doi.org/10.15446/esrj.v28n4.114856 Record Manuscript received: 04/06/2024 Accepted for publication: 19/11/2024 ABST R AC T Cao Bang province and adjacent areas in northern Vietnam are affected by modern tectonic fault zones, leading to frequent landslides and earthquakes. Understanding the deep structural characteristics of the Earth’s crust is crucial for advancing the study of geodynamics and geological hazards. Most previous studies on the deep structural characteris- tics of this region used only a few analysis methods, and their results were not highly reliable. Aiming to provide highly reliable information about the deep structural characteristics of this region, this study we apply a calculation process with a combination of gravity data analysis techniques and earthquake data. Specifically, the depth results obtained from the 2D modeling problem are utilized as inputs to construct an initial 3D model, thereby improving the accuracy of the depth of the crystalline basement and Moho surfaces. By integrating the results from these techniques with earth- quake data, a more accurate fault system map can be developed. Maximum horizontal gradient, tilt angle, constrained two-dimensional inversion, Euler deconvolution, 2D modeling and 3D inversion are the main techniques used. The results of the study area show that the fault system is quite complex. Large faults can be easily recognized such as F1 (Pho Day River), F2 (Yen Minh – Phu Luong), F3 (Sub-meridian fault) and F4 (Cao Bang – Tien Yen). The depth of the crystalline basement surface varies from 0 to 4.32 km, with the basement outcropping in Bac Kan province and reach- ing deep depths in the Cao Bang and That Khe basins. The depth of the Moho surface varies from 28.77 to 38.01 km. The shallowest depth of the Moho surface is located in the northern area of Cao Bang province, the southern area of Bac Kan province, the southeastern area of Lang Son province, and the Cao Bang and That Khe basins. The deepest depth value of the Moho surface is located in Bac Kan province, where the crystalline basement is outcrop. The study results have shown the effectiveness and high reliability of applying calculation processes and combining analytical methods. Keywords: Gravity; 2D modeling; 3D inversion; Fault; Crystalline basement surface; Moho surface. Palabras clave: Gravedad; Modelamiento 2D; Inversión 3D; Fallas; piso cristalino; superficie de Moho Estudio de las características estructurales de la corteza terrestre en la provincia de Cao Bang y de áreas adyacentes a través del uso de información de la gravedad La provincia de Cao Bang y las áreas adyacentes del norte de Vietnam están afectadas por zonas con fallas tectónicas modernas que le ocasionan deslizamientos de tierra y terremotos constantes. Conocer las características estructurales profundas de la corteza terrestre es determinante para el avance de la geodinámica y de los peligros geológicos. La ma- yoría de los estudios previos sobre las características estructurales profundas en esta región usaron solo algunos métodos de análisis, y sus resultados no son altamente confiables. Con el fin de aportar información altamente confiable sobre las características estructurales profundas de la región, en este estudio los autores emplean un proceso de cálculo combinado con técnicas de análisis de información de la gravedad e información de terremotos. Específicamente, los resultados de profundidad obtenidos del modelamiento del problema en dos dimensiones se utilizaron como referentes para construir un modelo inicial tridimensional, y de esta forma mejorar la precisión de la profundidad del piso cristalino y de la super- ficie Moho. Al integrar los resultados de estas técnicas con información de terremotos se puede desarrollar un mapa del sistema de fallas más preciso. Las principales técnicas usadas son el gradiente horizontal máximo, ángulo de inclinación, inversión bidimensional residual, la deconvolución de Euler, modelado 2D e inversión tridimensional. Los resultados muestran que el sistema de fallas es complejo. Grandes fallas se pueden reconocer fácilmente como la F1 (río Pho Day), F2 (Yen Minh-Phu Luong), F3 (falla submeridiana) y F4 (Cao Bang-Tien Yen). La profundidad del piso cristalino varia de 0 a 4.32 kilómetros, con un afloramiento en la provincia de Bac Kan y con mayor profundidad en las cuencas de Cao Bang y Tath Khe. La profundidad de la superficie Moho varia de 28,77 a 38,01 kms. La parte menos profunda de la superficie Moho se ubica en el norte de la provincia de Cao Bang, en el sur de la provincia de Bac Kan, en el sureste de la provincia de Lang Son, y en las cuencas de Cao Bang y That Khe. La parte más profunda de la superficie Moho se ubica en la provincia Bac Kan, donde todo el piso cristalino es un afloramiento. Los resultados del estudio han mostrado la efectividad y alta confiabilidad de aplicar procesos de cálculo combinados con métodos analíticos. RESU M EN https://doi.org/10.15446/esrj.v28n4.114856 https://doi.org/10.15446/esrj.v28n4.114856 398 Thanh Duong Van, Hung Pham Nam, Trieu Cao Dinh, Trong Cao Dinh, Tuan Thai Anh, Bach Mai Xuan, Luc Nguyen Manh 1. Introduction Study of the deep structure of the Earth’s crust is crucial for enhancing our understanding of geological structures, providing essential information for mineral exploration, and serving as a reference for further studies on geodynamics and geohazard risks. Gravity data are commonly used in large- scale studies due to its advantages. In Vietnam, there have been numerous studies on the Earth’s crust structure based on gravity data. These studies have produced depth diagrams of basic surface boundaries and fault system diagrams. However, there are still many shortcomings in the studies, such as sparse gravity measurements, limited analytical methods, and not high reliability of the results (Hai, 2003; Que, 1979, 1982, 1985; Que and Trieu, 1980, 1984; Trieu and Vuong, 1985; Trieu and Long, 2002; Trieu, 1983, 2005). In addition to using two-dimensional (2D) processing and analysis methods, numerous studies has also applied three-dimensional (3D) methods to study deep structures and determine the depth of basic Earth crust boundary surfaces. The effectiveness of the 3D inverse method has been demonstrated in determining the depth of basic boundary surfaces (Hong et al., 2016, 2017; Huong, 2006; Minh, 2014; Trung, 2014; Trung and Huong, 2005; Trung et al., 2004, 2018). However, these studies only used the 3D inverse method without combining it with other methods to build the initial model, leading to low reliability of the results. This observation is supported by the following studies, which highlight the effectiveness of combining 2D and 3D modeling approaches (Gaud, 2017; Kanthiya et al., 2019; Kebede et al., 2021; Laith, 2016; Trong, 2023). Cao Bang is a mountainous province on the Northeastern border of Vietnam, playing an important role in economic and diplomatic development, with a 314 km long border with China. The province has rich and diverse mineral resources and is a key economic area for mineral exploitation in the North of Vietnam. In particular, natural disasters such as landslides and earthquakes often occur in the Cao Bang area. One of the main reasons is that the terrain is strongly fragmented and affected by modern tectonic fault zones such as the Cao Bang – Tien Yen fault zone, Yen Minh – Phu Luong fault zone, Pho Day River fault,... In Cao Bang, a high-level earthquake with a strength equal to 5.4 occurred on November 25, 2019, causing ground collapse and severe cracking in the areas of Dinh Phong – Dam Thuy communes, Trung Khanh district, creating confusion among the population of the province, where 28 ethnic groups live. Although research on the structure of the Earth’s crust in the Northern region of Vietnam has yielded a number of published works, there have been no detailed studies of Cao Bang province and adjacent areas. Therefore, research to determine the structural characteristics of the Earth’s crust in Cao Bang province and adjacent areas is one of the important tasks, serving the assessment of earthquake severity and earthquake hazards. The distribution and structural characteristics of deep faults are the scientific basis for studying the prediction of seismic zones and earthquake activity mechanisms. In addition, the horizontal layering characteristics of the Earth’s crust also play an important role in the research and prediction of the active layer (Earthquake-generating layer) of earthquakes. This study presents a combination of 2D and 3D processing and analysis methods with gravity data. Constrained 2D inversion and Euler deconvolution methods are used as the basis to build the initial 2D model. From the depth results of the crystalline basement surface and Moho surface obtained from the 2D modeling problem, they will be used as input for the 3D inverse problem. The maximum horizontal gradient method, tilt angle method, upward continuation method, and earthquake location data are used to determine the fault system in the study area. Combining the above processing and analysis methods will help increase the reliability of the initial 3D model, thereby obtaining high-resolution results. Additionally, the combined use of different processing methods, along with data from earthquakes that have occurred, will help identify fault systems with high accuracy. The workflow diagram for calculating, processing, and analyzing gravity data for constructing the depth diagram of the basic boundary surfaces is shown in Figure 1, and constructing the fault system diagram is shown in Figure 2. Figure 1. A workflow diagram for 2-D forward and 3-D inversion gravity modeling (The blue rectangular boxes represent the input data, and the process proceeds sequentially from left to right) Figure 2. A workflow diagram for faults identification (The blue rectangular boxes represent the input data, and the process proceeds sequentially from left to right) 2. Geological and tectonic context The study area chosen is Cao Bang province and adjacent areas, limited by longitude 105°00’–107°00’; latitude 22°00’–23°30’ with geographical and X (500000–704229); Y (2432827–2600293) with the UTM perpendicular coordinate system (Figure 3). The study area is a limestone plateau mixed with dirt mountains, with an average altitude of over 0.2 km. The border area has an altitude of 0.6–1.3 km above sea level, with the highest point at 1.98 km and the lowest point at less than 0.2 km. The average altitude is 0.6–1 km above sea level, resulting in three distinct regions: the east, with many rocky mountains; the west, with dirt mountains mixed with rocky mountains; and the southwest, mostly dirt mountains with dense forests. The geology of the study area is shown in Figure 3, the 1:500,000 scale geological map established by the Vietnam General Department of Geology and Minerals (Luong et al., 1988). Cao Bang province and adjacent areas have quite complex geology, including many formations and complexes, influenced by major fault zones such as Cao Bang – Tien Yen, Yen Minh – Phu Luong, Pho Day River, etc. In the tectonic context, the Cao Bang terrane and adjacent areas are located on the edge of the Vietnam-China continental plate. The regional geological structural framework was shaped immediately after the Caledonian orogeny (Late Silurian). This structural formation system was later partly modified by the Indosinian orogenic phase (Late Triassic) and extensional deformation during 399Study of the structural characteristics of the Earth's crust in Cao Bang province and adjacent areas using gravity data the Cenozoic. The Cao Bang basin is situated at the northernmost point in the Cenozoic basin chain along the Cao Bang – Tien Yen fault zone (Figure 3). The basin has a diamond shape, nearly 20 km long in the northwest – southeast direction, with the widest point being about 6–7 km. The area is approximately 80–90 km2. Additionally, the basin is also segmented by high-order faults in the northwest – southeast and northeast – southwest directions. The Khe basin has a parallelogram shape formed in the northeast wing of the Cao Bang – Tien Yen fault along the fault direction with a length of about 6 km (Figure 3). Figure 3. Geological map of Cao Bang province and adjacent areas (Luong et al., 1988) The basement of the basin consists of Mesozoic and Paleozoic sediments (Conglomerate, sandstone, siltstone, and limestone) and igneous rocks (Triassic granite, Permian basalt, and tuff, etc.). The sediment thickness in the basin is about 0.9–1 km, including the Cao Bang formations of Eocene age (Coarse- grained sediments in the lower part) and the Na Duong formations of Eocene- Oligocene age (Fine-grained sediments in the upper part). They are distributed in uplifted marginal strips exposed in the areas of Coc Lung – Na Dong – Na Gia (Southwest range) and Na Ngan – Gia Cung (Northeast range). Late Cenozoic sediments are almost absent (Viet, 2005). 3. Materials and methods 3.1. Gravity anomaly and rock density values The Bouguer gravity anomaly data in the study area was established by the Geophysics Division in 2011 at a scale of 1:500,000 (Figure 4) and stored at the General Department of Geology and Mineral Resources (Bouguer gravity anomaly map, 2011). Some gravity measurements were compiled at a scale of 1:200,000 as part of the project “Research to identify source areas and evaluate maximum earthquakes at risk of occurring in Cao Bang province and adjacent provinces” under code VAST05.06/21-22. To build the initial model, in addition to information about boundaries and body edges, parameters about the density of rock and soil layers are also required. In this study, the soil and rock density values used are referenced from the studies of Trieu (2005) and Phong (2016). Table 1 shows the rock density values of the sedimentary crust, crystalline crust, and upper mantle. Table 1. Density parameters for soil and rock layers (Trieu, 2005; Phong, 2016) No. Layer Density (g/cm3) Average density (g/cm3) 1 Sedimentary crust 2.54–2.65 2.62 2 Crystalline crust 2.80–2.85 2.83 3 Upper mantle 3.10–3.23 3.15 Figure 4. Bouguer gravity anomaly in Cao Bang province and adjacent areas. (1) Bouguer gravity anomaly; (2) Boundary; (3) River 3.2. Terrain elevation and location of earthquakes Terrain elevation data are taken from the ASTER-GDEM source with a resolution of 30 meters (Figure 5). This is a collaborative program between NASA and the Japanese government that created the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) system. As part of this project, the ASTER Global Digital Elevation Model (GDEM) product suite was developed, providing a global-scale digital elevation model. Figure 5. Terrain elevation in Cao Bang province and adjacent areas. (1) Terrain elevation; (2) Boundary; (3) River This study utilizes additional earthquake data, specifically statistical data on earthquake locations, which plays an important role in providing more information to help identify tectonic movements and fault systems in Cao Bang province and adjacent areas. An earthquake is a shaking on the Earth’s surface resulting from the sudden release of energy in the Earth’s crust, generating seismic waves. There are many causes of earthquakes, including endogenous, exogenous, and human causes. Endogenous causes, particularly due to the movement of tectonic plates in the Earth’s crust, lead to faulting activities. Earthquakes often occur at or around active faults. Therefore, knowing the distribution of earthquakes that have occurred can provide valuable information to identify geological faults. 400 Thanh Duong Van, Hung Pham Nam, Trieu Cao Dinh, Trong Cao Dinh, Tuan Thai Anh, Bach Mai Xuan, Luc Nguyen Manh Figure 6 shows the locations of earthquakes that occurred in Cao Bang province and adjacent areas, established by the Institute of Geophysics in 2020 (Institute of Geophysics, 2020). Figure 6. Location and magnitude distribution diagram of earthquakes that occurred in Cao Bang province and adjacent areas. (1) M < 3.5; (2) M = 3.5–3.9; (3) M = 4.0–4.4; (4) M = 4.5–4.9; (5) M = 5.0–5.4; (6) Boundary; (7) River 3.3. Euler deconvolution The Euler deconvolution method was proposed by Thompson (1982) and further developed by Reid et al. (1990) to determine the boundary and depth to the source of the anomaly. we consider the function f (x,y,z) in the perpendicular coordinate system x, y, z. The function f (x,y,z) is called n degree homogeneous if: f(tx,ty,tz) = tn.f(x,y,z) (1) In this case the function f satisfies the following Equation (2): f f fx y z nf x y z ∂ ∂ ∂ + + = ∂ ∂ ∂ (2) Equation (2) is called Euler’s identity equation. Consider the set of functions f (x,y,z) of the form (Durrheim and Cooper, 1998; Reid et al., 1990; Thompson, 1982): ( , , ) N Gf x y z r = (3) Where 2 2 2r x y z= + + , N = 1, 2, 3,…; and G is a constant that does not depend on x, y and z. In practice, we take into account the deviation of the local anomalous field caused by the field with a deeper origin, which we consider constant and denote as B. We then have the Euler identity Equation (4): 0 0 0( ) ( ) ( ) ( )f f fx x y y z z N B T x y z ∂ ∂ ∂ − + − + − = − ∂ ∂ ∂ (4) Where (x0, y0, z0) is the coordinate of the magnetic source (or gravity), T is the observed field value at point (x, y, z), B is the background field value, and N is the structural index. For magnetic sources (such as uniformly magnetized spheres, magnetic dipoles, etc.), the value of the structural index varies from 0 to 3, and for gravity materials, the value of the structural index varies from 0 to 2. When calculating, one can use the actual value of the structural index. The main principle of the Euler Deconvolution method is to solve Equation (4) with the given structural index N and the need to determine 4 unknowns: x0, y0, z0, and B. 3.4. Constrained two-dimensional gravity inversion The constrained two-dimensional (2D) gravity inverse problem was applied in gravity data analysis by Guillen and Menichetti (1984) to determine the boundaries between geological structures with different density values. The structural cross-section containing the anomalous sources is divided into rectangular prisms in the x and z-axis system (the length of the object in the y-axis is considered infinite, making it a 2D case). The density contrast is constant within each prism and can vary. Using matrix notation, the gravity anomaly vector d = [di]; i = 1, 2,…, N (N is the number of data), is given by the function: d = G.m (5) Where G = [gij]; i = 1, 2,…, N; j = 1, 2,…, M is the Kernel Matrix, with gij being the contribution of the jth prism to the gravity value at the ith observation point. The model parameter vector representing the density contrast of the prism is m = [mj]; j = 1, 2,…, M; M is the number of model parameters (the number of prisms in the x and z-axis directions). The gravity component of an elementary prism is the gravitational force of a polygon in a 2D problem. In this case, the polygon has four perpendicular sides that coincide with the model grid. The value m can be determined according to the formula (Mendonca et al., 1994, 1995): m = GT [GGT + λI]-1 d (6) Where λ is the Damping Factor, I is the Unitary Matrix, and T is the Transposition Matrix. The iterative solution is carried out in 2D conditional inversion for the purpose of determining the density parameter (Hendra and Darharta, 2014; Spector and Granti, 1970). The maximum positive and minimum negative density values are of interest because they represent structural blocks. This method shows the linear residual density cross-section, helping to determine the residual density of geological objects compared to the surrounding environment. Thereby, it helps identify vertical boundaries that separate geological structural blocks with different density values (ZondGM2D, 2001). 3.5. Maximum horizontal gradient The maximum horizontal gradient method was first proposed by Blakely and Simpson (1986) to determine the boundaries between geological blocks of different densities. The horizontal gradient (G) on the Bouguer (or magnetic) gravity anomaly data grid is calculated according to the formula: 22( , ) ( , )( , ) g x y g x yG x y x y  ∂∆ ∂∆ = +   ∂ ∂    (7) With Δg(x, y) representing the Bouguer variation, the horizontal gradient value at the center point is compared with the value of the eight closest points in four directions (along the row, column, and both diagonals). If there is a maximum value, it must satisfy at least the following two inequalities: Gi-1, j < Gi,j > Gi+1,j, Gi, j-1 < Gi,j > Gi,j+1, Gi+1, j-1 < Gi,j > Gi-1,j+1, (8) Gi-1, j-1 < Gi,j > Gi+j, j+1, The maximum value of Gi,j is calculated by a quadratic polynomial function, for example: If Gi-1, j < Gi,j > Gi+1,j, then the position of the maximum value of Gi,j la xmax = 2 b a , with 1, , 1, 1, 1, 2 2 , 2 2 i j i j i j i j i jG G G G G ba d d ; d is the distance to the two nearest points. The maximum value at position xmax is calculated according to the formula: Gmax = ax2 max + bxmax +Gi,j (9) 401Study of the structural characteristics of the Earth's crust in Cao Bang province and adjacent areas using gravity data The distribution of Gmax in a linear shape reflects the existence of faults (Trieu, 2000; Van Kha et al., 2018). 3.5. Tilt angle The tilt angle method is presented in the study of Miller and Singh (1994). We can easily determine the boundary of an object from potential field data according to the arctan formula of the vertical field derivative on the horizontal field derivative: 1tan T z T h 1/222T T T h x y (10) Where: 1tan T z T h 1/222T T T h x y (11) With T x , T y , T z are the first derivatives of the potential field in the x, y and z directions; T h is the horizontal derivative, T is the potential field. Miller and Singh (1994) demonstrated the boundaries of the anomaly source using the tilt angle method. The part with positive tilt angle values represents the source location, while the part with negative tilt angle values represents the area outside the source. The part with a zero value represents the boundary of the source. 3.6. 2D gravity modeling In the case of a two-dimensional (2D) model with an n-sided polygon as an object structure, the vertical component of gravity anomaly is determined according to the formula based on the study of Talwani et al. (1959) and is simply expressed as equation (12): 1 2 m z j j g G zρ = ∆ = ∑ (12) Talwani et al. (1959) considered the case of an n-sided polygon and broke the line integral into n-contributions, each associated with a side of the polygon. They derived expressions for zj that make extensive references to trigonometric functions, as clearly stated in their paper (Talwani et al., 1959). The model involves creating a hypothetical geological model and calculating the gravity response to that Earth model. The process of solving the 2D model problem is performed using the GMSYS-2D tool, which implements the algorithm of Talwani et al. (1959) and is integrated into the Oasis Montaj software (Northwest Geophysical Associates, 2006). 3.7. 3D gravity inversion The theoretical basis of the method is based on the relationship between the Fourier transform of the gravity anomaly and the sum of the Fourier transforms of the density dividing boundary depth powers. According to Parker (1973) the Fourier transform of the gravity anomaly caused by the density dividing boundary h is determined by the expression: (Δ ) = −2 ( ) ∑ ! ∞ [ℎ ( )] (13) Where: F[] is the Fourier transform of the variables in square brackets, Δg is the gravity anomaly, G is the gravitational constant, 𝜌 is the contrast density (𝜌 = 𝜌dưới – 𝜌trên), k = (k 2 x + k 2 y) 1/2 with kx and ky are the wave numbers along the x and y axes respectively, 𝑧0 the average depth of the density separation boundary depth. In the study Oldenburg (1974) rearranged expression (13) to calculate the depth for the boundary surface from the gravity anomaly value by an iterative process and is given by: ( ) 1 2 [ ( )][ ( )] [ ( )] 2 zok n n n F g x e kF h x F h x G nπ ρ − −∞ − ∆ = − − !∑ (14) This expression allows us to determine the surface depth using an iterative back-solving process through the GMSYS-3D tool integrated in Oasis Montaj software, making the calculation and processing process more convenient (Northwest Geophysical Associates, 2006). In this procedure, we assume the average surface depth is zo and the density contrast is ρ. Gravity anomalies are determined before calculating the Fourier transform. The first term of equation (14) is then calculated by assigning h(x) = 0. The inverse Fourier transform of equation (14) will provide a first approximation of the surface, h(x). This value h(x), is then used in equation (14) to calculate a new depth value of h(x). This process is continued until a reasonable value is reached and will stop when the least squared error between the model’s gravity anomaly value and the actual gravity anomaly value drops below a given value or the allowed number of iterations is exceeded. 4. Results and discussion 4.1. Results of 2D gravity modeling According to geological documents of Cao Bang province and adjacent areas, it shows that the strike and geological faults direction are mainly in the northwest – southeast direction. Therefore, to facilitate the analysis and interpretation of geological structures in the study area, the analyzed profiles are cut perpendicular to the northwest – southeast structural direction. Based on the size of the study area as well as the map scale of the gravity anomaly and the analysis lines reflecting the regions in the study area, the author used 7 analysis profiles with a distance between points of 0.2 km. The analysis profile diagram is shown in Figure 7. Figure 7. Profile analysis diagram in the study area. (1) Profile; (2) Bouguer gravity anomaly; (3) Boundary; (4) River After obtaining the data from the 7 analysis profiles, the initial 2D model of the 7 profiles was built based on the results of constrained two-dimensional gravity inversion, Euler deconvolution, density values in Table 1, and the Bouguer gravity anomaly of analytical profiles. After completing the initial 2D model, solve the 2D forward modeling problem using the GMSYS-2D module to determine the depth to the basic boundaries, and thereby determine the appropriate density value for the geology of the study area. The results of solving the 2D forward modeling problem are shown in Figure 8, with the error between the model’s Bouguer anomaly curve and the actual Bouguer anomaly curve being ≤ 5% of the anomaly range. Blaikie et al. (2014) showed that with an error of ≤ 5% of the anomaly range, the results are highly reliable. The results of the 2D forward modeling for 7 profiles indicate that the appropriate density values for the sedimentary crust layer are 2.61–2.65 g/cm³, for the crystalline crust layer are 2.80–2.85 g/cm³, and for the upper mantle layer are 3.10–3.23 g/cm³. The crystalline basement is outcropped at profile 5. 402 Thanh Duong Van, Hung Pham Nam, Trieu Cao Dinh, Trong Cao Dinh, Tuan Thai Anh, Bach Mai Xuan, Luc Nguyen Manh Overall, the depth variations of the surfaces are quite small in profiles 1, 2, and 7, while there are significant up and down variations in profiles 3, 4, 5, and 6. These depth variations reflect important geological and tectonic characteristics of the area. They also indicate that past geological activities may have been more intense in these areas. The structural differentiation of the Earth’s crust at these locations is also manifested through different physical properties such as density and hardness. These up and down fluctuations are associated with deep fault systems, where tectonic activities may have altered the surface and depth of these layers. Another important aspect is the influence on earthquakes. The depth variations are related to the distribution of stress and energy in the Earth’s crust, impacting the formation and development of earthquake zones. In areas with significant depth changes, where the surfaces experience considerable uplift and subsidence, these regions are more prone to earthquakes due to the accumulation of stress that can easily lead to the release of energy in the form of earthquakes. Figure 8. Results of solving the 2D forward modeling problem of 7 analysis profiles, corresponding to each profile in turn from top to bottom: Actual gravity anomaly, gravity anomaly calculated from the model, error between actual gravity anomaly and calculated gravity anomaly; Results of 2D modeling problem; Geological- geophysical cross-section. (1) Euler deconvolution result; (2) Vertical boundary; (3) Horizontal boundary; (4) Constrained gravity inverse results; (5) Sedimentary crust; (6) Crystalline crust; (7) Upper mantle 403Study of the structural characteristics of the Earth's crust in Cao Bang province and adjacent areas using gravity data 4.2. Results of 3D inverse gravity modeling 4.2.1. Initial 3D model and comparison of calculation results from 2D and 3D From the results of the depth of crystalline and Moho surface, obtained after solving the 2D model problem. Proceed to create the mesh through Kriging interpolation, to obtain the depth of crystalline basement and Moho surfaces according to the area. The initial 3D model was built based on the following data: Terrain elevation, crystalline basement surface depth, Moho surface depth. Soil and rock density value data are used in combination in Table 1 and from the results of solving the 2D model problem, specifically: Sedimentary crust (2.61–2.65 g/cm3); Crystalline crust (2.80–2.85 g/cm3); Upper mantle (3.10–3.23 g/cm3). The initial 3D model is shown in Figure 9. The initial 3D model was built based on the following data: Terrain elevation, crystalline basement surface depth, Moho surface depth. Soil and rock density value data are used in combination in Table 1 and from the results of solving the 2D model problem, specifically: Sedimentary crust (2.61–2.65 g/cm3); Crystalline crust (2.80–2.85 g/cm3); Upper mantle (3.10–3.23 g/cm3). The initial 3D model is shown in Figure 9. Figure 9. The initial 3D model consists of 3 surfaces: The top layer is Terrain surface, the middle layer is Crystalline surface, the bottom layer is Moho surface After constructing the initial 3D model, the next step is to perform 3D inversion using the GMSYS-3D module to standardize the depth of the crystalline basement surface and Moho surface. This inversion process is carried out automatically. After 5 iterations with a convergence limit of 0.1 mGal, the process will stop, providing the results for the depth of the crystalline basement surfaces and Moho surfaces. Additionally, the 3D inverse calculation process will return the error result between the Bouguer gravity anomaly of the model and the actual Bouguer gravity anomaly. The results of the depth for the crystalline basement surface and Moho surface after the 3D inversion calculation are illustrated in Figure 10. According to the 2D model problem, the depth of the crystalline basement surface ranges from 0 m to 4.03 km. The location with a value of 0 m is where intrusive magmatic rocks appear, as indicated on the geological map of the study area. With the 3D modeling problem, the depth of the crystalline basement surface ranges from 0 m at the magma outcrop location to the deepest position of 4.32 km. The depth difference error of the crystalline basement surface ranges from -0.19 km to 0.12 km (Figure 10). Although this depth difference is not too large, thanks to the advantage of the 3D inverse method, the calculation and processing are carried out on grids, eliminating the need for interpolation, and at the same time standardizing the depth at positions on the grid. The results returned by the model are more consistent with the actual gravity anomaly value. For the Moho surface, the 2D model results range from 29.02 km to 37.77 km, while the 3D model results range from 28.77 km to 38.01 km. The difference between the two depth results ranges from -0.32 km to 0.38 km (Figure 11). In general, there isn’t much difference in the morphology of the Moho surface between the 2D and 3D modeling problems. However, as mentioned earlier, the 3D inverse calculation results help standardize the depth at the location of points on the grid. The final model shows a small error between the actual gravity anomaly and the calculated gravity anomaly from the model. Figure 10. Comparison of crystalline basement surface depth. (a) Calculation result from 2D model solution; (b) Calculation result from 3D model solution; (c) The crystalline basement surface depth difference between calculation results from 2D and 3D model solutions Figure 11. Comparison of Moho surface depth. (a) Calculation results from 2D model solution; (b) Calculation results from 3D model solution; (c) The Moho surface depth difference between calculation results from 2D and 3D model solutions 4.2.2. Gravity misfit between observation and calculation data Figure 12 shows the error between the actual Bouguer gravity anomaly and the feedback gravity anomaly from solving the 2D model problem, while Figure 13 shows the error between the actual Bouguer gravity anomaly and the feedback gravity anomaly from solving the 3D modeling problem. From the error results, we can see that when solving the 2D model problem, the gravity anomaly error value ranges from -1.6 to 1.6 mGal. In contrast, when solving the 3D model problem, the gravity anomaly error value ranges from -0.5 to 0.65 mGal. Therefore, it is evident that the results of the 3D inverse problem are more reliable than the results of solving the 2D model problem (Skalbeck, 2014). 404 Thanh Duong Van, Hung Pham Nam, Trieu Cao Dinh, Trong Cao Dinh, Tuan Thai Anh, Bach Mai Xuan, Luc Nguyen Manh Figure 12. Gravity misfit between observation and calculation data from the 2D model. (1) 2D anomaly error; (2) Boundary; (3) River Figure 13. Gravity misfit between observation and calculation data from the 3D model. (1) 3D anomaly error; (2) Boundary; (3) River After obtaining the results of the depth to the crystalline basement surface and the Moho surface from both the 2D and 3D model problems, it is evident that the results of the 3D model problem show higher reliability than those of the 2D model problem. Additionally, the 3D model results help standardize the depth to the crystalline basement and Moho surfaces. Figure 10 and Figure 11 illustrate the difference in the depth of the crystalline basement surface and the Moho surface between the results obtained from solving the 2D model and the 3D model problems. 4.2.3. Crystalline basement surface The crystalline basement surface is often considered the boundary between the poorly consolidated sedimentary rock layer in the upper part of the structural section and the underlying granite layer. This boundary often experiences significant changes in rock density and wave velocity. The study results show that the crystalline basement surface in the study area exhibits a complex structural morphology, ranging from basement outcrop in some locations to a depth of 4.32 km (Figure 14). The results reveal three distinct value ranges: the outcrop location, the location with medium depth, and the location with deep depth. According to Figure 14, in Cao Bang province, the depth of the crystalline basement surface varies from 2.01 to 4.31 km from west to east. The change in depth clearly reflects the influence of the Cao Bang – Tien Yen fault zone. The depth of the crystalline basement surface in Cao Bang reaches a high value, indicating a concave structure in the Cao Bang basin and its surroundings. In Lang Son province, the depth of the crystalline basement varies from 1.51 to 3.42 km, with the higher depth values representing the concave structure at That Khe basin. In Bac Kan province, the crystalline basement location ranges from an exposed basement with a value of 0 km to a depth of 2.92 km. In Ha Giang province, convex and concave structures are intertwined, but the difference in depth is not significant, ranging from 1.51 to 3.01 km. As for Tuyen Quang province, there is a small change in the depth of the crystalline basement, ranging from 1.41 to 2.51 m. Figure 14. Diagram of crystalline surface depth in Cao Bang province and adjacent areas. (1) Depth of crystalline basement surface; (2) Boundary; (3) River Comparing the results of crystalline basement depth in Cao Bang province and adjacent areas with previous studies such as Hai (2003) and Phong (2016) shows that most of the results are similar, with the difference in depth not being too significant. In particular, Hai’s study (2003) does not show the location of crystalline basement outcrop, although the geological map of Vietnam at a scale of 1:500,000 (Luong et al., 1988) indicates exposed magmatic rocks (Figure 3). In the eastern part of Cao Bang province, a concave structure is evident with a depth approximately 0.3 to 0.4 km greater than that reported by Phong (2016). 4.2.4. Moho surface The Moho surface is the boundary between the Earth’s crust and mantle. The analysis results in Figure 15 show details of the change in Moho surface depth in Cao Bang province and adjacent areas, with the depth ranging from 28.77 km to 38.01 km for the entire study area. The structures tend to stretch mainly in the northwest – southeast direction. In Cao Bang province, the Moho surface depth varies from 28.77 km to 36 km. The northern part of Cao Bang province has the shallowest and moderately shallow Moho basement depth in the Cao Bang basin at 32.6 km. In Lang Son province, the depth of the Moho surface varies from 31 km to 35 km. Around That Khe basin, the depth is 32.75 km. In Bac Kan province, the depth of the Moho surface varies from 32 km to 36.51 km. The location with the deepest depth value is where the crystalline basement surface is outcrop, with a depth of 36.51 km. The depth of the Moho surface in Tuyen Quang province varies from 32.01 km to 37.01 km, with the majority being deep depths. In Ha Giang province, which has the highest topographic elevation in the study area, the depth of the Moho surface varies from 32.51 km to 38.01 km. This is also the location with the deepest Moho surface depth of 38.01 km in the study area. Comparing the Moho surface depth of this study with previous studies such as Hai (2003) and Phong (2016), the results show similarities in Moho surface depth with the studies by Hai (2003) and Phong (2016). The deep depth values are all located in the north of Ha Giang province, with the depth gradually decreasing towards the south of the study area, including Bac Kan and Lang Son provinces. In this study, the deepest depth of the Moho surface is 38.01 km compared to the deepest depth of studies by Hai (2003) and Phong (2016), which is 37.2 km, showing a difference of about 0.8 km. The shallowest Moho surface depth in this study is 28.77 km, while the shallowest 405Study of the structural characteristics of the Earth's crust in Cao Bang province and adjacent areas using gravity data Moho surface depth in the studies by Hai (2003) and Phong (2016) is 30 km, showing a difference of about 0.12 km. In addition to these differences, the results of the study provide more detailed information about the depth of the Moho surface in the study area, which can be clearly seen as interwoven convex and concave structural blocks (Figure 15). This difference in results may be due to previous studies covering a larger area, along with limitations in data and analysis methods, resulting in a significantly reduced detailed reflection of study results for a smaller area. Figure 15. Moho surface depth diagram in Cao Bang province and adjacent areas. (1) Earthquake seismic recording station; (2) Moho surface depth; (3) Boundary; (4) River On the other hand, in determining the depth to the basic boundary surfaces, data obtained from earthquake seismic recording station networks reflect this very well at the station locations. However, because the number of earthquake seismic station networks is quite sparse, the results showing the depth to the basic boundaries for the entire area only reflect the average. In the study area of Cao Bang province and adjacent areas, there are two points of Moho surface depth results through the use of the Receiver Function analysis method of seismic data from the study by Duong et al. (2013). These two Moho surface depth points are shown in Figure 15, with symbols CBVB (Cao Bang) indicating a depth of 32.5 km and HGVB (Ha Giang) indicating a depth of 32.6 km. Comparing these results with the Moho surface depth results of this study, at the CBVB point, the depth is 32.563 km, resulting in a difference in depth compared to the results of the study by Duong et al. (2013) of about 0.063 km. Similarly, at the HGVB point, the depth in this study is 32.513 km, resulting in a depth difference of 0.087 km. With the very small depth differences that have been verified between this study and the earthquake seismic recording station network data, it has shown high reliability in the Moho surface depth calculation results of this study. 4.3. Identify the geological fault system Geological faults have become the main focus of study for geologists and geophysicists studying the deep structure of the Earth’s crust. Additional information is added to geological and geophysical documents, combining various analysis and processing methods, based on the analysis of gravity data along with the location data of earthquakes that have occurred in the study area. Some criteria for identifying faults include: • Local anomalous bands (negative or positive) forming elongated linear chains. • Sudden changes in gravity anomaly contour lines. • Locations with maximum horizontal gradient values of the gravity field extending in a straight line or crescent line with a large radius. • Large, deep faults often appearing at locations where maximum horizontal gradient values exist according to upward continuation levels. • In addition, with the tilt angle method, the fault location will pass through the position where the value is 0. • With the data on earthquake epicenters that have occurred in the study area, the locations of these events or their surroundings may indicate the presence of faults or tectonic plate boundaries. Figure 16 shows a diagram of the fault system in Cao Bang province and adjacent areas, determined based on geological data, earthquake data, results of gravity data processing, and analysis methods such as the Maximum Horizontal Gradient method with upward continuation levels of 10 km and 20 km, and the tilt angle method. From the results obtained through processing and analyzing gravity anomaly data combined with the location data of earthquakes in the study area, 17 large and small faults have been identified and labeled from F1 to F17 (Figure 16). Among these, there are four deep and large faults: F1, F2, F3, and F4 (F4.1, F4.2). Clearly identifiable signs of these faults are shown through the maximum horizontal gradient values at upward continuation of 10 km and 20 km, as well as the contour lines with a tilt angle value of zero. Additionally, earthquakes have occurred in that location. The remaining faults, labeled from F5 to F17, are identified as small faults. Most earthquakes in this study area are concentrated along major faults, while smaller and isolated earthquakes occur along minor and shallow faults. It is important to note that the epicenter locations of recorded earthquakes represent surface positions, not the hypocenters at the depths where the events occurred. This is only one of the indicators used to identify faults or tectonic boundaries. Therefore, some epicenters may be offset from the determined fault locations. According to the study results in Cao Bang province and adjacent areas (Figure 16), the sub-meridian arc fault zones in the north are F1 (Pho Day River), F2 (Yen Minh – Phu Luong), and F3, which tend to lean towards the west. F4.1 and F4.2 (Cao Bang – Tien Yen) and F9 tend to tilt clearly to the Southwest. Faults F7, F10, F12, F14, F15, and F17 all develop in the northwest – southeast direction. F16 is an arc-shaped, sub-parallel fault that develops from northwest to southeast. Additionally, the remaining faults, such as F5, F6, F8, F11, and F13, develop in the northeast – southwest direction. Figure 16. Diagram of the fault system in Cao Bang province and adjacent areas. (1) M < 3.5; (2) M = 3.5–3.9; (3) M = 4.0–4.4; (4) M = 4.5–4.9; (5) M = 5.0–5.4; (6) Provincial boundary; (7) River; (8) Gmax value at 0 km; (9) Gmax value at upward level 10 km; (10) Gmax value at upward level 20 km; (11) Value of tilt angle method; (12) Fault Comparing with fault documents of previous studies such as Luong et al. (1988), Toan et al. (2010), and Phong (2016), the results of determining the fault system in this study are highly similar to those of the aforementioned studies. The clearest similarity observed is that the large, deep faults F1, F2, and F4 (F4.1, F4.2) correspond to the Pho Day River fault, Yen Minh – Phu Luong fault, and Cao Bang – Tien Yen fault. Additionally, fault F3 also shows similarities with the results of fault identification in Toan et al. (2010). Small faults such as F5, F6, F7, F9, F10, F13, F15, F16, and F17 were also identified in the studies by Luong et al. (1988) and Phong (2016). However, some small 406 Thanh Duong Van, Hung Pham Nam, Trieu Cao Dinh, Trong Cao Dinh, Tuan Thai Anh, Bach Mai Xuan, Luc Nguyen Manh faults identified in this study, such as F8, F11, F12, and F14, did not appear in the above research. The manifestation of these faults through the analysis results of the methods is quite clear. Their absence in the above studies may be due to the fact that the data used for analytical processing of the above works are not detailed, and the analytical methods used are limited. Furthermore, research covering a large area may lead to small faults being overlooked and undetected. 5. Conclusions This study has rigorously and logically applied the process of analyzing gravity data from 2D to 3D, using the combined results of methods such as constrained gravity inversion and Euler deconvolution, to build an initial 2D model with high suitability. The calculated results from the 2D model were then used as input for building the initial 3D model. The initial 3D model obtained has parameters, such as density and depth, that are closer to the actual values. The final results obtained after solving the 3D model problem demonstrate high reliability, with small error values between the calculated and actual measured gravity anomalies, ranging from -0.5 to 0.65 mGal. At the same time, it shows the advantages of time, fast processing speed, and high accuracy of the 3D inverse method compared to the 2D model method. Through the combined use of a number of processing and analysis methods such as upward continuation, maximum horizontal gradient, tilt angle, and location data of earthquakes that have occurred, it has helped to identify the fault system in Cao Bang province and adjacent areas. Shows high efficiency in using a combination of the above methods and data to identify the fault system. The results showed that the fault system in the study area is quite complex, the main direction of development is northwest – southeast. A total of 17 faults, both large and small, have been identified, large faults are easily recognizable such as F1 (Pho Day River), F2 (Yen Minh – Phu Luong), F3 (Arc, sub-meridian faults) and F4 (Cao Bang – Tien Yen). In particular, this study identified additional faults F8, F11, F12, and F14, which were not detected in previous research. In addition, there are a number of small arc and sub-parallel faults, along with a number of other small faults developing in the northeast – southwest direction. The depth of the crystalline basement surface in Cao Bang province and adjacent areas varies from 0 km to 4.32 km. The location of the basement outcrop is in Bac Kan province, the location has deepest depth values in the Cao Bang, That Khe basins and the location around basins. The depth of the Moho surface in Cao Bang province and adjacent areas varies from 28.77 km to 38.01 km. Locations with shallow Moho surface depth values range from 28.77 km to 30 km, appearing in the northern region of Cao Bang province, the southern region of Bac Kan province, and the southeastern region of Lang Son province. In the Cao Bang and That Khe basins, the depth of the Moho surface is relatively shallow, ranging from 32.51 km to 33.01 km. The deepest Moho surface depth is found in Ha Giang province, at 38.01 km. The uplift and subsidence of the crystalline basement and Moho surfaces reflect significant geological and tectonic features of the study area. Locations where alternating convex and concave structures are observed indicate variations in surface depth, resulting in higher stress accumulation compared to other areas. When combined with fault zones intersecting these locations, energy release occurs, leading to earthquakes. From the fault map and the depth profiles of the crystalline basement and Moho surfaces derived from the study, areas with a high potential for seismic activity have been identified, such as the deep fault zones F1, F2, F3, and F4. These faults show significant changes in surface depth where they pass through. Similarly, the regions traversed by faults F8, F9, F10, F11, F15, F16, and F17, which have already experienced some earthquakes, also exhibit notable variations in the depth of the crystalline basement and Moho. All this information provides deeper insights into the Earth’s structure and its dynamic processes, thereby supporting the assessment of geological risks, resource planning, and the study of regional geodynamics. The results of this study were specifically evaluated and compared with documents and results from previous studies, demonstrating the high effectiveness of our approach. Applying this analysis process to other study areas shows promise and is feasible. References Aitken, A. R. A., Salmon, M. L., & Kennett, B. L. N. (2013). Australia’s Moho: A test of the usefulness of gravity modeling for the determination of Moho depth. Tectonophysics, 609, 468–479. Alessandra, B. (2022). Moho depths for Antarctica Region by the inversion of ground-based gravity data. Geophys. J. Int, 231, 1404–1420. Bach, L. D. (2017). Fault in the Northeast. Bilim, F. (2017). Investigating Moho depth, Curie Point, and heat flow variations of the Yozgat Batholith and its surrounding area, north central Ana- tolia,Turkey, using gravity and magnetic anomalies. Turkish journal of earth sciences, 26, 410–420. Blaikie, T. N., Ailleres, L., Betts, P. G. & Cas, R. A. F. (2014). Interpreting sub- surface volcanic structures using geologically constrained 3-D gravity inversions: Examples of maar-diatremes, Newer Volcanics Province, southeastern Australia. J. Geophys. Res. Solid Earth, 119. Doi:10.1002/ 2013JB010751. Blakely, R. J. & Simpson, R. W. (1986). Approximating edges of source bodies from magnetic or gravity anomalies. Geophysics, 51, 1494–1498. Blakely, R. J. (1995). Potential Theory in Gravity and magnetic application. Cam- bridge: Cambridge University Press. Böhme, M., Jérôme, P., Simon, S., Hung, N. V., Quang, D. D. & Tran, D. N. (2011). The Cenozoic on-shore basins of Northern Vietnam: Biostratigraphy, vertebrate and invertebrate faunas. Journal of Asian Earth Sciences, 40(2), 672–687. Bouguer gravity anomaly map (2011). Cao Bang and adjacent newspaper, Publi- shed at the Federation of Geophysics, Hanoi. Chinh, V. V. & Thinh, N. H. (2001). Research on tectonic faults in Ha Giang Town and adjacent areas. Journal of Earth Sciences, 23(4), 378–389. Dung, L. V. & Trieu, C. D. (2012). Representation of lithospheric activity in the Northwest region. Geological Journal, series A, 331, 111–123. Dung, T. T., Kha, T. V. & Ha, N. T. H. (2012). Application of new documents on seismicity, seabed depth and gravity to build a structural model of the Earth’s crust in the East Sea and adjacent areas. Dung, T. T., Kulinich, R. G., Sang, N. V., Que, B. C., Dai, N. B., Dung, N. K., Duong, T. T. & Lap, T. T. (2019). Improving Accuracy of Altimeter-deri- ved Marine Gravity Anomalies for Geological Structure Research in the Vietnam South-Central Continental Shelf and Adjacent Areas. Russian Journal of Pacific Geology, 13(4), 364–374. Durrheim, R. J. & Cooper, G. R. (1998). Euldep: A Program for the Euler decon- volution of magnetic and gravity data. Computers & Geosciences, 24(6), 545–550. Ebadi, S., Safari, A., Barzaghi, R. & Bahroudi, A. (2020). A recovered Moho mo- del by integrated inversion of gravity and seismic depths in Iran. Heli- yon, 6(3), e03636. Gaud, P., Gerry, C., Hannah, K., Mohamed, G. & Douglas, P. (2017). Public domain satellite gravity inversion offshore Somalia combining laye- red-Earth and voxel based modeling. First Break. Gozzard, S., Kusznir, N., Franke, D., Cullen, A., Reemst, P. & Henstra, G. (2018). South China Sea crustal thickness and oceanic lithosphere distribution from satellite gravity inversion. Petroleum Geoscience, 25(1), 112–128. Grigoriadis, V. N., Tziavos, I. N., Tsokas, G. N. & Stampolidis, A. (2015). Gravity data inversion for Moho depth modeling in the Hellenic area. Pure and Applied Geophysics, 173(4), 1223–1241. Guillen, A. & Menichetti, V. (1984). Gravity and Magnetic Inversion with Mini- mization of a Specific Functional. Geophysics, 49(8), 1354–1360. Haase, C., Ebbing, J. & Funck, T. (2016). A 3D regional crustal model of the NE Atlantic based on seismic and gravity data. Geological Society, London, Special Publications, 447(1), 233–247. 407Study of the structural characteristics of the Earth's crust in Cao Bang province and adjacent areas using gravity data Hai, D. T. (2003). Research on some deep structural features of the Earth’s crust and seismic tectonic zoning in Northern Vietnam. Ph.D. Thesis, Institute of Geophysics, Hanoi. Hendra, G. & Darharta, D. (2014). Constrained Two-Dimensional Inversion of Gravity Data. ITB Journal Publisher, ISSN: 2337–5760, J. Math. Fund. Sci., 46(1), 1–13. Hong, P. T., Trung, N. N. & Nam, B. V. (2016). Determining the depth of the Con- rad basement surface in the northern area of the Red River basin and adjacent areas based on gravity data analysis. Journal of Mining Science and Technology – Geology, 57, 1–13. Hong, P. T., Trung, N. N., Nam, B. V. & Phuong, D. M. (2017). Deep structure of the Earth’s crust in the Hanoi sag area and surrounding areas based on gravity anomaly field analysis. Science Magazine Mining and Geology Engineering, 58(5), 325–334. Hung, P. N. & Dung, L. V. (2011). Deep geological structure in Hanoi and ad- jacent areas based on gravity data analysis. Journal of Earth Sciences, 185–190. Hung, P. N., Trong, C.D., Dung, L. V., Tuan, T. A., Bach, M. X. & Duong, T. S. (2019). Study on structure of the Earth’s crust in Thua Thien-Hue pro- vince and adjacent areas by using gravity and magnetic data in combina- tion. Vietnam Journal of Marine Science and Technology, 19(4), 517–526. Huong, N. T. (2006). Geological effectiveness of gravity method to study the deep structure of the Earth’s crust and tectonics of sedimentary basins in the continental shelf of Vietnam. Ph.D. Thesis, Hanoi. Institute of Geophysics (2020). Updated earthquake notices and list of earth- quakes in Vietnam until the end of 2020, Documents archived at Insti- tute of Geophysics, VAST. Kanthiya, S., Mangkhemthong, N. & Morley, C. K. (2019). Structural interpre- tation of Mae Suai Basin, Chiang Rai Province, based on gravity data analysis and modeling. Heliyon, 5(2), e01232. Kebede, H., Alemu, A., Nedaw, D. & Fisseha, S. (2021). Depth estimates of ano- malous subsurface sources using 2D/3D modeling of potential field data: implications for groundwater dynamics in the Ziway-Shala Lakes Basin, Central Main Ethiopian Rift. Heliyon, 7(4), e06843. Khoan, P. & Ixaev, A. N. (1971). Some basic features of the deep structure of the Earth’s crust in Northern Vietnam. Journal of Biology – Geology, 9(1–2), 55–63. Laith, S. & Ali, M. A. R. (2018). Constrained 3D Gravity Inversion for the Me- sopotamia Basin – A Case Study. International Journal of Science and Research (IJSR), 7(4), 2319–7064. Luong, T. D., Bao, N. X. et al. (1988). Geological map at scale 1:500,000. Published at the Vietnam General Department of Geology and Minerals, Hanoi. Mendonca, C. A. & Silva, J. B. (1994). The Equyvalent Data Concept Applied to the Interpolation of Potential Field Data. Geophysics, 59(5), 722–732. Mendonca, C. A. & Silva, J. B. (1995). Interpolation of Potential Field Data by Equyvalent Layer and Minimum Curvature: A Comparative Analysis. Geophysics, 60(2), 399–407. Milano, M., Kelemework, Y., La Manna, M., Fedi, M., Montanari, D. & Iorio, M. (2020). Crustal structure of Sicily from modeling of gravity and magnetic anomalies. Scientific Reports. Miller, H. G. & Singh, V. (1994). Potential field tilt a new concept for location of potential field sources. Journal of Applied Geophysics, 32, 213–217 Minh, N. Q. (2014). Research on the characteristics of pre-Cenozoic basement structures in the Truong Sa and Tu Chinh – Vung May archipelagos. Mas- ter’s Thesis. Northwest Geophysical Associates, Inc. (2006). GM-SYS-2D Gravity and Mag- netic Modeling Software User’s Guide Revision 4.10.03. Northwest Geophysical Associates, Inc. (2006). GMSYS-3DTM 3D Gravity and Magnetic Modeling for Oasis montaj™ User’s Guide release 1.3. Oldenburg, D. W. (1974). The inversion and interpretation of gravity anomalies. Geophysics, 39, 526–536. Parker, L. (1973). The rapid calculation of potential field anomalies. Geophysical Journal of the Royal Astronomical Society, 31, 447–455. Phong, L. H. (2016). Model of Earth’s crust structure in Northern Vietnam based on seismic and gravity data. Ph.D. Thesis, Vietnam Academy of Science and Technology. Prasanna, H. M. I., Chen, W. & İz, H. B. (2013). High resolution local Moho de- termination using gravity inversion: A case study in Sri Lanka. Journal of Asian Earth Sciences, 74, 62–70 Quang, P. T. (2016). Research on the structural characteristics of the earth’s crust in the Song Da hydropower terrace area using gravity method. Master’s Thesis. Que, B. C. (1979). Method for synthesizing exploration geophysical data using multiple correlation analysis. Research works of the Institute of Earth Sciences 1977–1978, Physics Volume Earth, 80–96. Que, B. C. (1982). Geological effectiveness of the gravity method in studying the deep structure of the Earth’s crust in the territory of Vietnam. Journal of Earth Sciences, 4, 23–31. Que, B. C. (1985). Some deep structural and tectonic features of the southern part of Vietnam according to geophysical documents. Collection of works on Geophysics 1984, 4, 156–165. Que, B.C. & Trieu, C. D. (1980). Learn more about some structural features within the Earth’s crust in Northern Vietnam. Results of Geophysics research in 1979, 239–245. Que, B. C. & Trieu, C. D. (1984). Using the synthesis of geophysical data to study the deep structure of the Earth’s crust in the territory of Vietnam. Scien- tific Notice 1, 79–81. Quyet, P. Q. (1985). Application of gravity exploration method to study geological structures in the Mekong Delta. Associate Doctor of Science thesis, Ha- noi University of Mining and Geology. Reid, A. B., Allsop, J. M., Granser, H., Millet, A. J. & Somerton, I. W. (1990). Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics, 55, 180–191. Saada, S. A., Eldosouky, A. M., Abdelrahman, K., Al-Otaibi, N., Ibrahim, E. & Ibrahim, A. (2021). New insights into the contribution of gravity data for mapping the lithospheric architecture. Journal of King Saud Univer- sity – Science, 33(3), 101400. Salem, A., Green, C., Campbell, S. D. & Fairhead, J. (2012). A Practical Approach to 3D Inversion of Pseudo-gravity for Depth to Basement Mapping – A Test Using the Bishop Model, 74th EAGE Conference and Exhibition Incorporating EUROPEC. Salem, A., Green, C., Fairhead, D., & Aboud, E. (2012). Mapping basement relief of Abu Gharadig Basin, Western Desert of Egypt using 3D inversion of pseudo-gravity data. ASEG Extended Abstracts, 2012(1), 1–4. https:// doi.org/10.1071/ASEG2012ab385. Sanchez Rojas, J. (2012). New Bouguer Gravity Maps of Venezuela: Represen- tation and Analysis of Free-Air and Bouguer Anomalies with Empha- sis on Spectral Analyses and Elastic Thickness. International Journal of Geophysics, 1–15. Skalbeck John, D., Koski Adrian, J. & Peterson Matthew, T. (2014). Estimation of Precambrian basement topography in Central and Southeastern Wis- consin from 3D modeling of gravity and aeromagnetic data. Journal of Applied Geophysics, 106, 187–195. doi:10.1016/j.jappgeo.2014.04.019. Spector, A. & Granti, F. S. (1970). Statistic model for interpreting aeromagnetic data. Geophysics 35(2), 293–302. Steffen, R., Steffen, H. & Jentzsch, G. (2011). A three-dimensional Moho depth model for the Tien Shan from EGM2008 gravity data. Tectonics, 30(5). doi:10.1029/2011tc002886. Talwani, M., Worzel, J. L. & Landisman, M. (1959). Rapid gravity computations for two dimensional bodies with application to the Mendocino subma- rine fracture zone. Journal of Geophysical Research, 64, 49–59. https://doi.org/10.1071/ASEG2012ab385 https://doi.org/10.1071/ASEG2012ab385 408 Thanh Duong Van, Hung Pham Nam, Trieu Cao Dinh, Trong Cao Dinh, Tuan Thai Anh, Bach Mai Xuan, Luc Nguyen Manh Thai, H. M. (1990). Determining crystalline basement in the Red River basin according to gravity data using the adjusted selection method. Journal of Earth Sciences, 3, 20–25. Thompson, D. T. (1982). EULDPH-A new technique for making computer-assis- ted depth estimates from magnetic data. Geophysics, 47, 31–37. Toan, D. V. (1988). Experimental research applying a number of new methods to process and analyze gravity data in the Hanoi lowland area. Ph.D. Thesis, Institute of Geophysics, Hanoi. Toan, D. V. (1993). Ability to use two-dimensional modeling method of gravity anomalies to study geological structures. Journal of Earth Sciences, 15, 92–96. Toan, D. V. et al. (2010). Report on the results of studying the deep structure of the Earth’s crust in Northern Vietnam using deep seismic and tellurian mag- netic detection to improve the reliability of geological disaster forecasts, code: KC.08.06/06-10. Archives of the Institute of Geology, Hanoi, 256 pp. Trieu, C. D. (1983). Deep structure of North Vietnam’s territory. Ph.D. Thesis, Vietnam Academy of Sciences, Hanoi, 140 pp. Trieu, C. D. (2000). Gravity and gravity probe methods. Science and Technology Publishing House, Hanoi, 276 pp. Trieu, C. D. (2005). School of Geophysics and lithospheric structure of Vietnam’s territory. Science and Technology Publishing House, Hanoi, 330 pp. Trieu, C. D. & Vuong, H. V. (1985). Fault research method based on the Earth’s crust structure model. Collection of works on Geophysics in 1984, 4, 185–197. Trieu, C. D. & Long, P. H. (2002). Fault tectonics in Vietnam’s territory. Science and Technology Publishing House, Hanoi, 200 pp. Trinh, D. D., Phuong, B. T. A. & Liet, D. V. (2009). Determining the Moho surface in Southern Vietnam. Journal of Earth Sciences, 31(4), 403–409. Trong, C. D. (2022). Report on the results of research to identify source areas and evaluate maximum earthquakes at risk of occurring in Cao Bang province and adjacent provinces, code VAST05.06 /21-22. Archives of the Institute of Geophysics. Trong, C. D., Hung, P. N., Thanh, D. V., Luc, N. M., Bach, M. X., Trieu, C. D. & Hung, L. V. (2023). Chromite Ore Modeling Based on Detail Gravi- ty Method in Pursat, Combodia. Advances in Geospatial Technology in Mining and Earth Science, Selected Papers of the 2nd International Conference on Geo-spatial Tecnologies and Earth Resources 2022, pp. 259–275. Tru, V. et al. (2013). Synthesizing the results of geological – geophysical research and proposing a research plan for pre-tertiary sedimentary basins to serve oil and gas exploration in Vietnam. Summary report of industry-level scientific research topic, Vietnam Petroleum Institute. Trung, N. N. (2014). Determining the depth of sedimentary basement in the Bac Bo basin area according to 3D direct inverse analysis of satellite gravity data. Journal of Earth Sciences, 315–320 Trung, N. N. & Huong, N. T. T. (2005). The structure of the Earth’s crust in the East Sea region according to satellite gravity anomaly data and deep seismicity, Collection of reports from the Petroleum Institute’s Scien- ce and Technology Conference: 25 years of construction and growth, 336–354. Trung, N. N., Hong, P. T., Van Nam, B., Huong, N. T. T. & Lap, T. T. (2018). Moho depth of the northern Vietnam and Gulf of Tonkin from 3D inverse interpretation of gravity anomaly data. Journal of Geophysics and Engi- neering, 15(4), 1651–1662. Trung, N. N., Lee, S. M. & Que, B. C. (2004). Satellite Gravity Anomalies and Their Correlation with the Major Tectonic Features in the South China Sea. Gondwana Research, 7(2), 407–424. Van Duong, N., Bor Shouh, H., Tu Son, L., Van Toan, D., Lupei, Z. & Kuo Liang, W. (2013). Constraints on the crustal structure of northern Vietnam based on analysis of teleseismic converted waves. Tectonophysics, 601, 87–97. Van Kha, T., Van Vuong, H., Thanh, D. D., Hung, D. Q. & Anh, L. D. (2018). Im- proving a maximum horizontal gradient algorithm to determine geolo- gical body boundaries and fault systems based on gravity data. Journal of Applied Geophysics, 152, 161–166. Viet, L. T. (2005). Formation of some Cenozoic depressions in the Northeastern region, Vietnam. Geological Journal. Wysocka, A. (2009). Sedimentary Environments of the Neogene Basins Asso- ciated with the Cao Bang – Tien Yen fault, NE Vietnam. Acta Geologica Polonica, 59(1), 45–69. Zhang, J., Yang, G., Tan, H., Wu, G. & Wang, J. (2021). Mapping the Moho depth and ocean-continent transition in the South China Sea using gravity inversion. Journal of Asian Earth Sciences, 218, 104864. Zhao, G., Liu, J., Chen, B., Kaban, M. K. & Zheng, X. (2020). Moho beneath Ti- bet based on a joint analysis of gravity and seismic data. Geochemistry, Geophysics, Geosystems. Doi:10.1029/2019GC008849. ZondGM2D, (2001). Program for 2D interpretatino of magnetic and gravity data. Saint-Petersburg.