Building the fracture network model for the Okuaizu geothermal field based on microseismic data analysis

Understanding flow behavior in a geothermal reservoir is important for managing sustainable geothermal energy extraction. Fluid flow in geothermal reservoirs generally occurs in complex existing fracture systems in which the reservoirs are situated in highly fractured rocks. To simulate a discrete fracture model, the location and orientation of the fracture were computed using statistical processes and observational data. In many cases, estimating the location and orientation of fractures from 1D borehole logging data is challenging. In this study, we used microseismic data to build the fracture network systems and extract the detailed positions and dimensions of the fractures. We used the microseismic data recorded at the Okuaizu Geothermal Field, Fukushima Prefecture, Japan, from 2019 to 2021. First, we located the hypocenters, removing the effect of uncertainty in the velocity structure of the geo-thermal fluids. We relocated and clustered the seismic events based on waveform similarity. We analyzed each cluster to define the fracture orientation using principal component analysis (PCA) and focal mechanism (FM) analysis. We used the P polarity with the S/P ratio as a constraint for a better fault-plane solution. With PCA, we can extract the fracture dimension of each cluster. Our cluster analysis showed that the clusters were not always planar fractures, and we interpreted them as fracture zones. Based on the consistency between PCA and FM, each cluster/fracture zone was classified into three conceptual models to characterize the fracture network system in this field. This model showed variations in the orientation of small fractures within the fracture zone. We characterized the spatial variation in fracture distribution and orientations in the reservoir and demonstrated the fracture network system of this field. The fracture zone near the injection well has a N–S strike, and the dip is above 80°; however, the fracture zone in the northeastern part of the injection well has a NW–SE strike with a dip between 60° and 80°. The fracture network system estimated in this study is crucial for robust reservoir modeling because our model is more realistic, observa-tion-orientated, and includes local anomalies of reservoir properties.


Graphical Abstract 1 Introduction
Geothermal energy harnesses the natural heat emanating from the Earth's interior, unlike conventional energy sources such as fossil fuels or nuclear power, which rely on finite resources and degrade the environment.Geothermal energy offers multifaceted advantages.First, geothermal energy can provide a constant power source, if appropriately used, as Earth's heat is continuously produced.In contrast to solar or wind power generation, geothermal power generation is not dependent on weather or time of day.Additionally, geothermal energy is an environmentally friendly energy source that emits significantly fewer greenhouse gasses than fossil fuels.Geothermal power plants have a relatively small carbon footprint, with minimal impact on ecosystems.As the world seeks to transition toward greener and more sustainable energy solutions, geothermal energy will increasingly play an important role in the future, as we strive for carbon-neutral energy solutions (DiPippo 2012).
Geothermal fluid flows into existing permeable fractures in a geothermal reservoir.Many geothermal reservoirs consist of fractured rock, in many cases granite.These existing fractures have large variations in orientation, dimensions, and spatial distribution and are complexly interconnected, forming a fracture network.A fracture network explicitly represents the spatial distribution of fractures and their dimensions, orientations, and connectivity (Adler and Thovert 1999).Understanding the fracture network system in geothermal reservoirs is crucial for assessing and developing geothermal resources, as it serves as the primary pathway for geothermal fluids.Understanding the interactions and variations of fractures within fracture networks is fundamental for characterizing the fluid flow in fractured reservoirs and determining reservoir properties, such as bulk permeability (Berkowitz 2002;Zhao et al. 2011).
Various methods have been used to estimate fracture network systems, such as outcrop mapping, fracture geometries from borehole image analysis, stochastic modeling, and geomechanical modeling (Lei et al. 2017).Outcrop mapping and borehole imaging can provide accurate fracture distributions; however, fracture information from both methods provides a limited view of fracture in 1D or 2D (Wu and Pollard 2002).Stochastic modeling generates a 3D discrete fracture model (Dershowitz and Einstein 1988) based on a random process of location, orientation, and power law on fracture sizes (Davy 1993;Bonnet et al. 2001).Geomechanical modeling can add a preferential bias according to the state of stress (Dershowitz and Einstein 1988;Nick et al. 2011;Ishibashi et al. 2012).A 3D fracture network can also be generated by extrapolating 1D/2D geological data and partially utilizing observational data (Wu and Pollard 2002).However, all of these methods for 3D modeling rely on stochastic processes to some degree to generalize the spatial distribution of fracture systems, and the fracture network models using these methods cannot really employ local anomalies of fracture distribution.The fluid flow in the reservoir occurs very heterogeneously; therefore, using the fracture network model based on real observed information (realistic) is favorable.Obtaining a complete picture of a fracture network is challenging owing to the complexity of the geological processes involved in fracture formation.
To address this problem, we estimated a fracture network system by fully utilizing microseismicity analysis, which has the advantage of directly providing 3D information on each fracture in the reservoir (Koepke et al. 2020).Microseismic monitoring is typically used to monitor reservoirs in both hydrothermal and enhanced geothermal systems (EGS) cases.In general, microseismicity occurs owing to shear slip on existing fractures, regardless of the triggering process.We can estimate the locations of microseismic events and regard them as the locations of existing fractures (Shelly et al. 2013).Obtaining the precise location of microseismic events is important for estimating the real location, density, and orientation of fractures (Eaton 2018).Fractures can be extracted from seismic clusters or nearby microseismic distributions with similar waveform characteristics.We can then estimate the fracture orientation from the shape of the seismic clusters or directly from focal mechanism (FM) results.Thus, by appropriately utilizing microseismic data and analysis, we can build a fracture network that is more realistic and based on observational data and that can also employ local anomalies of fracture network characteristics.
This study focuses on the Okuaizu Geothermal Field (OGF), one of Japan's largest geothermal fields (Fig. 1).Several microseismic studies have been conducted on injection test activity (Okamoto et al. 2018(Okamoto et al. , 2020(Okamoto et al. , 2021)).However, the fracture system at the OGF is not well known, and only known fault system information (Kato et al. 2021).The hydraulic connection between wells and the effect of water injection as recharge are not well understood because of the lack of fracture network information (Okamoto et al. 2018(Okamoto et al. , 2020)).This study elaborates on microseismic data analysis to identify the fracture network system in the OGF.First, we located precise hypocenters by selecting proper seismic phases based on station correction and travel time residual error distribution.We then relocated and clustered microseismic events.We estimated each cluster dimension and orientation using principal component analysis (PCA) and compared the orientation with the FM results.We integrated all the analysis results and characterized the fracture network system of the OGF, consequently proposing a conceptual and equivalent model of the fracture network system.

