Monitoring seismicity triggered by geothermal site shutdown with a surface DAS array at Brady Hot Springs

SUMMARY

Seismicity induced by fluid injection including wastewater injection, hydrofracking and enhanced geothermal system (EGS) site production draws public attention. Dense arrays have been deployed to improve monitoring capability. In 2016 March, the PoroTomo experiment deployed an 8.6-km-long fibre-optic cable for distributed acoustic sensing (DAS) in the geothermal field at Brady Hot Springs, Nevada, covering an area of 1.5 by 0.5 km. The goal was to assess an integrated technology for characterizing and monitoring temporal changes in the rock mechanical properties of an EGS reservoir in three dimensions. We applied a neural network designed for earthquake detection called ADE-Net2 to the DAS data set to detect seismic events in continuous records. We were able to detect a total of 90 earthquakes, which included 21 events that had not been reported by a previous template-matching study. Additionally, we were able to successfully detect almost all of the active source signals, with only seven events being missed. We used the STA/LTA (short-/long-term average) method to pick arrivals and a clustering method to remove outliers. We initially tried a standard event location algorithm, but the low signal-to-noise ratio resulted in significant picking uncertainty that is up to ∼0.5 s, leading to large location uncertainty. Therefore, we developed a new location method based on the similarity between the theoretical traveltime curve and picked moveout. A grid search scheme was adopted to find the optimal point at which the traveltime curve is most similar to the picked one. Most newly detected earthquakes locate southwest of the DAS array, where five earthquakes were reported by a local seismic network. The plant began shutting down at 19:15 UTC on the March 13, and most earthquakes occurred on the March 14, indicating a relationship between the seismicity and the pressure changes caused by the shutdown of the plant. The pressure changes at epicentres obtained from a simplified model range from 71 to 157 kPa, exceeding a typical earthquake trigger threshold of 10 kPa.

INTRODUCTION

Human-induced earthquakes have been studied for a long time. Following a large number of earthquakes caused by wastewater injection at Rocky Mountain Arsenal in the 1960s (e.g. Bardwell 1966; Healy et al. 1968), induced earthquakes have been reported in several industries, including mining (Li et al. 2007; Qian et al. 2019), reservoir (Talwani 1997; Cheng et al. 2012), wastewater injectioechnology satisfies the requirements of both high-density observation and real-n (Horton 2012) and oil and gas development (van Eck et al. 2006; Sheng et al. 2020). Recently, shale gas and enhanced geothermal systems (EGS) have gained more investment to provide more energy options. However, both of them rely on fluid injection, which is supposed to be related to increased seismicity at the operation sites (e.g. Ellsworth 2013; Schultz et al. 2020; Wang & Manga, 2021) and raised public concerns (e.g. Keranen & Weingarten 2018). Studies on induced earthquakes are relevant not only to industrial operations but also provide new opportunities for earthquake science (Ellsworth et al. 2016; Galis et al. 2017; Schill et al. 2017).

Induced seismicity is frequently reported during EGS operations as low-temperature fluid is continuously injected into the reservoir and high-temperature fluid is extracted from production wells (Olasolo et al. 2016; Jung et al., 2013; Yin et al. 2021). One recent example is the 2017 M_w 5.5 Pohang earthquake, which occurred during geothermal power generation (Grigoli et al. 2018; Kim et al. 2018). Studies of the 2017 Pohang earthquake show that it occurred on an unmapped fault that was activated by injected high-pressure fluid (Ellsworth et al. 2019 William et al. 2019). Such a strong earthquake leads to strong negative community impacts and may result in the termination of an EGS project. For example, the EGS project in Basel, Switzerland was terminated after a series of moderate earthquakes (Diehl et al. 2017). Therefore, the industry has proposed the establishment of a seismic monitoring network to guide pumping based on seismic activity, such as the Traffic Light System (e.g. Kwiatek et al. 2019; Schultz et al. 2021), and this idea has been initially validated in several geothermal fields (e.g. Ader et al. 2020).

