Petrophysical parameters are a bridge between geology and geophysics, and are often used as the foundation for forward modeling and inversion (Eric, 1999; Wang et al., 2015; Xu et al., 2019). Density of about 180 samples was measured by a domestic digital DH-300 density meter with an accuracy of ±0.001 g/cm³. Information on the specific number, lithology, and sample locations is listed in Table 2. Samples include all lithology the shied-forming basalt, volcanic rocks, Late Triassic intrusion, Paleoproterozoic stratigraphic units, and Archean gneiss in study area, which help us decipher the gravity anomaly.

The samples from the Late Triassic intrusion exhibited the lowest density. The shied-forming basalt and volcanic rocks were characterized by a low to intermediate density. Archean gneiss showed an intermediate density, while the Paleoproterozoic stratigraphic units were characterized by highest densities. These differences in density allowed the need for geophysical research.

A total of 2,000 large-scale high-precision ground-based gravity survey points were collected using a CG-5 type Gravimeter (Scintrex Ltd., Canada) at a scale 1/50,000. The surveyed area was ~400 km² with a density of 4 to 6 measurement points per square kilometer. Subsequently, the data were gridded using the Gauss–Kriging method to provide detailed information of the geological bodies and structures (Han et al., 2020).

Corrections for Free-air and Bouguer anomalies were done using the International Association of Geodesy (1971) 1,967 formula considering 2.67 g/cm³ as the density of the middle layer (Morelli et al., 1972; Sundaralingam, 1987; Luo and Yao, 2007; Xu et al., 2022a). Corrections for the gravimetric terrain in near (0–20 m), intermediate (20–2,000 m) and remote areas (2–20 km) were performed using the Trupulse 200 type laser height instrument monitoring data, 1:10,000 scale elevation database and 1:50,000 scale elevation database, respectively. Finally, based on cone-shaped column formula, the topographic correction value of each gravity points are calculated (Laramie and Ralph, 2003; Xu et al., 2022a). The theoretical formula is listed as follows.

Where, Δg-Gravity anomaly (m/s²), ρ-Density of topographic correction (2.67 g/cm³), R-Topographic correction radius, i-Slope angle, n-Azimuth number.

Upwards continuation and first vertical derivative are two efficient methods for processing initial gravity data and each has its own merits (Lu et al., 2019; Xu et al., 2022b). The advantage of upwards continuation is that it can segment the deep Bouguer anomalies from the background. Various heights in the range of 200–1,000 m were adopted to highlight the various depth anomalies. In contrast, the advantage of the first vertical derivative is that, it enhances the local anomalies caused by subsurface geological bodies (Cooper and Cowan, 2004). The present study maximized the advantages of both methods for processing the original gravity data.

Further, two point five dimensional man-machine interactive inversion technology is a perfect compromise between quantitative interpretation of the complicated anomalies and varying the theoretical parameters (Pinto and Casas, 1996). It is used to decipher the four elements of deep geothermal resource (heat source, geothermal reservoir, caprock, and geothermal migration pathway) (Ahumada et al., 2022). The first step for the prediction model was established according to the estimated geological bodies. The theoretical gravity curve was subsequently generated based upon the estimated geometry of the density distribution. The shape of the model is terminated at the least difference between the calculated and measured anomalies (Yao et al., 1998). Finally, hypothetical geological bodies were accomplished by two dimensional interactive inversion. The theoretical formulas and computational model are listed as follows (Figure 3) (Luo and Yao, 2007; Fan et al., 2014).

Where, Δg-gravity anomaly (m/s²), P-position, G-constant value of earth's gravity (N•m²/kg²), σ-prism density, i-prism corner number, N-number of prism, I-incidence, y-Coordinate in y direction, w_i-Center angle (rad/s).

The interpretation of gravity data involves i) the division of sub-regions with striking discrepancies in the gravity field, ii) geological inference of Bouguer anomalies, and iii) the determination of the spatial distribution of anomalies (De Castro et al., 2014). Figure 4 presents the distribution of the Bouguer anomalies in the north of Tianchi caldera. The range of remarkable lowest Bouguer anomaly A zone is from −57 × 10⁻⁵m/s²−66 × 10⁻⁵m/s², the B anomaly zone is from −47 × 10⁻⁵m/s²−56 × 10⁻⁵m/s², and the C zone is from −34 × 10⁻⁵m/s²−46 × 10⁻⁵m/s².

