Microtremor Survey Method: A New Approach for Geothermal Exploration

Geothermal resources are a type of sustainable and green energy, which can play an important role in emission peaks and carbon neutrality. Determining the best development target area is key to resource development and geophysical methods are commonly used for this purpose. Owing to serious human and industrial interference, the microtremor survey method is often adopted for geothermal exploration in urban areas. It is a passive source method, which is non-invasive and environmentally friendly. In this method, the Rayleigh wave dispersion curve is extracted using spatial autocorrelation based on the vertical component signal at the observation station. A genetic algorithm is used to invert the dispersion curve of one survey point to obtain strata parameters such as layer thickness, S-wave velocity, and density. It provides critical parameters for the cap layer and reservoirs for geothermal exploration. For a chain microtremor measurement, a two-dimensional (2D) apparent S-wave velocity section can be generated. The apparent S-wave velocity is calculated from the phase velocity using the following empirical method: the 2D apparent S-wave velocity section helps to identify the buried channel for heat flow and track the irregular shapes of the reservoirs or cap layers. It has been verified that the microtremor survey method is reliable and accurate compared with borehole materials. As a newly developed non-invasive geophysical method, it can be widely used in geothermal exploration.


INTRODUCTION
High-frequency climate and environmental disasters push mankind to pay more attention to protecting the environment. One of the best methods is energy conservation and emission reduction. Switching from traditional fossil fuels to renewable and green energy decreases the emission of carbon dioxide (CO 2 ). Renewable and green energy mainly include solar, wind, biomass, and geothermal energy.
Geothermal energy comes from the interior of the Earth, and is stable because it does not depend on the weather. It originates from the decay of radioactive elements contained in the Earth. As a result, the amount of geothermal energy is large, equal to 170 million times that of coal with a potential exploitable resource of 4,948 trillion tons of standard coal equivalent. China is located in the circum-Pacific geothermal belt and the Mediterranean-Himalayan geothermal belt, where geothermal energy is abundant. Geothermal types include shallow geothermal energy, mediumdeep depth hydrothermal resources, and hot dry rock geothermal resources. It is used for various proposed, such as generating power, heating, and cooling based on the different types of geothermal resources. Identifying an excellent target for exploiting geothermal resources is the most important factor, which typically relies on geophysical methods to ensure geothermal reservoirs and channels for heat flow (Kana et al., 2015).
Geophysical methods can be divided into invasive and noninvasive methods. Seismic methods, magnetic methods, and gravity methods are non-invasive methods. Based on the wave velocity difference of the strata, seismic methods can identify strata with different rock properties and determine the location and parameters of buried faults. Magnetic methods such as the magnetotelluric method (MT), audiomagnetotelluric method (AMT), and controlled source audio-frequency magnetotelluric method (CSAMT), can explore the distribution characteristics, properties, and geological fault occurrences and determine the thickness and depth of cap layers. Gravity methods help to depict the bedrock conditions, spreading of structures, and location of faults. For example, a 3D seismic survey was conducted in Saxony, Germany for geothermal exploration (Lüschen et al., 2015), which is an indispensable tool for geothermal exploration, even in crystalline basement rocks. To identify the geothermal potential, seismic reflection imaging can provide amplitude contrasts that may be related to volcanic rocks and possible hydrothermal alteration fronts (Sena-Lozoya et al., 2020). MT is an important method for the exploration of geothermal systems (Rosenkjaer et al., 2015). Uchida and Sasaki (2006) developed a stable inversion technique for the 3D interpretation of MT data and applied it to a large-volume MT dataset obtained in the Ogiri geothermal area, in southwestern Japan. CSAMT provides an efficient means of delineating the shallow resistivity pattern above a hydrothermal system (Sandberg and Hohmann, 1982). Di et al. (2006) successfully applied the CSAMT method to geothermal exploration in the Beijing area and accurately located geothermal reservoirs deeper than 2 km. Gravity methods can help to delineate the structural features that control the geothermal system, such as the research on the Eburru geothermal system in Kenya (Maithya et al., 2020). To understand the subsurface structure and its relation to the observed geothermal phenomenon, a land gravity survey was carried out in the Kinigi geothermal field, Northwest Rwanda (Uwiduhaye et al., 2018). In addition, joint inversion is typically adopted for geothermal exploration. Zaher et al. (2018) used airborne gravity and magnetic geophysical data to preliminarily explore the geothermal potential in the Siwa Oasis, Western Desert, Egypt. Joint inversion of gravity and surface wave data constrained by MT highlighted unconventional geothermal prospects in felsic basement (Ars et al., 2019).
Ambient noise tomography has also been adopted for geothermal exploration (Planès et al., 2020). The microtremor survey method (MSM) is a new, non-invasive geophysical method used for geothermal exploration, which also uses ambient noise as a signal. An observation array is typically adopted to collect field data. The frequency-wavenumber spectrum method F-K, (Capon, 1969), and spatial autocorrelation method SPAC, (Aki, 1957), are commonly adopted to extract the Rayleigh wave dispersion curve. Compared with F-K, SPAC has more accurate results (Ohori et al., 2002) with circular arrays. Bettig et al. (2001) presented a modified SPAC, which allows the use of irregular, almost carbitrarily shaped arrays. Okada and Ling (1994), Okada (2006) proposed an extended SPAC method with arbitrarily shaped arrays. Asten (2006) used the SPAC method to process the passive seismic data from finite circular array. Owing to the fundamental assumption of plane wave propagation of surface waves, attention must be paid to the near-source effect in array-based microtremor surveys (Roberts and Asten, 2008). Asten et al. (2019) developed a spatially averaged coherency (krSPAC), which allows spatial averaging of spectra from multiple pairs of sensors, irrespective of differences in spatial separation of the pairs. Cho and Iwata (2021) examined the limits and benefits of the SPAC microtremor array method. Microtremors have been applied to various aspects based on different processing methods. Roberts and Asten (2005) used microtremor array (SPAC) measurements to estimate the S-wave velocity profile of Quaternary silts. Chávez-García et al. (2007) adopted microtremors to achieve the S-wave velocity structure around the Teide Volcano. Stephenson and Odum (2010) applied the SPAC to characterise the S-wave velocity in the upper 300 m of Salt Lake Valley, Utah. It is also possible in SPAC observation to identify a 2D structure and interpret S-wave velocity profiles in a deep and narrow valley environment (Claprood et al., 2011). MSM has been used to determine soil characteristics (Ozer et al., 2017), estimate the near-surface S-wave velocity structure at a site exhibiting low to high impedance contrast (Setiawan et al., 2019), and identify the sediment (Raptakis and Makra, 2010;Tian et al., 2019;Tian et al., 2020a). Tian et al. (2016; applied the MSM method to identify the parameters of different strata and ensure the geothermal channel information. In addition, Tian et al. (2020b) combined MSM with geothermal assessment to improve precision. Using the regression method, the relationship between S-wave velocity, density, and thermal conductivity was identified (Tian et al., 2020c).
In this study, we investigated MSM, which can be widely applied to geothermal development and utilisation. First, we briefly introduce the observation array and the SPAC method. For different exploration targets, MSM is classified into two specific techniques: microtremor sounding and 2D microtremor profile methods. The inversion results of the dispersion curve provide layer parameters (depth, thickness, S-wave velocity, and density) for the geothermal system. To identify the geothermal channel, a 2D microtremor profiling method was adopted, which contained several survey points in a survey line. Our preliminary results appear to be useful for studying geothermal systems. MSM is a passive surface-wave geophysical method. Microtremor signals mainly contain surface waves and body waves. The low-frequency waves (<1 Hz) arise from changes in pressure, sea waves, and tides, while the high-frequency (>1 Hz) waves originate from human sources and living activities. Based on the SPAC requirements, the microtremor signal must be collected by an array. Commonly, a regular triangular array ( Figure 1) is selected to record the microtremor signal. The Rayleigh wave dispersion curve is extracted from the vertical component signal using the SPAC method (Aki, 1957), which assumes that the signal is stationary and random. Meanwhile, the surface wave is mainly the fundamental mode (Okada, 2003;Okada, 2006).
Usually, the frequency spectrum of the microtremor signal for two observation points A (0,0) and B (r,θ) can be expressed as; where ω is the angular frequency, θ is the azimuth angle, k is the wave number, and r is the distance (radius) between the two stations.
The SPAC function between two stations is defined as Here, the spatial covariance function is defined as For the circular array, the azimuth average of the spatial covariance function for the SPAC function and the phase velocity can be obtained.
Here, J 0 is the zero-order Bessel function and h is the average power.
Then, the SPAC coefficient ρ(ω, r) is defined in the angular frequency ω.
As a result, we can obtain Using Eq. 7, the phase velocity can be calculated in distance r and frequency f.
When the Rayleigh wave dispersion curve is extracted, it is used to invert the S-wave velocity structure using the Monte Carlo algorithm. In our research, a genetic algorithm was selected to invert the geological structure. During the inversion, a initial inversion model with layer's parameters is assumed to begin the inversion which is one horizontally stratified geological model. If the observation array is located at the fractured strata, the results of inversion has a relatively lower accuracy.
In some cases, it is not necessary to invert the S-wave velocity structure during the detection of geological structures. The relative change trend such as low velocity anomaly in the S-wave velocity section can reveal the location and depth of the structure accurately. The 2D microtremor profiling method is also based on the Rayleigh wave dispersion curve. It transforms the phase velocity (Vr-f) to the apparent S-wave velocity (Vx-H) using an empirical Eq. 8 for each survey point (Tian et al., 2017;Tian et al., 2022). Then, the apparent S-wave velocity section is obtained by horizontal interpolation and smoothing between every survey point. The workflow of MSM is shown in Figure 2.
Here, the Vx represents the apparent S-wave velocity of the ith layer, t represent the time, Vr represents the phase velocity of the ith layer.

DATA ACQUISITION, PROCESSING, AND DISCUSSION
According to the exploration target, different techniques were selected to determine geothermal geological parameters. This section introduces three objectives for geothermal exploration. Different observation parameters and inversion methods were presented specifically to satisfy each item.

Determining the Geothermal System Parameters
Geothermal cap layers and reservoirs are known for some geothermal fields. Therefore, it is critical to determine the corresponding parameters. Accurate parameters of all strata are important for the different stages of geothermal development. For example, during the feasibility step, the estimation of the geothermal deposit and exploitation amount depends on the reservoir thickness. In addition, the reservoir Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 817411 depth is a vital parameter for drilling, which may decide whether the geothermal development is successful or not.

Data Acquisition
JingJiang, China is taken as an example where the geothermal reservoir is Triassic limestone. Based on the probing depth requirement, we deployed one survey point with four overlay circular observation arrays (Figure 1) to record the microtremor signals. The radii were 80-160-320-640 m. At every observation point, a three-component MTKV-3C seismometer and Datamark LS-8800 were used to record the microtremor signals. The main frequency of the observation seismometer was 1 Hz. At the same time, one cycle delay and amplifier circuit was adopted to achieve lower frequency signal. At each survey point, a 30 minute continuous signal was obtained.

Data Processing and Discussion
First, the original microtremor signal was processed by removing the data mean and trend, normalization and spectral whitening. Based on the original field data, the SPAC coefficients were calculated for every two stations. Then, they were fitted to the zero-order Bessel function for different frequencies and distances. After the dispersion curve was extracted, a genetic algorithm was adopted for inversion. The dispersion curve and inversion S-wave velocity structures are shown in Figure 3. Simultaneously, one borehole was drilled at this microtremor survey point. A comparison between the inversion results and borehole materials is shown in Table 1. The MSM distinguished several subsurface layers based on the wave velocity differences. In the same geologic age, it also distinguished different layers if there was a wave velocity difference. The error was small when the wave velocity difference was large. For example, it accurately identified the interference between the Quaternary Neogene and the Zhouchongcun Formation of the Triassic. However, the probing error was slightly larger when identifying the interference between two layers with small wave velocity differences, such as the layer of cream mudstone, gypsum, and layer of gypsum, creating mudstone. In general, the MSM can provide accurate parameters for layers with different rocks.

Identifying the Buried Channel for Heat Flow
Convective geothermal resources are one resource type, which is frequently chosen for development and exploitation, which relies on channel formation for heat flow. Therefore, locating the buried channel and determining the channel parameters are the primary tasks for geothermal exploration. Faults (fracture zones) are channels for heat flow. To identify buried faults, one 2D section was the most intuitive and accurate method. In the horizontal direction of the section, the lateral position was determined. This section provided an accurate understanding of the width and depth of the channel. It also showed the fracturing degree of the fault by comparing the background values.

Data Acquisition
A research area in Zhejiang Province, China was selected to depict the ability of the MSM to identify the buried channel for heat flow. To obtain a 2D section, one microtremor survey line was deployed, including six survey points. The observation system is shown in Figure 4. For each observation array, the radii are 125-250-375-500 m using a 1 Hz main frequency for the observation seismometer. The time was automatically corrected using the GPS system among all seismometers. The microtremor signal was recorded for 30 min in each observation array.

Data Processing and Discussion
The field data were processed by removing the data mean and trend, normalization and spectral whitening. For every survey point, the dispersion curve was extracted using the SPAC method. All the dispersion curves are shown in Figure 5A. To achieve an apparent S-wave, the dispersion curve was transformed to the curve of the apparent S-wave velocity (Vx) and depth (H) using Eq. 8. Using the results of Vx-H, the 2D apparent S-wave velocity Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 817411 section ( Figure 5B) was calculated using Kriging interpolation and smoothing among all the survey points. At the survey point located in the fault, the dispersion curve showed a low phase velocity at a low frequency compared with other dispersion curves. The apparent S-wave velocity section clearly indicated a buried fault. The strata can be divided into several layers. The depth of the first layer was 0-100 m, with a Vx of 300-400 m/s. The Vx of the second layer varied from 600 to 900 m/s at depths of 100-350 m. The lithology barely changed in this layer, and fracture zones developed in the third layer. An obvious low-velocity anomaly was observed from 800 to 1,500 m at distance of 600-650 m. The Vx values of the low-velocity anomaly and normal layer were approximately 1,600 m/s and 3,000 m/s, respectively. This anomaly is named Fault F1, and its width is approximately 50 m. Compared with the fracture zone, the surrounding rock is intact with high wave velocity. Geothermal resources are strongly related to this fault, which breaks and cuts off the fresh bedrock, resulting in fracture zones. This provides an excellent channel for heat flow.
The 2D microtremor profiling method has relatively high precision for probing geological structures. In the section of apparent S-wave velocity, the strata that are stably deposited (uniform colour) could be continuously traced laterally. Areas with obvious lithology changes, such as fault fracture zones, showed low-velocity anomalies because of the underestimation of phase velocity using SPAC. The obvious low-velocity anomaly  FIGURE 4 | Observation system for a microtremor survey line. In this figure, R represents the observation radii and D represents the distance between two survey points.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 817411 can be used to interpret and identify buried faults (fracture zones) and other geological structures in the apparent S-wave velocity profile.
According to the comparison study and analysis of the detection results of the 2D microtremor profiling method and drilling results, it can be seen that the buried fracture zone displayed an obvious low-velocity anomaly in the 2D apparent S-wave velocity profile, which is an important indicator explaining the buried geothermal structure. Using the 2D micromotion profiling method to detect the thermal control structure and exploit the convection-type geothermal resources, the detection accuracy was high, and the measured effect was good. There were obvious lithologic differences on both sides of the normal fault plane, and the detection effect was good.

Tracking Irregular Strata Shapes
Estimating the amount of geothermal energy and exploiting geothermal resources requires a comprehensive understanding of the upward and downward trends of strata. The diagenetic environment of different geological periods and the function of tectonic movement in the test area resulted in complex geological structures. As a result, the strata were not homogeneous. Geophysical sections imaged the subsurface to generate a prospective stratigraphic structure. For conductive geothermal resources, information on reservoirs helps to calculate the total reserves of geothermal resources and locate the best drilling position.

Data Acquisition
Here, one research area in Shandong Province, China, is taken as an example to track the shape of different layers using a 2D apparent S-wave velocity section. This research area mainly exploits shallow geothermal resources above 150 m. The microtremor survey line included four survey points. The observation radii were 5-10-20-40 m. The main frequency of the observation seismometer was 2 Hz. The microtremor signal was recorded for 50 min in each observation array.

Data Processing and Discussion
The processing workflow for extracting the Rayleigh wave dispersion curve was the same as that for the above two sections. All the dispersion curves of the survey points are shown in Figure 6A. In addition, the dispersion curve of every survey point was transformed into Vx-H. The section was plotted Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 817411 6 using lateral Kriging interpolation and smoothing of the survey points. Figure 6A shows that the phase velocity of survey points RW1-4 is larger than that of the others. The dispersion curves of RW1-1 to RW1-3 exhibit the same spreading trend. The apparent S-wave velocity section shows that the strata are complex and variable ( Figure 6B). Soil and mudstone develop from 0 to 10 m near survey point RW1-3. It only develops around survey point RW1-3 and stretches approximately 100 m horizontally. The Vx of the soil and mudstone was approximately 350-390 m/s. In the depth range of 25-40 m, a low-velocity layer was developed, which disappeared at a horizontal distance of 270 m with a Vx of 300-350 m/s. In this research area, weathered granite was developed to different degrees. The strongly weathered granite had a Vx of approximately 600 m/s leading to wide dispersion. Below this layer, medium-weathered granite also developed widely. Its thickness is approximately 40-50 m and the Vx is approximately 900-1,000 m/s. The Vx of the fresh granite was approximately 1,100-1,200 m/s. Below a depth of 125 m, the bedrock is intact, and the Vx was greater than 1,400 m/s. At survey point RW1-3, drilling to a depth of 102 m was emplowed to verify the effectiveness of the MSM. The rock core of this borehole included soil, mudstone, strongly weathered granite, medium-weathered granite, and granite. Based on the correlation with these drilling materials, the 2D microtremor profile method could distinguish weathered granite with different degrees. In addition, it accurately determined the bedrock. Using the 2D apparent S-wave velocity section, we could track the developing trends of different rock layers, including geothermal reservoirs and cap layers. Meanwhile, it provided geological information on the bedrock, which plays an important role in the development of shallow geothermal resources.

CONCLUSION
Geothermal resources have received much attention because of the requirement for energy conservation and emissions reduction. As one kind of green energy that comes from the Earth's interior, it is stable and sustainable. Geothermal amounts above 3 km of, subsurface and generating power, can be widely used for geothermal resources, such as warming, cooling, and generating power. Geophysical methods are usually adopted to provide essential information to search for an excellent developing target area.
The MSM is a new approach that has been successfully applied to geothermal exploration. It adopts ambient noise as a signal to avoid the inconvenience of an artificial signal source. In addition, the observation array for recording the field data is independent. Compared with traditional geophysical methods, the MSM is flexible and economical. In particular, it can be deployed under the interference of human and industrial activities. Because most of the development and utilisation of geothermal resources is located in cities or towns, this method is regarded as one of the most effective.
According to different exploration purposes, MSMs can be classified into two types: microtremor sounding methods and 2D microtremor profiling methods. The microtremor sounding method uses a genetic algorithm for inversion. It provides the parameters of the different rock layers and the S-wave velocity, layer thickness, and density are obtained from the inversion. In some research areas, cap layers and reservoirs are known. For geothermal exploration, it is necessary to provide the depth, thickness S-wave velocity, and density. The microtremor sounding method can accurately identify the related parameters of the different layers. The 2D microtremor profiling method must deploy several survey points to achieve an apparent S-wave velocity section. It can be used to identify the location and depict the developed characteristics of the buried channel. In addition, it can track the irregular shapes of cap layers or reservoirs. Based on the comparison between the results of the MSM and borehole, it is verified that this method is accurate and reliable for geothermal exploration. The development of clean and green energy is a new approach for geothermal exploration.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.