To improve real-time monitoring capability, dense surface arrays (e.g. Zeng et al. 2014; Roux et al. 2014) or downhole arrays (e.g. Maxwell et al. 2010) are often employed. A downhole array offers better monitoring capabilities but with high cost, and limited azimuth coverage affects location precision (e.g. Eisner et al. 2010). A surface array provides better azimuth coverage, but it suffers from stronger noise as well as weaker signal due to attenuation. Array techniques have been proposed to improve signal quality via stacking traces, which helps to improve monitoring capability. A larger aperture array with a smaller station spacing generally yields better detections and locations. (Li et al. 2018) However, it requires deploying more geophones and high-speed data transmission for big data volumes, which significantly increases the cost. Distributed acoustic sensing (DAS) technology satisfies the requirements of both high-density observation and real-time data transmission to a certain extent (Song et al. 2021). It has been applied to temperature change monitoring in permafrost (Ajo-Franklin et al. 2017), physical quantities measuring (Becker et al. 2017; Becker & Coleman 2019), measuring changes in liquid levels, depths and fluid infiltration pressures in soil and hydrological environments. It has also been used for EGS induced seismicity monitoring (Lellouch et al. 2021), exploration and monitoring (Henninges et al. 2019). There are also challenges in the application of DAS. For instance, the ground coupling effect of DAS is often weaker than that of traditional seismometers, and the noise level of DAS is usually higher (McClellan et al., 2018). Additionally, the factors that affect waveforms and amplitudes are not as clear (Hornman 2017). In the PoroTomo experiment, the deployed fibre optic cable was close to a highway, which introduced more noise from traffic signals. Moreover, as large-N arrays, DAS generates data that is orders of magnitude larger than that of traditional seismometers, making automatic detection algorithms essential. Seismic detection algorithms developed by using artificial intelligence technology can achieve efficient real-time processing on large volume data sets (e.g. Strok et al. 2020). Applications of AI combined with DAS for signal detection are widespread, ranging from microseismic detection (Huot & Biondi 2020; Stork et al. 2020; Huot et al. 2022), to traffic detection (Martin et al. 2018) to footstep detection (Jakkampudi et al. 2020). This study uses DAS surface array data collected by the PoroTomo team at Brady Hot Springs, Nevada (Feigl et al. 2017) to study the earthquake sequence during the experiment (Cardiff et al. 2018). Li & Zhan (2018) applied a template matching algorithm to DAS data for event detection, which is based on multi-channel cross-correlation. By utilizing templates from five catalogued events included in Table 1, they successfully detected a total of 116 events, 68 of which were significantly below the noise level, demonstrating the potential of using this approach to investigate seismicity response in detail. In this study, a newly developed earthquake detection method called ADE-Net (Lv et al. 2022) was used to scan the 14 d of continuous data. Dozens of earthquakes were newly detected and located with a novel location method and the location results are used to assess their relationship to the pore pressure changes.

BACKGROUND ON THE GEOTHERMAL OPERATION AT BRADY HOT SPRINGS, NEVADA

In 2016 March, the PotoTomo team conducted a series of surveys in the geothermal field at Brady Hot Springs, Nevada to investigate the rock mechanical properties and temporal changes during changing plant operations (Feigl et al. 2017). The injection wells are located in the northern part of the field and the faults with high permeability are presumed to act as fluid channels to the deep reservoir (Cardiff et al. 2018). High-temperature (420 K) brine is extracted from six production wells in the southern part of the field (Fig. 1). The 14-d PoroTomo experiment period consisted of four phases (Feigl et al. 2017), including (1) all injection and production wells working normally; (2) power plant was shutdown, no pumping occurred during this phase; (3) all production and injection wells were reactivated with varying injection rate, except for the far-field injection well and (4) return to normal operation as Phase 1. The suspension of long-term extraction can potentially lead to an increase in ambient pore fluid pressure near the production wells, as observed in previous studies that have shown a clear correlation with the seismicity at Brady Hot Springs (Davatzes et al. 2013). The various operating conditions during the experiment provide an opportunity to capture seismicity response to pore pressure changes.

The local seismic network (Foxall 2016) reported five local earthquakes (Table 1) during the shutdown and a previous study (Cardiff et al. 2018) explored their relation to pore pressure change due to plant shutdown. As part of the PoroTomo seismic array, an 8.6 km long fibre-optic cable was deployed in shallow trenches and a DAS interrogator was used to continuously record the seismic wavefield with a channel spacing of 1 m. Using data from the five catalogue earthquakes as templates, Li & Zhan (2018) detected more than 100 earthquakes, demonstrating the feasibility of using these data to investigate seismicity response in detail. However, the signal quality of most of the detected events in their study is too poor to allow events to be located and the limited number of templates may also lead to missed detections.

DETECTION

Neural network