These results indicate the existence of density difference in north of Tianchi caldera. According to the descriptions of geological setting and petrophysical parameters in Section Petrophysical parameters, the low anomaly is related to volcaniclastic rocks, and the high anomaly is related to the Paleoproterozoic stratigraphic units or Archean gneiss (Figure 4).

However, the boundaries between volcaniclastic rock and Paleoproterozoic stratigraphic units or Archean gneiss are faintly discernible. The results of upwards continuation with different heights (200 m to 1,000 m) indicate that the value of low and high gravity anomaly decreased rapidly with increasing depth, leaving two dome shaped low anomalies (Figures 5a–c). The first vertical derivative results highlight the scope of gravity gradient (Figures 5a'–c').

There is a close link between gravity value and geological factors, such as petrophysical properties and geological structures (Yao et al., 1998; Xu et al., 2022b). Residual gravity anomaly has been separated from the Bouguer gravity for carrying out two point five dimensional man-machine interactive inversion technology (Figure 6). Then, two orthogonal gravity sections (P8 and P9) were extracted from the residual gravity anomaly.

P8 and P9 showed a similar trend with a C-type gravity distribution. The gravity value varied from −6 × 10⁻⁵ m/s²−4 × 10⁻⁵ m/s² (Figures 7, 8). This trend demonstrates that the density in middle are lower than that in both ends and there exits the geological structures crowd zones.

On the basis of surface exposure and statistical distribution in Section Petrophysical parameters, the theoretical gravity curves (red lines) are matching to the field surveyed gravity curves (black lines) and are deducing deep geological bodies and structures (Figure 7). From west to east, the high density geological body (~2.77 g/cm³) is the Archean metamorphic crystalline basement, the low ones are Late Permian monzogranite (~2.55 g/cm³) and Middle to Late Jurassic pyroclastic rock (~2.60 g/cm³), and the high one is Paleoproterozoic marble (~2.79 g/cm³). Because of the large porosity, loose cementation and high permeability, the lower pyroclastic rocks are generally regarded as the geothermal reservoir. The gravity gradient zones (sties 2,500 and 4,000) indicate that there exit geological structures which provide migration pathway for deep geothermal resource. In addition, the shied-forming basalt shield (~2.59 g/cm³) in the shallow overlying the geothermal reservoir acts as a caprock reducing heat loss function (Figure 8).

As the only active volcano in Changbaishan volcanic field, Tianchi volcano has shown huge potentials in geothermal resources (Zhang et al., 2018). More than 120 hot springs congregate around the Tianchi caldera (Xu et al., 2021). However, due to different control conditions of geothermal resources, uplift mountain or sedimentary basin geothermal systems show different geophysical anomaly (Zhang et al., 2018). Taking the Magmatic type (I₁ type in uplift mountain geothermal system) as example, continuous advection of heat from Tianchi volcano which provides powerful thermal energy and thermal cycle to directly heat the hydrothermal reservoirs is the hot source (Table 3). The geothermal gradient or radioactive element decay with the geothermal gradient >3.0°C/100 m provides heat for medium to low geothermal resource in the sedimentary basin geothermal system.

Effective assessment of potential of different types of geothermal resources is crucial and difficult. Results from the gravity planar overall indicate that the low gravity anomaly should be corresponding to the sedimentary basin and the high anomaly is related to the uplift mountain. Then, the advanced data processing technologies (upwards continuation and first vertical derivative) attempts to precisely assess geothermal resource potential in sedimentary basin and uplift mountain (Figure 9). In addition, based on the two point five dimensional man-machine interactive inversion technology, gravity sectional results indicate that the cataclastic marble and deep fault are geothermal reservoir and migration pathway of uplift mountain geothermal system, respectively. To the sedimentary basin geothermal system, pyroclastic rock is the primary geothermal reservoir and can migrate geothermal energy without fault. It is the formation thickness that determines the success of sedimentary basin geothermal system. Combing with regional geothermal well, temperature, lithology and water yield data, the geothermal resource in North Tianchi caldera can be catagorized into the sedimentary basin and uplift mountain geothermal system (unpublished data) (Figure 10).