The Gonghe Basin, which is located in the northeast of Qinghai-Tibet Plateau in Qinghai Province, has characteristics of strong tectonic activity, relatively complex stress state and uneven crustal structure and has become the third largest basin in Qinghai Province (Zhang et al., 2018). The Gonghe basin is surrounded by several faults(Fig. 1). As a typical Cenozoic faulted basin bounded by Ela Mountain and Yellow River valley in the south and Qinghai Nan Mountain in the north, it is tectonically controlled by Kunlun Orogen in the south and Qilian Orogen in the north (Fang et al., 2005). Surrounded by Qinghainan fault in the north, Duohemao Fault in the east and East Kunlun Fault in the southwest, Gonghe Basin is about 280km in length, 95km in width, 2.1×104 km2 in area and distributed in a diamond shape (Zhang et al., 2006).
The upper base of the Gonghe basin is mainly composed of the Early-Middle Triassic Longwuhe Formation, Middle Triassic Gulangdi Formation and Middle-Late Triassic granites. The middle and lower basement is presumed to be Paleozoic-Proterozoic metamorphic rock series(Feng et al., 2003). The sediments above the basement are mainly Neogene and Quaternary, and the surface layer is Quaternary. As indicated in regional geological results, these has been no magmatic activities around Gonghe Basin since Cenozoic (Zhang et al., 2007;Zhang et al., 2020). The igneous and metamorphic rock sequences of the East Kunlun and West Qinling Orogen developed around the basin is composed of the middle and late Triassic metamorphic sedimentary rocks and Indosinian granites and granodiorites (Yan et al., 2015). The Cenozoic deposits in the basin generally show a trend of thin in the east and thick in the west. And the thickness of Cenozoic deposits is generally about 500m to 1500m in the eastern part of Guide area, about 1500m in the central and eastern parts of Qiabuqia area and can reach 6000m in the west (Tang et al., 2020).
There are many hot springs exposed in Gonghe Basin with a relatively obvious concentrated distribution along the fault zone (Feng et al., 2018). In addition to hydrothermal geothermal resources, Gonghe Basin is also one of the most powerful areas for exploration and development of hot dry rock geothermal resources in China. The research shows that the average geothermal gradient in Gonghe Basin is twice as big as the normal geothermal gradient (Liu., 2017). According to the report of China Geological Survey and Department of Natural Resources of Qinghai Province in August 2017, the bottom hole temperature of five wells with depths between 3000m and 3705m in basin is between 180℃ and 236℃(Wang and Kang, 2017).
3-d Inversion Of Gravity Data
Details of the gravity data
Satellite gravity observation technology has the characteristics of high coverage, high precision and high resolution, which greatly compensates for the lack of ground gravity measurement. Topographic data and satellite gravity data can be obtained from open source data access websites such as the International Center for Global Earth Models, the European Space Agency, and the Geoscience Research Center in Potsdam, Germany. Three professional gravity satellites CHAMP (Challenging Minisatellite Payload), GRACE (Gravity Recovery and Climate Experiment) and GOCE (Gravity Field and Steady-state Ocean Circulation Explorer) are all polar low-orbit satellites with high observation accuracy, and their observation range basically covers the world, and play an important role in the detection of the earth's gravity field. Based on the massive satellite gravity data and in-situ gravity data collected, major research institutions have established numerous earth gravity field models. Among them, the EIGEN-6C4 gravity field model was released in 2014, with a spatial resolution of 9km and the accuracy of gravity anomaly is up to 2.73mGal, it has a high accuracy in each ultrahigh order gravity field model. Therefore, in this article we use the gravity data provided by the EIGEN-6C4 gravity field model with a resolution of 1 arc minute, and the range of data obtained by satellite gravity is 99°36′E- 100°48′E longitude and 35°36′N- 36°36′N latitude (Ince et al., 2019). In order to facilitate the data processing, the Gauss-Kruger projection method is used to transform the satellite gravity data. Then the data after the coordinate conversion is gridded by Kriging interpolation, and the grid spacing is set to 1.5km. The gridded data is used as the research object and the interpolated satellite gravity Bouguer anomaly map in the research area is shown in Fig. 2.
3-D gravity focusing inversion method
In this section, we will introduce the 3-D gravity inversion method briefly and the method was proposed by Gao et al. (2019). The inversion method is realized by matlab language. The 3-D gravity inversion is to calculate the density or geometric parameters of the underground structure or anomaly based on the anomalous data on the observation surface. For linear inversion problems, the underground space needs to be discretized.
A typical method is to divide the earth region of interest into cells with constant density and fixed interfaces. With \(N\) observations and \(M\) rectangular prisms, the discrete forward modeling operator for the potential field problem can be written in matrix form
$$\begin{array}{c}d=Am (1)\end{array}$$
Where \(d\) represents ground observation data with length N, \(m\) is the vector of the model parameters with length M. In 3-D gravity inversion \(m\)represents the density of uderground medium. \(A\) is the sensitivity matrix, which is a rectangular matrix of size \(N\times M\). Usually, \(A\) is not a square matrix and the measured data is often much less than the parameters to be solved. In order to solve the serious multiplicity of solution and non-uniqueness problems in the solution process, the Tikhonov regularization method is usually used (Tikhonov and Arsenin, 1977). This problem is solved by minimizing the objective function. The objective function is composed of data fitting items, model constraints and regularization parameters, and its form is as follows:
$$\begin{array}{c}\varphi \left(m\right)={\varphi }_{d}\left(m\right)+{\alpha \varphi }_{m}\left(m\right) (2)\end{array}$$
Where \(\varphi \left(m\right)\) is the objective function to be solved, \({\varphi }_{d}\left(m\right)\) is the fitting items, \({\varphi }_{m}\) is the model constraints and \(\alpha\) is regularization parameters, The calculation formula of the objective function is as follows
$$\begin{array}{c}\varphi ={\varphi }_{d}+{\varphi }_{m}={‖Am-d‖}_{2}+\alpha {‖m‖}_{2} (3)\end{array}$$
In order to suppress the versatility and instability of the inversion as much as possible, it is necessary to add a model fitting term to the objective function. Model weighting matrix (\({W}_{m}\)) and data weighting matrix (\({W}_{d}\)) need to be introduced, The data weighting matrix has less influence on the inversion result, while the model weighting matrix has a greater influence on the inversion result. Since the sensitivity matrix decays with depth, gravity data faces a serious skin effect problem. As a result, the inversion results will be concentrated on the surface. The depth weighting function proposed by Li and Oldengburg(1996,199, )solves the serious skin effect problem of inversion, and its form is as follows
$$\begin{array}{c}w\left(z\right)=\frac{1}{{\left(z+{z}_{0}\right)}^{\beta /2}} (4)\end{array}$$
Where \(z\) is the depth of the center point of the underground space geological body, \({z}_{0}\) depends on the shape of the underground division unit and the observation altitude. In the inversion process, β takes a value of 2–3. Due to the low vertical resolution of the potential field, the inversion results often have serious tailing at the bottom. By using the improved depth weighting function proposed by Gao et al. (2019) makes the inversion result more concentrated,ma, es the vertical resolution of the inversion result higher. Its form is as follows
$$\begin{array}{c}w\left(z\right)=\frac{1}{{\left(z+{z}_{0}\right)}^{\beta /2}}\frac{1}{{\left(H-z-{z}_{0}\right)}^{\beta /2}} (5)\end{array}$$
Where \(H\) is the total vertical range of the inversion space.
Different norms can be used for model constraints. The model constraint term using the L2 norm constraint can produce a smooth solution, its suitable for the situations where the geological boundary is not obvious. However, the objective function using the focus constraint function can often better indicate the boundary of the ore body.In this paper we adopt the idea based on the minimum anomalous volume introduced by Last et al. (1983), and introduces the minimum support function into the objective function expression. Its form is as follows
$$\begin{array}{c}{W}_{e}\left(m\right)=\frac{{m}^{2}}{{m}^{2}+{e}^{2}} (6)\end{array}$$
In summary, the objective function form used in the inversion is as follows
$$\begin{array}{c}\varphi ={\varphi }_{d}+{\varphi }_{m}={‖{W}_{d}\left(d-Am\right)‖}_{2}+\alpha {‖{W}_{m}{W}_{e}m‖}_{2} (7)\end{array}$$
The model constraint term \({W}_{m}\) adopts the improved depth weighting function form, and uses the conjugate gradient algorithm to solve the above objective function. Various methods, such as quasi-Newton method, Gauss-Newton method(GN method) and conjugate gradient method(CG method), have been proposed to solve inverse problems. The CG method is widely-used in geophysics and is adopted in this paper because it has advantages of low dependence on the initial model, less memory, rapid convergence and can avoid matrix decomposition and matrix inversion. The calculation process is that the partial derivative form is set to \({g}_{k}\). The first iteration direction \({p}_{k}\) of the conjugate gradient algorithm is the negative direction of \({g}_{k}\), and the direction of each iteration thereafter is determined by the previous iteration direction and the new derivative function. In the first iteration process, the iteration direction is the negative direction of the partial derivative of the objective function, the search step size satisfies \({t}_{0}\frac{{p}_{k}^{T}{g}_{k}}{{p}_{k}^{T}{p}_{k}}\), and the model parameter \(m={m}_{k-1}+{t}_{k-1}{d}_{k-1}\). When the inversion objective function is less than the cut-off error or the number of iterations is less than the preset value of the number of iterations, the next iteration process is entered until the cut-off error value or the preset maximum number of iterations is met. And the factor k represents the number of iterations.
3-d Gravity Inversion Results
The gravity anomaly value is a comprehensive superposition result of anomalies generated by different depths and different density contrast sources underground. It is generally composed of two parts, the local field and the regional field. The regional field is caused by deeply buried and widely distributed rock masses, and the frequency of anomalies is low, and the amplitude of anomalies is large. The local field is caused by the rock mass with relatively shallow depth and small distribution, and its anomaly amplitude is small and the frequency is high. Anomaly separation is an important part of the work of gravity anomaly interpretation. In this paper, the most commonly used upward continuation method will be used to process the satellite gravity data in the study area, and then the local anomaly field in the study area will be obtained. The original gravity anomaly field is separated by the result of 30km upward continuation, and the local anomaly field obtained is shown in Fig. 3.
The residual Bouguer gravity anomaly shows that there is a wide range of low gravity anomalies in the middle of the study area, which may be related to the low-density rock masses existing underground. The Bouguer gravity anomaly in the central region is significantly lower than that in other regions. The 3-D physical property inversion for the residual gravity anomaly is performed in a Cartesian coordinate system, and the entire spatial range of the inversion is 111km×114km×40km. The underground grid is divided into 74×76×40 points, and the grid unit size is 1.5km×1.5km×1km. Selecting appropriate regularization coefficients, after 20 iterations of inversion, a more credible underground density distribution model is obtained.
Figure 4 shows the horizontal slice of the three-dimensional gravity inversion results in the study area at 2km, 5km, 8km 10km, 15km, 20km, 25km and 30km depth slices, respectively. There are widely distributed low-density anomaly areas below 2km underground of the study area, which has a good corresponding relationship with Quaternary and Neogene sediments. Obvious low-density banded anomalies appear on 5km and 8km sections. Combined with geological data, these banded low-density anomalies are interpreted as three proven concealed faults in the area, which are recorded as F1, F2 and F3, respectively. With the increase of horizontal slice depth, there is an obvious low-density anomaly area in the study area within 15km-25km. The low-density anomaly area is consistent with the distribution position and strike of the three known faults. According to the horizontal slice of the inversion results, the buried depth of the low-density abnormal body may be between 15km and 30km. Through the known hidden faults, there is a certain connection relationship between the low-density body and the surface.
We have selected four sections along the north-south direction (Fig. 5b-e), which can clearly see the widespread low-density anomalies that exist underground. The depth range is about 15km-30km, and the trend of hidden faults in the inversion results can be clearly seen above the low-density body. These faults can serve as channels for underground heat sources of hot dry rock geothermal resources in the Gonghe Basin.