To perform seismic event detection from continuous data, we designed a neural network named ADE-Net2 that is modified from ADE-Net (Lv et al. 2022). ADE-Net2 uses three convolutional blocks for feature extraction and two fully connected layers for classification (Fig. 2). In the Lv et al. (2022) study, ADE-Net demonstrated its ability to detect events with low signal-to-noise ratio (SNR) signals. However, the size of the ADE-Net is so large that it leads to slow training convergence. In addition, neural networks are more likely to overfit with a small-size training data set. Since the features, which are shown as continuous gradient change, extracted by the ADE-Net are relatively simple, it is possible to simplify the network to speed up training and avoid overfitting. The number of filters in each convolution layer is reduced by a factor of six. Inspired by the Fully Convolutional Network (FCN), a global pooling layer is added at the end of the convolutional block, which downsamples the input representation by taking the maximum value (Fig. 2) to replace the flattening layer. There are only 16 870 trainable parameters in the refined network (ADE-Net2), which is about 2 per cent of the original ADE-Net. We tested with 100 windows, and the original network took 752 ms, whereas the streamlined network averaged 292 ms, resulting in a reduction of over 60 per cent in processing time.

As in ADE-Net, each convolutional block is composed of one convolutional layer, one pooling layer and one activation layer. We employed the backpropagation method for updating the weights and biases of the neural network. This method involves calculating the gradient of the loss function with respect to the model parameters and adjusting the values accordingly during the training process. To prevent gradient instability problems that can arise from the chain rule, ReLu (Rectified Linear Unit) is commonly used as the activation function in neural networks. This function increases the network's nonlinearity by setting all negative input values to zero while maintaining positive values, thus preventing gradient disappearances and explosions (Agarap 2018). Additionally, ReLu is preferred for its simplicity and computational efficiency. The input layer of our neural network consists of a 512 × 512 matrix, where the first dimension represents different channels in the spatial domain, and the second dimension represents the time-series, which is sampled at 100 Hz. The outputs are probabilities of earthquake and noise, which is the same as the original ADE-Net. If the probability of an earthquake given by the neural network is greater than that of noise, a potential event is reported for that time window.

Pretreatment

The quality of the input directly affects the performance of the network. Since the coupling and background noise level across the whole array varies strongly, we pre-processed the raw data to improve the SNR and for equalization (Lv et al. 2022). Following Lv et al. (2022), bandpass filtering and the short-/long-term average (STA/LTA) method were applied to the raw data. The bandpass filtering frequency band is chosen to be 20–25 Hz, which is effective in suppressing the local noise. We find that with a short window of 0.2 s and a long window of 1 s, the earthquake signal can be enhanced via the STA/LTA method.

We used only four earthquakes on March 14 in the catalogue (Table 1) to build the training data set (Fig. 3), which may lead to overfitting due to the small training set size. Therefore, events detected in the Tangshan data set studied by Lv et al. (2022) were also added to the training set, which enriches the set of positive samples (earthquakes). This is reasonable because the earthquake signals share similar morphology and similar energy distribution with similar pre-processing. Data augmentation was employed to further increase the number of positive samples by the time-shifted slicing method (Zhu et al. 2023). Specifically, we fixed the time window length to 512 pixels (5.12 s) and shifted the sampling window along the time axis to generate new earthquake signal samples. We ensured that the first arrival time of the earthquake signal was always within the time window. By using this method, we were able to obtain multiple effective samples from a single earthquake event. For negative samples (noise), we selected time windows that do not contain any earthquakes. The selection of these noise windows was based on the network catalogue, operation log file of the vibroseis truck, and Li's catalogue, which were used to exclude known events. Additionally, the selected windows were manually checked. We noted the existence of a lot of noise with infinite apparent velocity in the data set, which was rarely reported in the Tangshan data set, and we inferred that it was caused by the vibration of the container lab that housed the interrogator unit. We added such vibration noise into the training set, so the network can learn the characteristics of this noise to recognize it. In total, the training set was constructed from 30 Tangshan earthquakes, 4 Brady catalogue earthquakes on March 14 and 600 noise samples. The data set was divided into a training set and a validation set with a ratio of 4:1. The fully trained ADE-Net2 achieved an accuracy of over 99 per cent on both sets, which indicates that it has an excellent ability to detect earthquakes of both the Tangshan and PoroTomo data sets.

Detection results