Field overview
The OGF is located in Fukushima Prefecture, Japan, approximately 50 km west of the volcanic front of the Northeast Japan Arc (Fig. 1).The OGF started commercial operation with a capacity of 65 MW in 1995, which was the maximum capacity of a Japanese geothermal power plant.However, its capacity was reduced to 30 MW in 2017 because of steam depletion.The OGF reservoir was dominated by steam, with a temperature reaching 313 °C before the water injection test in 2015 (Okabe et al. 2016).The areas in the Chinoikezawa footwall fault and its southeastern fault have gradually become superheated because of a drop in pressure due to production activity.
Since 2015, Okuaizu Geothermal Co. Ltd. (OAG) has been conducting gravity-driven water injection tests in the reservoir to mitigate the reduction in steam production as part of the "Technology to Evaluate and Manage Geothermal Reservoirs" project of the Japan Organization for Metals and Energy Security (JOG-MEC) (Okabe et al. 2016).The total injected volume in 2015 was 1.7 × 10 5 m 3 with an average wellhead pressure of 0.17 MPa (Okamoto et al. 2018).The project was divided into two stages: Stages I and II.Stage I ran from 2013 to 2015 with the objective of locating the recharge well and initiating the injection experiment.Stage II ran from 2016 to 2019 with the intention of maintaining the injection, calibrating the simulation model, and creating a guideline to wrap up the project (Okabe et al. 2015).In 2019, the water injection test resumed, which injected 6.14 × 10 5 m 3 until 2022 as part of reservoir management (JOGMEC 2022).
The principal production zone occurs at a 1.0-2.6 km depth within fractured Neogene formations along two NW-SE trending faults, Chinoikezawa and Sarukurazawa Faults, with a dip steeply toward the north (Okabe et al. 2015) (Fig. 1).Borehole logging found steep northward-dipping fractures, and lost circulation occurred especially at the fractures with a NW-SE strike steeply north dipping, which is consistent with the characteristics of the Chinoikezawa footwall Fault (Kato et al. 2021).An unknown fault with a NE-SW strike trend was proposed based on drilling, logging, and tracer tests.This fault may hydraulically extend from the Chinoikezawa Fault and other surrounding areas (Kato et al. 2021).The OGF reservoir was dominated by steam, with a temperature reaching 313 °C before the water injection test in 2015 (Okabe et al. 2016).
The maximum horizontal stress (S Hmax ) trend of this region is WNW-ESE or NW-SE, and the stress state is suggested to be a strike-slip or normal fault type according to the Japanese crustal stress database (AIST 2024).However, the FM from the Japan Meteorological Agency (JMA) unified catalog around this region shows a normal fault-type stress state with a NW-SE strike direction (Asanuma et al. 2012).This information suggests that the stress state in this field has a strike-slip fault and a normal fault with an orientation of WNW-ESE to NW-SE.