During the experiment, a vibroseis truck (T-Rex) was operated at the site, which excited strong swept-frequency signals of 20 s duration with three vibration modes at each site for seismic tomography (Parker et al. 2018). Since the active source signals share similar features with earthquake signals, they are also detected by ADE-Net2 (Fig. 4a). Seven sweeps were missed by ADE-Net2 and no clear common feature was observed. Advances in ADE-Net2 visualization and interpretability may help to understand such rare outliers. The detected active source events were filtered out according to the operation log file.

ADE-Net2 detected a total of 90 earthquakes, of which 69 were reported by Li & Zhan (2018) and 21 were newly detected (Fig. 4b, Fig. 5a). The other 43 events in the Li & Zhan (2018) catalogue were missed in our detection result, including 10 events that overlap with the active source signals and others that are too weak. Considering the union set of earthquakes detected by both methods as the complete catalogue, ADE-net2 attains a recall rate of 67.7 versus 84.2 per cent for Template Match Filter (TMF). These results indicate that TMF has a stronger detection capability, whereas ADE-Net2 exhibits lower requirements in terms of template variety.

ADE-Net2 can also detect the signals related to T-Rex operation (Fig. 4c). They have a similar energy distribution pattern as the active source sweeps indicating similar location, but the amplitude is much lower and their origin time is about 10–15 s after the vertical mode or about 20–30 s after the longitudinal mode operation time listed in the log file. Therefore, we speculate that they were mechanical signals when T-Rex changes between modes of vibration. Besides, traffic noise by heavy trucks can lead to false detection by ADE-Net2, which occurred about 10–15 times per day. However, we found that by removing channels 3000–3500 that records record strong traffic noise, most of such traffic-related false detections were eliminated, but the detection of earthquakes was not affected.

ARRIVAL PICKING AND LOCATION

Arrival picking

The STA/LTA method (Allen 1978) has been widely used to automatically pick arrivals, and we used this method to pick the first arrival on each trace (Fig. 6a). The arrival was picked at the peak on the STA/LTA trace (Figs 6a and b). Fig. 6(b) shows the picking result of an M 0.05 catalogue earthquake. However, the picked arrivals show strong variations across the whole array and the picking scatter is up to 0.5 s, which is likely due to the low SNR on records of part of the array. Since all catalogue earthquakes that we used occurred to the south of the DAS array, leading to a large azimuth gap, such large picking scatter causes relocation error up to kilometres due to a strong trade-off among origin time, depth and epicentre. To improve location precision, quality control on raw picks is needed.

To remove the outliers, we tried to separate picks into clusters. The Density-Based Spatial Clustering of Applications with Noise (DBSCAN) method developed by Ester et al. (1996) was adopted. This method divides data points into specific groups based on the density of the continuous region. The operation starts with a picked point on the time-channel figure (Fig. 6b), and if there are more than ρ points located within a distance R, the point will be set to be a core point and all points within distance R are pulled out as border points. In our case, the distance is defined as sqrt((Δx/dx)^2+(fs*Δt)^2), where Δx is the space difference, dx is the channel spacing, Δt is the time difference and fs is the sample rate. If one border point gets more than ρ neighbours within a distance R, it will become a new core point. This process is repeated until all points have been visited, then the points that are neither a core point nor a border point will be identified as outliers. The R value is chosen based on trials. Fig. 6(c) shows an example in which each cluster includes at least 100 points. Aside from the outliers, clusters far away (>1 s) from the average arrival time were often found (Fig. 6b), which were also removed. As shown in Fig. 6(d), the standard deviation of the picks decreases from 0.775 to 0.296 s after removing outliers and the excess Kurtosis of the distribution (Fig. 6d) decreases from 3.68 to −0.43.

Location

We located the earthquakes by comparing the similarity of the picked arrival time curves with theoretical travel time curves using a grid search algorithm based on the residual shown in eq. (1), which is similar to the Geiger (1912) method. We found that the picking results are highly consistent with the synthetic traveltimes, which indicates the picks are P-wave arrivals or a shallow converted phase. A 1-D P-wave velocity model (Fig. 7d) was generated by interpolating the average velocity model obtained from local earthquake tomography (Foxall 2016) for location. The search area is 10 km by 10 km (Fig. 7a), with a grid spacing of 50 m.

where N is the number of picks used for location, $T_{o_i}$ is the observed picked arrival time, $T_{c_i}$ is the theoretical arrival time, and ${\overline{T}}_{c_i}$ and ${\overline{T}}_{o_i}$ are $T_{c_i}$ and $T_{o_i}$ minus array minimum value. However, the uncertainty of the location result is up to hundreds of metres using eq. (1). Therefore, we introduced another misfit function: scale misfit (SM, eq. 2).

SM reflects the difference between scales of theoretical arrival times and picked arrivals, which is sensitive to distance from the epicentre to the DAS array (Fig. 7). We designed an experiment to demonstrate SM's sensitivity to location (Fig. 7). In Fig. 7, theoretical arrivals at four possible epicentres (P1-4) at different distances are compared with the picks and the theoretical arrivals for the catalogue location (P0). The value of SM of each point in the whole search area also suggests that SM is sensitive to the epicentral distance.

Therefore, we chose to use a joint misfit function to reduce location uncertainty. The joint misfit function (eq. 3) is defined as the sum of individually normalized MSM and SM (eq. 4):

where:

The improvement of using the joint misfit is demonstrated in Fig. 8. In Fig. 8, the possibility of the hypocentre can be expressed as a function of the normalized misfit. For example, the joint possibility is defined as:

Fig. 8(a) shows relocation results for one catalogue earthquake to demonstrate the sensitivities of the three misfit functions. The large picking scatter results in significant location uncertainty with MSM, whereas SM provides a reasonable epicentral distance estimation but with an incorrect backazimuth. Combing MSM and SM yields a better result with less uncertainty, which is about 100 m away from the catalogue epicentre.

Several detected events were not located due to different reasons. Some detected event's signal is contaminated by the active source signal, which makes it difficult to pick the arrival. The clustering also reduces available picks, when the events with available picks less than 4000 after clustering were also discarded. Events with non-overlap of high probability regions of the MSM (90 per cent confidence level) and SM (80 per cent confidence level) misfits after location. In total, we obtained locations for 36 events (Fig. 9), with most of them occurring in the southern part. As shown in Fig. 9, most located events occur out of the array, which means the depth is not well constrained, so we only discussed epicenters in the following. Compared to the catalogue result, the four relocated catalogue earthquakes shift in the Y-direction, which indicates limited constraint in this direction.

DISCUSSION AND CONCLUSIONS

Pressure modelling

Cardiff et al. (2018) reported a correlation between the occurrence of earthquakes and the cessation of pumping at Brady. In the case of cessation of production pumping, an increase in pore pressure may be the cause of the earthquakes as a result of reduced normal stress on fault planes close to failure. As shown in Fig. 5, a rapid increase in seismicity emerged on March 14 when the pressure in well 56A-1 increased as a result of the shutdown.

The pore pressure change at the epicentre of each located event was estimated by utilizing a simplified model, as investigated by Cardiff et al. (2018), who assumed that fluid flow occurs dominantly along a single faulted zone within largely impermeable host rock at the Brady site. The simple model represents that pressure changes with distance and time, and can be calculated using parameters including long-term volumetric flow rate, fault zone storativity and density of the extracted fluid.

where W is the well function:

where Δp represents the change in pressure after t in seconds from when pumping has stopped for a point at a distance of r in metres from the pumping well. Q is the long-term volume flow rate of the pumping well, ρ is the fluid density and g is the acceleration of gravity. S is the fault zone storativity and T is the transmissivity of the fault zone. The latter two are the largest source of uncertainty in this function. Cardiff et al. (2018) fit the pressure data recorded in well 56A-1 during the test period to give estimates of T and S. We compared the predicted pressure changes with the catalogue locations and our relocation results for the four catalogue earthquakes that occurred on March 14. The location difference results in less than 20 per cent difference in pressure. As shown in Figs 9 and 10(a), the newly located events provide better coverage in both the temporal and spatial domains. The maximum pressure change is up to 157 kPa, which corresponds to the most northerly event. The lower bound is about 71 kPa, which is within the range of the typical published threshold (0.01–0.1 MPa) for earthquake triggering (e.g. Gomberg et al. 2001). However, it should be noted that this simplified model ignores many factors that may result in significant uncertainty of pressure change.

Location uncertainty