Microseismic data
Since 2015, OGF seismic monitoring networks with sampling rates of 1000 Hz, consisting of four borehole stations and five surface stations, have been installed to cover the seismic activity during the water injection operation as part of the project.The surface stations (YAE1-YAE5) were buried several meters apart, and the borehole stations (YAE6-YAE9) were installed at 280-390 m depths (Fig. 1).Continuous seismic record data were transmitted in real-time to the Fukushima Renewable Energy Institute-National Institute of Advanced Industrial Science and Technology (FREA-AIST) server and stored there (Okamoto et al. 2018).FREA-AIST recently developed an automatic machine learning picker for geothermal microseismicity analysis (Okamoto et al. 2024).FREA-AIST employs a uniform velocity structure for P (Vp = 3900 m/s) and S-waves (Vs = 2250 m/s) for hypocenter determination (Okamoto et al. 2020).The same velocity model was used in this study.Hypocenters were routinely determined using FREA-AIST with manually selected P-and S-wave arrival times.
Several natural earthquakes occurred in the area prior to the water injection tests and were detected by the JMA network.The microseismicity activity did not correlate with the change in the production rate, even during the pressure build-up operation (Asanuma et al. 2011).Microseismic events generally occur in a minimal area between two NW-SE faults, Chinoikezawa and Sarukurazawa, and between two NE-SW faults, which are not exactly the same region as the main production zone (Okamoto et al. 2018;Kato et al. 2021).The microseismicity focal depth was mostly 500-2000 m below sea level (mbsl).Two significant seismic clusters were activated by the water injection test in the vicinity of the injection well and in the northeastern region of the injection well, which is associated with the Sarukurazawa Fault (Okamoto et al. 2018(Okamoto et al. , 2020)).
In addition to routine monitoring by FREA-AIST, several studies on microseismicity have been conducted at the OGF (Okamoto et al. 2018(Okamoto et al. , 2020(Okamoto et al. , 2021)).Timedependent microseismic activities occurred near the injection point and in the northeastern part of the injection well, some of which were correlated with operational activity (Okamoto et al. 2018).Geophysical imaging analysis has depicted the migration of the injected fluid, and its contribution to local microseismic activation has been elucidated (Okamoto et al. 2020).Okamoto et al. (2021) monitored the horizontal-to-vertical (H/V) spectral ratios and connected them to fluid activity that may trigger a major seismic event during water injection.However, previous studies have not focused on the fracture system in the OGF, and the details of the fracture system in this field remain unknown.Therefore, this study aims to build fracture network systems using microseismic data analysis from January 2019 to September 2021.

Hypocenter determination
Hypoinverse, which employs a nonlinear least-squares algorithm, has been widely used to determine the absolute location of earthquakes (Klein 2002).In this study, we first determined the absolute hypocenters using P-and S-wave arrivals that were picked manually by FREA-AIST analysts for routine analysis.A homogenous elastic wave velocity structure (Vp = 3900 m/s and Vs = 2250 m/s) for the OGF (Okamoto et al. 2020, p. 20) was used to determine the hypocenter.From all the events determined by Hypoinverse, we extracted reliable seismic events based on root mean square (RMS) of travel time less than 0.2 s, which corresponds to ~ 400 m spatial error, and the hypocenter depth < 5 km.We successfully located 8732 absolute locations from the 10,242 detected events (Fig. 2).All microseismic hypocenters show a cloud-like distribution and are generally denser near the injection well.
To compensate for the heterogeneity in real velocity structures, we evaluated the travel time residual and station correction in the hypocenter determination.The travel time residual of the P-and S-wave arrivals is the difference between the observed and theoretical travel times from an earthquake source to seismic stations.The travel time residuals can be positive if seismic waves arrive at a station later than the forecasted phase arrival, suggesting a delayed arrival due to the slow velocity anomaly in the ray path and vice versa.By selecting the appropriate phases after station correction and travel time residual checking, the effects of velocity structure heterogeneity can be prevented.We assumed that the P-and S-wave travel time residuals followed normal distributions with a mean value close to 0 (Fig. 3a).In fact, we found a second peak at station YAE9 for the residual P-and S-wave travel time between − 0.35 s and − 0.15 s (Fig. 3a), which can be caused by different velocity structures.However, the standard deviation of P-and S-waves at stations YAE2 and YAE4 exceeding 0.4 for both P-and S-waves indicates significant local heterogeneity.
Then, we calculated the station correction for P-and S-wave arrival times by Joint Hypocenter Determination (JHD) (Pujol 1992), where a systematic pick error was added for a group of events and tested to minimize the RMS during an iterative process.Station correction was used to remove uncertainty from various factors, including instrumental drift, site effects, noise, and velocity uncertainty in the ray paths.In this context, we assumed that the station corrections were mainly due to velocity structure uncertainty; therefore, they were effectively equivalent to the aggregated travel time residuals.We found that the station correction has a range value of ± 80 ms and ± 200 ms for P-wave and S-wave, respectively.The station correction of P-wave was consistent and had values around − 30 to 0 ms for the majority of the stations.The two exceptions were station YAE5, which had a station correction close to 70 ms, and station YAE9, which had a station correction of approximately − 85 ms.The station correction of S-waves was large for the surface stations (YAE2, YAE4, YAE5) and borehole station (YAE9), which were more than 100 ms and less than − 100 ms, respectively (Fig. 3b).
Travel time residual and station correction values serve a similar purpose, considering the principles and differences in the process of analyzing them individually or aggregately.We observed a consistent trend between the mean value of the travel time residual and station correction for almost all stations for the P-and S-waves.At station YAE9, both the travel time residual and station correction values were negative, suggesting the early arrival of the seismic waves.Both methods suggested significant uncertainty in the velocity model-related stations YAE5 and YAE9.Therefore, we did not use the phase data from stations YAE5 and YAE9.We also did not employ the S-wave phase, which had large station correction values exceeding 100 ms at stations YAE2 and YAE4.Station corrections for other stations were used as inputs for the hypoinverse.We improved the RMS travel time residual of the hypocenter determination from 0.057 to 0.046 s on average (Fig. 2d).

Waveform similarity evaluation
We used cross-correlation to remove potential bias from human detection in phase picking and improve the consistency of the detection phase arrival.Cross-correlation was used to evaluate the similarity between two waveforms of different events at the same station and provided a lag time between two events.The lag time, also known as the differential time, represents the difference in the P-and S-wave arrival times between the two waveforms.The cross-correlation value and differential time were used as input data for a more accurate estimation of the hypocenter relocation using a double-difference equation-based method (Schaff 2005;Trugman and Shearer 2017).
We calculated the cross-correlation value and lag time in the time domain (Deichmann and Garcia-Fernandez 1992;Beyreuther et al. 2010) for all possible pairs of events in each year at the same station and component.We computed P-and S-wave differential times using the vertical and EW components, respectively, as the S-wave signals tend to be clearer than the NS components.Cross-correlation for P-and S-wave were calculated in frequency 5-20 Hz, with a time window of 0.05 s before P-and S-wave arrival time and 0.2 and 0.3 s after P-and S-wave arrival time, respectively.We set the minimum cross-correlation value for relocation to 0.6, which is often used (Bulut et al. 2011;Kwiatek et al. 2014).Consequently, we obtained P-and S-wave differential times of approximately 3.4 × 10 6 pairs for P-wave and 3.9 × 10 6 pairs for S-waves.Figure 4 shows an example of the waveforms from two different clusters with high similarity in terms of the cross-correlation value.

Hypocenter relocation and clustering
The hypocenters were relocated using GrowClust (Trugman and Shearer 2017).GrowClust relocates event pairs based on the waveform similarity defined by the crosscorrelation value and differential time to identify them in the same event cluster simultaneously, maintaining fixed relative locations of these events.GrowClust uses the initial hypocenter, station, cross-correlation value, and differential time as almost the same inputs as another known relocation method, the double-difference method (hypoDD) (Waldhauser 2000).GrowClust integrates cluster analysis and event relocation within each cluster using a hybrid hierarchical clustering algorithm.This algorithm clusters events based on waveform similarity and relocates each event relative to others within its respective cluster.Unlike hypoDD, GrowClust relies exclusively on the differential cross-correlation times of pairs of events.The accurate location of clusters was crucial in this study to obtain the structure of each cluster.
In the clustering and relocation processes, the event pairs were set to have more than five differentials crosscorrelation results.We used the same 1D velocity model that was applied to estimate the absolute locations: Vp = 3900 m/s and Vs = 2250 m/s (Okamoto et al. 2020).The maximum station distance from an event pair was limited to 5 km.A seismic cluster is defined as the presence of at least five events.We set a maximum RMS travel time residual of 0.2 s after experimenting with other thresholds and following Trugman and Shearer (2017) and Skoumal et al. (2019), then finally obtained optimal results without losing more events.
We successfully relocated 5041 events (Fig. 5), about 58% of all located events.The percentage of events identified as clusters may be related to the complexity of the fracture system (Asanuma et al. 2007).For example, Basel, Switzerland case showed approximately 70% events (Asanuma et al. 2008;Häring et al. 2008;Mukuhira et al. 2013), and Soultz France cases showed 37% and 58% for shallow and deep reservoirs, respectively (Moriya et al. 2002(Moriya et al. , 2003)).In both fields, the microseismic clusters exhibited complex fracture systems.On the other hand, Cooper Basin, Australia case showed nearly 100% of events were identified as one cluster.Microseismic events occur in large single horizontal structures (Baisch et al. 2006;Asanuma et al. 2009).Owing to variations in geological conditions, there were different preexisting fracture systems in the Basel, Soultz, OGF, and Cooper Basins before stimulation.We believe that the variation in existing fractures makes the spatial characteristics of the fracture network systems different.Considering the relationship between the ratio of the clustered events and the complexity of the reservoir, the fracture system at the OGF was as complex as that of the reservoir in the Soultz case.
From the 5041 events relocated in this study, the relocation events encompassed 82 clusters (Fig. 6).The seismic cloud was tighter than the absolute locations, which were scattered near the injection well.Single events and small clusters (fewer than five events) did not relocate during the clustering.With these improvements, the RMS of travel time residual was also improved below 0.1 s (Fig. 5d).

Focal mechanism analysis
We estimated the FM with the first motion of P-waves and the S/P amplitude ratio using HASH (Hardebeck 2002(Hardebeck , 2003)).Since March 2021, FREA-AIST has begun to perform manual polarity picking of the P-wave first motion as routine analysis.However, the polarity of the data from January 2019 to February 2021 has not yet been determined.Therefore, we selected the P-wave polarity using a deep learning model (Uchide 2020) to reduce the processing time.The algorithm classified the polarity as 'up' and 'down' using the vertical component of seismograms around P-arrival times as input.
First, we checked the performance of the deep learning model (Uchide 2020) on our data from March to September 2021, which already contained manual polarity-picking data.The results of polarity picking by the deep learning model show satisfactory performance to the microseismic events with Mw ≥ 1, especially for surface stations YAE2 and YAE5, the precision is close to 0.8 (Supplementary file 1: Fig. S1-36).The precision was calculated as the number of correct predictions per accepted prediction.However, for microseismic events with Mw < 1, the precision is 0.5 on average for the borehole stations and above 0.6 for surface stations, except for station YAE4.Therefore, we focused on microseismic events with Mw ≥ 1 for polarity picking using a deep learning model.In total, we have 146 microseismic events with Mw ≥ 1 for FM estimation in data from January 2019 to September 2021.We did not include microseismic events with Mw < 1 because of the lower precision of polarity picking by the deep learning model and the lower signalto-noise ratio (SNR) owing to the smaller magnitude.
For the HASH application, at least four polarities were required, and the SNR between the P amplitude and noise level was > 3 (Hardebeck 2003;Skoumal et al. 2023).We determined noise levels by calculating the Fig. 4 Example of the waveform from clusters 5 and 9 at station YAE6.The waveform has been filtered using a frequency band of 5-20 Hz.The waveforms show a high cross-correlation value above 0.9.The blue and green vertical line shows P-wave and S-wave arrivals maximum absolute amplitude in selecting waveform windows − 0.5 to − 0.2 s relative to the P-wave arrivals.To calculate the S/P amplitude ratio, we filtered all waveforms of the three components between 5 and 30 Hz.We calculated the maximum absolute amplitude of P-and S-wave in waveform windows from − 0.1 to 0.1 s relative to the P-and S-wave arrival times.Subsequently, we computed the absolute amplitude of each P-and S-wave by summing the squares of their three components, and then taking the square root before calculating the ratio.
Seventeen FM solutions were successfully estimated from 146 microseismic events (Fig. 7).Although we obtained FM solutions for other events, we decided not to use them for further analysis because the fault-plane uncertainty was greater than 45°.Different fault types (strike-slip, normal, and reverse faults) were obtained; three of them had a fault-plane uncertainty below 35°, Despite their high uncertainty, the FM solutions presented in this study can still be used to estimate fracture orientation by combining the seismic cluster geometry.
For example, in cluster 3 (C3), which has an uncertainty below 45°, we obtained three consistent FM solution results that are also consistent with the PCA orientation, which is nearly N-S with a dip of approximately 90°.

Principal component analysis
We applied PCA to each microseismic cluster to estimate the cluster dimension (size) and orientation based on the hypocenter distribution (Michelini and Bolt 1986;Moriya et al. 2003;Vidale and Shearer 2006;Mukuhira et al. 2023).PCA is a statistical technique that helps to identify a dataset's principal components or axes of maximum variance.Using PCA, we calculated three principal components (values and vectors) representing the microseismic cluster dimensions as ellipsoidal.We modeled the ellipsoidal microseismic cluster using principal component vectors, and their lengths (principal component values) were increased by a factor of two.The defined ellipsoid should cover 95% of microseismic events in the microseismic cluster.The hypocenter is often a seismological quantity that can be determined with the greatest degree of accuracy.Thus, the cluster orientation estimated by PCA, which is used to estimate the fracture orientation, can stand alone, although no FM information is available.
We target clusters that have a minimum of 20 events for a detailed orientation evaluation because it is difficult to estimate a cluster's shape for a small population cluster.In this step, we obtained 30 of the 82 clusters for further analysis to estimate the fracture network system (Figs.8 and 9).We calculated the contribution ratios of the three principal components (Fig. 8a) to obtain an overview of the 3D aspect ratios of the cluster.The cluster should be planar if the third principal component is close to zero, as is the case in Basel (Asanuma et al. 2008;Mukuhira et al. 2017Mukuhira et al. , 2021Mukuhira et al. , 2023)).The contribution ratio of the three principal components shows that our cluster shape is mostly not planar but ellipsoidal, as the contribution ratio of the third principal component was 0.175 on average in the range of 0.05-0.3.A cluster was defined as a fracture if the third principal component ratio was < 0.1.Because we cannot disregard all the clusters assessing with the third principal component ≥ 0.1, we introduced another group as a fracture zone if the third principal component contribution ratio is between 0.1 and 0.3.We defined the fracture zone as comprising small fractures that could have the same orientation as the fracture zone or a different orientation (Fig. 10).
Of the 30 clusters, we identified only one fracture (cluster 8) and 29 fracture zones.We also estimated the strike and dip of the fracture/fracture zones using third principal component vectors.The third principal   (Trugman and Shearer 2017).Each cluster has at least five events.Clusters 1 to 5 have more than 300 events in two dominant areas near the injection well and in the northeastern part of the injection well component is the normal vector of the fault, and the strike can be converted from the normal vector of the fault (Shearer 2009).The strike orientation was generally NW-SE and NNW-SSE, with a dip > 60 • (Fig. 8b).The first large cluster, with a total of 1150 events, located in the northeastern part of the injection well, had a strike orientation of NW-SE with a dip of 82 • .The other four large clusters with more than 300 microseismic events have strike orientation NNW-SSE with a dip > 70 • and are located either near the injection well or the northeastern part of the injection well.

Conceptual fracture network model
We built a fracture network system model using the results of the microseismic data analysis, which can provide 3D information on the size, distribution, and orientation of the fracture system.We propose three conceptual models to characterize the fracture network system in this field using the FM and PCA results (Figs. 9 and 10).The proposed model showed variations in the orientation of small fractures in the fracture zone.By representing each fracture zone (cluster) using these conceptual models, we can characterize the spatial variation of the fracture distribution in the reservoir and consequently model the fracture network system of the entire reservoir.Descriptions of the three conceptual models are as follows.
Model 1 is an equivalent fracture system containing small conjugate fractures that are evenly distributed in both directions (Fig. 10b) with respect to the orientation of the principal cluster direction, which is represented with PCA 1st component orientation of the PCA.Model 1 is based on the fracture mesh model proposed by Hill (1977) and Sibson (1996), which is often used to describe the fracture system of geothermal reservoirs (Evans et al. 2005;Häring et al. 2008).The microseismic distribution within the fracture mesh model is considered to extend in the general direction of the maximum principal stress (Hill 1977;Häring et al. 2008;Mukuhira et al. 2023).However, the stress state in this field is somewhat uncertain between the strike-slip and normal fault types.We slightly changed the concept of the fracture mesh model and defined the small fractures as conjugate fractures with respect to the cluster orientation, that is PCA 1st component orientation.Because most PCA 1st component directions are consistent with the range of S Hmax (NW-SE and N-S), this is not an unreasonable assumption.
Model 1 is the base case of our fracture network model, which utilizes the PCA results to define the conjugate fracture orientation.We can update Model 1 to Model 2 or Model 3 using FM information to update one of the conjugate fracture orientations, which is more likely to exist.
Model 2 is the predominant fracture system, consisting of small fractures with nearly identical orientations (Fig. 10b).These small fractures are conjugate fractures with respect to the orientation of the PCA 1st principal component.The fracture zone is classified as Model 2 when the FM result is consistent with one of the conjugate fractures with respect to the orientation of PCA 1st principal component of the PCA.We place more weight on the probability of small fractures in this direction, which are the dominant small fractures.Small subsidiary conjugate fractures to the cluster orientation also exist.
Similar to Model 2, Model 3 is the predominant fracture system (Fig. 10b).The dominant small fractures can be described as fractures in one direction, and subsidiary conjugate small fractures to the cluster orientation also exist slightly.The fracture zone is classified as Model 3 when its orientation of the fracture zone estimated by PCA 1st component is consistent with that of the FM solution of the seismic events included in the fractured zone.The consistency between these results places more weight on the probability of a small fracture occurring in one direction.

Fracture network system characterization
We classified 21 fracture zones as Model 1 (Table 1).Each fracture zone in Model 1 had different strike orientations estimated by the PCA 1st component, mostly NW-SE and N-S (Fig. 8b).We found that the fracture zones of Model 1 located near the injection well had a strike orientation of nearly N-S with a dip of 80°-90°.In contrast, the fracture zones of Model 1, located in the eastern part of the injection well, have a NW-SE strike orientation with a dip ranging from 60° to 80°.We consider the fracture zones in the eastern part of the injection well to be an extension of the Sarukurazawa Fault, as both share a similar NW-SE orientation.
Under specific conditions, a fracture zone can provide FM solutions for different fault types and orientations.Fracture zone C1 (Fig. 7) in the northeastern part of the injection well is an example of exceptional Model 1 with multiple FM solutions.This fracture zone is more complex than other fracture zones that have different fault mechanisms and orientations, although the FM reliability remains an open question.With a variety of FM solutions, we still assume fracture zone C1 as Model 1 because adjusting the orientation of the conjugate fracture without a dominant type of FM solution is challenging.However, we can update the fracture network model based on the FM solutions if additional information is available to support different types of models.
Model 2 represents eight fracture zones (Table 1).Most of the fracture zones in Model 2 were located in the eastern part of the injection well.These fracture zones have a strike orientation of NW-SE and a dip range of 80°-90 o (Figs. 7, 8b).In addition, two fracture zones were located in different areas: fracture zone C7 was located near the injection well and fracture zone C10 was located in the southeastern part of the injection well, which had different strike orientations (Table 1).We observed that the fracture zones of Model 2, especially those located in the eastern part of the injection well, exhibited nearly identical strike orientations and dips, as determined by PCA.
We identified the fracture zones C3 and C28 as Model 3 (Table 1).These fracture zones showed consistency between the FM solution and the PCA 1st principal component.Fracture zones C3 and C28 were located near the injection well and in the eastern part of the study area, respectively.Both fracture zones show different fault types and strike orientations derived from the FM (Fig. 7).We found that the fracture zone in Model 3 had almost the same strike orientation and dip as those in the previous model, which was located near the injection well for fracture zone C3 and in the eastern part of the injection well for fracture zone C28.
In the study area, several local faults exist with orientations of NW-SE and NE-SW, which influence the microseismic distribution.Consequently, we categorized 29 fracture zones in the OGF estimated from 29 clusters comprising at least 20 events.Among these, we identified 21 fracture zones classified as Model 1, six as Model 2, and two as Model 3. Notably, two distinct groups of fracture zones were observed: one situated Fig. 10 Proposed conceptual model of the fracture network system a An example of a map view of a fracture network system for clusters having more than 50 events (to get a better view of small fractures illustration).Small fractures (yellow line) are just an illustration, not give a real size, and b proposed model of fracture network system near the injection well and the other in the eastern part of the injection well.These two groups exhibited different strike orientations and dips, as estimated using PCA.We found that the fracture zones near the injection well had a strike orientation of nearly N-S, with a dip of 80°-90°.We interpreted that the fracture zones near the injection well work as vertical structure zones and are not related to known local faults.Considering that the water injection pressure (gravity driven) was too low to form a new structure, a vertical structure zone should have existed near the injection well before injection.The fracture zones in the eastern part of the injection well are mostly NW-SE in orientation, with dips ranging from 60° to 80°.
The fracture zones are correlated with the Sarukurazawa Fault.

Implication of fracture network system for reservoir modeling
Building a realistic fracture network system is challenging for the geothermal community and other relevant fields (Lei et al. 2017).An ideal fracture network model derived from a microseismicity analysis should accurately represent the location, dimensions, orientation, density, and connectivity of each fracture within the reservoir.There are only a few cases in which such fracture network systems have been delineated by microseismic analysis, such as in Basel (Häring et al. 2008;Mukuhira et al. 2013).We used a similar approach to estimate the location, orientation, and dimensions of the fracture zone.However, our microseismic analysis encountered technical or inherent limitations associated with the microseismic data or monitoring conditions.Our relocation results show ellipsoidal microseismic clusters rather than planar fractures.Therefore, we introduced the concept of fracture zones.Additionally, uncertainties in the FM make it difficult to employ information from the FM to build the fracture system model.However, despite these technical limitations, our microseismic data analysis still offers advantages for characterizing a fracture network system in an OGF reservoir.We successfully built a fracture network system and estimated the detailed location, dimensions, and orientation of the fracture zone and the small fracture variation of each fracture zone based on microseismic data analysis.The number of events in each fracture zone indicated a small fracture density within the fracture zone.
Compared to other conventional methods for estimating fracture network systems (e.g., geological mapping, borehole imaging, stochastic modeling, and geomechanical modeling) (Lei et al. 2017), our approach can provide a more detailed and physics-based fracture network system.Information on the dimensions, orientation, coordinates of each fracture zone, and conceptual model variation can be used to build a fracture network model that can realize heterogeneity in the spatial and orientational distribution of fractures for reservoir modeling.Furthermore, the incorporation of fractures allows for a detailed examination of the individual contributions of fractures to the fluid flow (Berkowitz 2002;Ishibashi et al. 2019).

Conclusions
In this study, we determined the precise hypocenters by selecting P-and S-wave phases, considering the velocity structure anomaly, phase correction by cross-correlation, and the latest clustering and relocation techniques.Using detailed locations and clustering results, we evaluated the shape of the clusters using PCA and found that most of the clusters were ellipsoidal.We integrated the cluster analysis results and FM to characterize the fracture system of each fracture zone.The main findings of this study are summarized as follows: • We found that the fracture network system of the OGF is very complex, as shown by the diverse orientations of the fracture zones, FM, and aspect ratios of the seismic clusters.However, the fracture system is still considered as an extension of a known fault system.
• We proposed three conceptual models with different variations of small fracture distribution that comprised the fracture network and characterized the fracture network system of the OGF.• The fracture network system of the OGF has spatial anomalies in its characteristics, such as the orientation or variation in fractures.Fracture zones near the injection well strike N-S and dip above 80° (this area has three of the fracture network models), and those in the northeastern part of the injection well have a dominant strike NW-SE and dip between 60° and 80° (Models 1 and 2 are found in this area).
As shown in this study, microseismic data analysis can be beneficial for estimating fracture location, orientation, and detail in fracture systems compared to other conventional methods, which are mainly based on stochastic models.Microseismic events provide observational physical data on fracture networks when analyzed appropriately.Recent seismological analysis techniques, such as new relocation and clustering techniques (GrowClust) and machine learning P-wave pickers for FM estimation, can enhance the utilization of microseismic data for fracture network system estimation and other purposes, although careful employment is necessary.Our approach is considered as a standard procedure for building fracture networks based on microseismic analysis.

Fig. 2
Fig. 2 Hypocenter distributions by Hypoinverse (Klein 2002).a Map view, b RMS of travel time residual (observed-predicted), c A-A' cross-section, d B-B' cross-section, and e C-C' cross-section

Fig. 3 a
Fig. 3 a Histogram of travel time residual of P-and S-wave arrival with mean and standard deviation from hypocenter determination results by Hypoinverse (Klein 2002), b station correction values for P-and S-wave arrival estimated with JHD for all stations (Pujol 1992), and c number of picking phases of P-and S-wave arrival used for station correction calculation

Fig. 5
Fig. 5 Relocation and clustering result by GrowClust (Trugman and Shearer 2017).Events are color-coded by the cluster ID. a Map view, b RMS of travel time residual (observed-predicted), c A-A' cross-section, d B-B' cross-section, and e C-C' cross-section

Fig. 7
Fig. 7 Focal mechanism solutions with hypocenter distribution.The color hypocenter distributions indicate clusters that have focal mechanism solutions, and the gray color of the hypocenter indicates hypocenters without focal mechanism solution

Fig. 6
Fig.6Distribution of number of events per cluster obtained by GrowClust(Trugman and Shearer 2017).Each cluster has at least five events.Clusters 1 to 5 have more than 300 events in two dominant areas near the injection well and in the northeastern part of the injection well

Fig. 8 a
Fig. 8 a Contribution ratio of cluster shape from three components PCA for clusters having microseismic events more than 20, and b lower hemisphere plot of PCA third component (showing by circle plot) and pole plot of focal mechanism solution [showing by cross (X)]

Fig. 9
Fig. 9 Three-dimensional ellipsoidal plot of fracture network system.The color ellipsoidal indicates a different model of the fracture network system

Table 1
Summary of the fracture system information for each cluster