To evaluate location uncertainty, a series of tests based on synthetic arrival times was conducted. The theoretical arrival times of each search gridpoint were computed with the 1-D velocity model and Gaussian noise with σ = 0.1 s added. The location error is defined as the epicentre difference between the true location and the grid-search result and is shown in Fig 11. The results indicate that the location error generally increases with epicentral distance from the array centre. The location error includes contributions from picking error, velocity model error and systematic error due to array geometry. The location error shown in Fig. 11 includes picking error and systematic error since we do not know the true velocity structure, so it is a lower bound of real error. In the area of most of the newly detected events, the location errors are about 200 m, which will not result in significant error in pore pressure change computation. Such location errors indicate that it is difficult to identify the ruptured fault of an individual earthquake. However, most newly detected events occur at the edge of the area where most historic earthquakes have been located (Fig. 9). Since only the first arrivals were used to locate events, the uncertainty in the Z-direction is much larger than that in the horizontal plane. A similar result was also reported in a previous study (Piana Agostinetti et al. 2022), which recommended using both P- and S-wave arrivals to address this issue.

Detection methods and comparison

It is of interest to compare the detection results of our method to that of Li & Zhan (2018). We checked the raw data of several examples (Fig. 12). Earthquakes with high SNR can be detected by both methods (Fig. 12a), whereas earthquakes with very low SNR are only detected by the template matching method. For event detection, Li & Zhan (2018) calculated the detection significance as:

where CC is the mean cross correlation. All events with a mean absolute deviation (MAD) higher than 30 were detected in this study, whereas some of the weaker events are missed (Figs 12 and 13a). Such weak events are almost invisible on both STA/LTA diagrams and raw data. On the other hand, the ADE-Net2 detected 21 events that were missed by (Li & Zhan 2018, Fig. 12). The detection significances (MAD) of newly detected events are lower than 9 (Fig. 13), which indicates that the waveforms are different from templates. Presumably, the newly detected events are at a distance from the template events. However, some of the events detected by both methods are not close to the template events. This may indicate that the newly detected events are not only different in location but also ruptured with different focal mechanisms. Aside from weaker events, earthquakes that occurred near in time to the active source sweeps are mislabelled as active sources by ADE-Net2, whereas the template matching method can separate them from the active source sweeps (Fig. 12, only by TMF 2).

CONCLUSION

We designed a new neural network ADE-Net2 based on ADE-Net (Lv et al. 2022) for earthquake detection on continuous records. The neural network detected 90 events, of which 21 events were not reported by either the local permanent network (Foxall 2016) or a previous study with TMF (Li & Zhan 2018). The neural network detected 90 events, of which 21 events were not reported by either the local permanent network (Foxall 2016) or a previous study with TMF (Li & Zhan 2018). This result suggests that the detection capability of ADE-Net2 may be weaker than TMF, but it serves as a valuable supplement. Using the ADE-Net2 detection results as templates can increase the variety of templates, which could further improve detection capability. Using the ADE-Net2 detection results as templates can increase the variety of templates, which could further improve detection capability.

The STA/LTA method was used to pick arrivals in the detection time windows and the DBSCAN method (Ester et al. 1996) was applied to conduct quality control on the raw picks. Due to the large pick scatter, we proposed a new misfit function for locating epicentres combining the arrival time residual and a scale factor. Due to the large pick scatter, we proposed a new misfit function for locating epicentres combining the arrival time residual and a scale factor. In total, 36 events were located in the southern part of the field, partially outside the area with the most historical microearthquakes. By combing more measurements (e.g. S-wave arrival) on high qualify waveforms, the uncertainty of location is expected to be improved and it may be possible to provide depth information too.

The pore pressure changes at epicentres caused by the pumping shutdown were computed with a simplified model and the predicted values range from 71 to 157 kPa, which is larger than the typical earthquake-triggering threshold. In the future, the heterogeneity should be taken into account during the pressure change computation and the change of Coulomb failure stress on the fault plane could be used instead of the pore pressure change, which is expected to be more sensitive to earthquake occurrence.

ACKNOWLEDGEMENTS

We thank the team's effort in the field and ORMAT Technologies Inc. for providing site access. The authors thank the team for generously sharing this fantastic data set. This work is partly supported by the National Key Research and Development Program of China (2021YFA0716800) and the National Natural Science Foundation of China (U2239203).

CONFLICT OF INTEREST

Authors declare no conflict of interest.

DATA AVAILABILITY

The PoroTomo Team acquired the DAS data set used in this study (https://gdr.openei.org/submissions/980). The catalogue data were accessed through the Northern California Earthquake Data Center (NCEDC), doi:10.7932/NCEDC.

References