Freely available online through the SEG open-access option.

ABSTRACT

The forecasting and risk assessment of induced seismicity associated with fluid injection have considerable importance for subsurface energy development. We have developed a seismic evaluation method called the possible seismic moment (PoSeMo) model to assess the potential seismic moment that could be released in the future based on current seismic activity. The PoSeMo model assumes the existence of a representative parameter that can describe the seismic characteristics of a given field. This parameter is defined as the seismic moment density, which quantifies the seismic moment able to be released per rock volume. The rock volume presumed to be in critical condition because of stimulation is defined as the stimulated rock volume. The current stimulation condition for the PoSeMo model can be estimated from the product of these two parameters. The difference between the output of the PoSeMo model and the observed cumulative seismic moment corresponds to the cumulative seismic moment that could be released in the future. This value can be transformed into the possible maximum magnitude that has clear physical meaning and that can be used as feedback on the stimulation operation for seismic hazard assessment. We have applied this model to a microseismic data set from the Basel engineered geothermal system project. We have successfully estimated reasonable values for seismic moment density and stimulated rock volume. The PoSeMo model performed well, and it provided reasonable estimates of seismic moment. The maximum magnitude estimated by the PoSeMo model was almost identical to the largest event that had occurred previously. Thus, it was concluded that the PoSeMo model satisfactorily demonstrated its feasibility as a real-time seismic evaluation method, based on physical parameters derived from microseismic information.

INTRODUCTION

Injection-induced/triggered seismicity has become a crucial problem associated with the global growth of shale gas/oil production, enhanced/engineered geothermal system (EGS) development, and carbon capture and storage, where hydraulic/fluid stimulation is implemented (Evans et al., 2012). Such activities are considered to increase the risk of seismic hazard in the regions of active underground development (Ellsworth, 2013). Therefore, an increasing number of studies have addressed the development of methods to mitigate this seismic hazard, e.g., the “traffic light system” and to forecast the seismic risk, e.g., the seismogenic index (Shapiro et al., 2010).

In this area of research, seismostatistical approaches have become increasingly important (Bachmann et al., 2011; Mena et al., 2013; Hallo et al., 2014; Mignan et al., 2015) because seismostatistical models are principally designed for forecasting future seismicity. In a typical seismostatistical model, e.g., the RJ model (Reasenberg and Jones, 1989) and the epidemic-type aftershock sequence model (Ogata, 1999; Hainz and Ogata, 2005), which were used in the above studies, the input data comprise mainly occurrence time (occurrence rate) and magnitude, taken as the b-value from the Gutenberg-Richter relationship (Gutenberg and Richter, 1942). Therefore, it is clear that microseismic information is not fully used.

Several recent studies have addressed the estimation of the maximum seismic magnitude from a combined physical and statistical perspective, which have inspired the model proposed in this study. The seismogenic index (Shapiro et al., 2010) provides the maximum magnitude of fluid-induced seismicity by combining information on the progress of the microseismic cloud modeled by the diffusion process and a and b values. Their subsequent work (Shapiro et al., 2011) extended their original model and proposed a relation between the maximum magnitude and principal axis length of the ellipsoid, which represents the seismic cloud. In their series of studies, the propagation of pore pressure was modeled with a diffusion process; therefore, it was reasonable to model the seismic cloud as an ellipsoid. Using a model that was more physically orientated, McGarr (2014) finds a bond between the cumulative seismic moment and the volume of injected fluid, based on elastic deformation theory combined with the Coulomb failure criterion and moment tensor. By introducing the Gutenberg-Richter relationship, McGarr (2014) establishes that the maximum seismic moment release was also described as a function of the volume of injected fluid. Kwiatek et al. (2015) extend McGarr’s model and they propose that the maximum magnitude could be expressed as a function of the volume of the formation weakened by the injection and the average increase in pore pressure. Although McGarr (2014) does not expand on the hydraulic processes behind his model, it could be presumed that the pressure perturbation of the diffusion process was invoked for the conditions of his model. Therefore, it was entirely reasonable that Kwiatek et al. (2015) use the ellipsoid volume to represent the seismic cloud.

The review of these earlier studies revealed some common points of interest. For example, information related to the stimulated rock volume or fluid volume and the magnitude (i.e., the b-value) is linked, whereas the theory/concept to integrate this information and to derive the magnitude is different. It appears universal that the seismic moment released by induced seismicity due to fluid injection can be correlated with information related to the stimulated rock volume/injected fluid volume and to information related to the magnitude, such as the b-value.

Following the recent increasing demand for shale gas/oil development, the processing methodologies and monitoring systems of microseismicity have progressed rapidly, enabling precise determination of hypocenters and source parameters in real time, which can be exploited by procedures for the forecasting and assessment of seismic risk. Based on experience obtained from a large number of events from various fields and on lessons from previous research, this study proposed a physically based method to access seismic risk using physical parameters from microseismic information pertaining to the stimulation.

POSSIBLE SEISMIC MOMENT MODEL

Basic concept of the possible seismic moment model

The current research group of authors has worked on microseismicity in many EGS/hydrothermal fields (e.g., Soultz, Cooper Basin, Basel, and Nishiyama), and its experience is that the rupture area of the largest event is constrained by the seismically active region, as has been reported by other researchers (Baisch et al., 2009; Shapiro et al., 2010,, 2011; McGarr, 2014). However, there is a possibility that critical geologic faults within the stimulation area have much larger rupture areas than indicated by the seismic cloud, which could result in hazardous earthquakes (Gischig, 2015). Furthermore, we have observed in the study for Basel (Mukuhira et al., 2013) that the seismic cloud is extended most just after shut-in, and that many larger events occurred simultaneously within the extended region. Based on these observations, we have established empirically that the maximum seismic moment or magnitude can be constrained by the volume of the seismic cloud, and that the occurrence of large events is related to the expansion of the seismic cloud. Therefore, based on this experience and on these observations, the proposed model is intended to explain the seismic moment as a function of the seismic cloud.

Before the introduction of our model, we first explain the assumption of a single parameter that characterizes seismic activity within a given field. This parameter is defined as the seismic moment density D (Nm/m3 or Mo/m3), which represents the seismic moment that could be released per rock volume for each specified field. We also define the rock volume in which pore pressure has increased sufficiently to trigger seismic events as the stimulated rock volume Vstim (m3). Then, the seismic moment density D and stimulated rock volume Vstim can be correlated as in equation 1:  
ΣMo=Mopossible=D×Vstim,
(1)
where Mopossible (Nm) is the possible seismic moment (PoSeMo) indicating the cumulative seismic moment that could be released from the entire stimulated rock volume. This is the definition of the PoSeMo model. Therefore, by estimating appropriate values for D and Vstim from the microseismic data, we can obtain Mopossible. Then, Mopossible is compared with the observed cumulative seismic moment, and the difference between them corresponds to the seismic moment that could be released in the future, which is finally translated into a maximum magnitude.

Fluid injected into fractured rock in the process of hydraulic stimulation, as is the case with geothermal or shale reservoirs, propagates within the permeable flow path as fluid flow. Then, the increase in pore pressure due to the injected fluid induces shear slip on existing faults. The PoSeMo model attempts to estimate the volume of rock in a critical state, meaning that pore pressure has risen sufficiently to trigger induced seismicity from microseismic migration. In combination with the index representing seismic activity per rock volume, the total seismic moment that could be released from the current microseismic distribution can be estimated.

Our model is designed to perform in semi-real time using on-site microseismic information, and it is intended to evaluate seismic activity over the short term (approximately three days) and middle term (approximately several weeks). This model can provide quantities to operators and stakeholders that are intuitively comprehensible for decision making in semi-real time. Because our model is based entirely on physics and does not use statistical parameters, its output would be helpful as an alternative point of view regarding probabilistic seismic hazard analysis.

Methodology

The parameters of the PoSeMo model: Mopossible, D, and Vstim are all estimated from the microseismic information. As a prior condition, the hypocenter and magnitude of the seismic events must be estimated reliably. The methodology is described as follows.

Seismic moment density D

We attempted to estimate D from microseismic events that occurred from one local area of the reservoir. Figure 1a shows a schematic of the sample region and the seismic events used to estimate D. This sample region is one local part of a larger reservoir and it has sufficient size to contain the hypocenters of several seismic events. It is accepted based on logging analysis at KTB, a German scientific deep drilling project (Ito and Zoback, 2000), and microseismic analysis at Basel (Mukuhira et al., 2013) that various sizes of fractures with various orientations with respect to the in situ stress exist within the rock mass. Some of these fractures can experience shear slip if the imposed stress state meets the Coulomb failure criterion under pressurized condition. Note: We assumed that the seismic moment is released at the point of the hypocenter (point source model), even though the fault area is shown in Figure 1. Therefore, we extracted the events having shear slip within the sample region, and then we estimated D by summation of all the seismic moments released from the sample region and the volume of the sample region using equation 2:  
D=ΣMogrid/Vgrid.
(2)

In practice, the microseismic cloud at any given time can be divided into uniform grids and one of the grids chosen for the estimation of D. In the case study of the application of this method to field data from Basel, Switzerland, described later with reference to Figures 2 and 3, the chosen grid is shown in Figure 4a. The D value can depend greatly on the concept behind the selection of the grid for the estimation of D and on the length of the grid. For example, the chosen grid could be the one in which the maximum number of seismic events has occurred before the given time. The meaning for D could be interpreted as the best sampling of D because D is estimated from the densest grid. In the same manner, Mopossible calculated from this D could be interpreted as the best sampled Mopossible in terms of the number of samples. Hereafter, we refer to this as Mopossible (Max num). The concept used to determine D should correlate with the physical meaning of the output of this model. In another example of the concept used to determine D, the grid that includes the largest event by the given time could be chosen and this D could be interpreted as the worst scenario D, and therefore, the Mopossible based on this D could be considered the worst-case scenario.

Stimulated rock volume Vstim

In hydraulic stimulation, the microseismic events are triggered by the increase in pore pressure. We considered that pore pressure propagates as fluid flows throughout the permeable fracture network (Watanabe et al., 2009), rather than using the diffusion model, as used in other related studies. During or just after stimulation, the increase in pore pressure should be the dominant mechanism for triggering induced seismicity. Therefore, it should be noted that other possible trigger mechanisms, including thermal stress and the poroelastic effect, were not considered here. This is because the contribution to the occurrence of induced seismicity by these other mechanisms might be much smaller than that from the pore pressure increase in the short term, although these other mechanisms can be the greatest contributors in the long term (Jeanne et al., 2015; Kwiatek et al., 2015). Under this assumption, the progress of the seismic cloud can be presumed to reflect the propagation of pore pressure. Therefore, Vstim can be estimated from the temporal migration of microseismic events.

First, the study area was divided into grids, as in the case for the estimation of D, but the grid size for the calculation of Vstim does not have to be the same size as for the D calculation. Referring to the hypocenter distribution of the microseismic events that occurred by the given time, those grids that included at least one seismic event were counted as shown in Figure 1b. Then, Vstim was calculated by multiplying the number of grids n by the volume of the unit grid Vgrid as described in equation 3:  
Vstim=Vgrid×n,
(3)
where n is the number of grids in which at least one seismic event occurred. Here, appropriate grid size is also critical for the Vstim estimation. The occurrence of at least one seismic event means that the pore pressure increased to a certain level at the hypocenter of that event. Moreover, it suggests that the pore pressure around the hypocenter increased, or specifically, the pore pressure in the permeable fracture around the hypocenter increased. Therefore, it is important to set an appropriate grid size based on this concept. For example, consider the case in which a large grid includes two seismic events that belong to different fluid path systems. When one of the events occurs, the grid could be regarded as stimulated, even though the other event has not occurred. This would lead to an overestimate of the stimulated rock volume and information on the other flow path system would be lost. In such a case, an inappropriate grid size can lead to the overestimation of the stimulated rock volume. Therefore, a grid size that can detect single seismic events would be favorable.

Evaluation of seismic activity

To evaluate seismic activity, Mopossible is compared with the cumulative seismic moment, which is referred to as ΣMo (observed), calculated from all the seismic events that occurred by the given time. When Mopossible is larger than ΣMo (observed), the difference between them corresponds to the cumulative seismic moment that could be released in the future. If this difference is released in a single shear slip, this difference represents the maximum PoSeMo, which is eventually translated into the moment magnitude Mw. Meanwhile, when Mopossible is smaller than ΣMo (observed), it can be said that the model fails to replicate the seismic activity.

It should be noted that this model is unable to determine when the differential seismic moment may be released (occurrence time) or how it may be released (magnitude frequency distribution). Therefore, the maximum Mw estimated by this model stands under the assumption that the differential seismic moment is released by a single shear slip. It must be highlighted here that the appropriate overestimate of the seismic moment is important to obtain the plausible maximum Mw. The PoSeMo model derives the seismic moment that should be released when considering D. Therefore, an overestimate, which is the discrepancy between the model and observed values, can be interpreted as the residual seismic moment able to slip.

Comparison to relevant studies

The PoSeMo model has some analogies with the seismogenic index (Shapiro et al., 2010) and McGarr’s theory (McGarr, 2014), including its different versions (Kwiatek et al., 2015) in terms of using magnitude information and the extension of the seismic cloud. However, these three models are based on entirely different concepts with respect to the evaluation of the stimulated volume, magnitude information, and theory used to estimate the maximum magnitude, which are summarized as follows: It is meaningful to compare the results of three models and this is addressed in a later section.

  1. All three models include the stimulated volume or volume of formation weakened by injection, which are estimated from the seismic cloud. Furthermore, all three models agree that the increase in pore pressure induces the seismicity and that the seismic cloud represents the stimulated volume. The PoSeMo model considers fluid flow in the permeable fractures for the hydraulic model, whereas seismogenic index considers the pore pressure perturbation of the diffusion process and McGarr (2014) does not involve hydrologic modeling explicitly. This difference leads to different representations of the seismic cloud, i.e., a complex of grids or an ellipsoid.

  2. McGarr’s theory initially estimates all seismic moments and then it estimates the seismic moment of the largest event with the b-value. The PoSeMo model also estimates all the seismic moments first, but the maximum magnitude is estimated by comparison with observational values. Shapiro et al. (2011) estimate the maximum magnitude directly from a statistical method based on the Gutenberg-Richter relationship and use the axis of the ellipsoid to represent the seismic cloud.

  3. In McGarr’s theory, deformation of the rock mass because of the injection and the estimated maximum magnitude is connected with the b-value. Shapiro et al. (2011) estimate the progress of the stimulated volume using the diffusion model and then estimate the maximum magnitude with the b-value. The PoSeMo model obtains a sample of the index having magnitude information from one local region, which is then applied to the entire stimulated volume.

APPLICATION TO FIELD DATA

Microseismic data from Basel, Switzerland

We applied the concept outlined above to the microseismic data set of the Basel EGS project (Häring et al., 2008) to test the feasibility of our model. The Basel EGS project was conducted mainly by Geothermal Explorers Ltd. (GEL, now Geo Explorers Ltd.). In December 2006, hydraulic stimulation was implemented via injection well Basel-1 (true vertical depth: 5000 m) in Basel’s urban area. A total volume of 11,500 cubic meters of water was injected more than six days into the open-hole section (4603–5000 m), and the wellhead pressure rose to 30 MPa (Häring et al., 2008). Hydraulic stimulation activated a significant volume of enhanced artificial permeability within the geothermal reservoir causing several microseismic events. Because the predetermined level of event magnitude was exceeded under almost maximum wellhead pressure, the stimulation was halted according to the traffic light system (Häring et al., 2008). Just after the shut-in, which was performed shortly after stimulation stopped, several felt events occurred, including the largest event of ML 3.4, which ultimately led to the cancellation of the EGS project (Häring et al., 2008).

The microseismic monitoring network, which consisted of six borehole stations, was deployed by GEL. Members of our research group have worked on the Basel microseismicity project, as well as for GEL and the Swiss Seismological Service. Asanuma et al. (2008) determined the precise hypocenter locations of microseismicity using microseismic wave data provided by GEL with cluster analysis and the double difference (DD) method (Waldhauser and Ellsworth, 2000) as shown in Figure 2. The hypocenters of the microseismic events delineated the north–northwest/south–southeast subhorizontal reservoir showing consistency with the regional stress state of the strike slip regime (Valley and Evans, 2009). The largest event (Mw 3.51) was located in the bottom part of the seismic zone. Figure 3 shows the hydraulic records and occurrence of induced seismicity with magnitude information. We were able to find that the largest event and several other large events occurred at the shut-in phase. Our previous study has revealed several characteristics related to some selected large events (Mukuhira et al., 2013).

The Mw used here was an updated version of Mukuhira et al. (2013), originally provided by GEL and the Mw used in this study was estimated from the waveform amplitude, as in Dyer et al. (2008). On this scale, the largest event that occurred just after the shut-in was Mw 3.51, which is known as the ML 3.4 event on the local magnitude scale of the Swiss Seismological Service (Deichmann et al., 2014). The value of Mw is transformed into the seismic moment using  
Mw=23logMo6.07.
(4)

Parameter determination

Seismic moment density D

To estimate a robust value for D, it should be noted again that the grid size used to estimate D is crucial. When the grid size becomes bigger, it could be expected that D would be small because the denominator of equation 2 becomes bigger and thus, many grids might not include event hypocenters, except for a few grids containing most of events. Therefore, it is important to determine the appropriate grid size for each case study. We first conducted a random selection test to investigate the effects of grid size and choice of grid on the D calculation. In Figure 5a, ƒ shows the frequency distribution of D resulting from the random test for various grid sizes. For smaller grids (20, 40, and 60 m), there is a clear peak in histogram, but we were unable to determine a clear peak for larger grid sizes (approximately 100 m). Figure 5b shows the average of randomly estimated D values with two standard deviations (2σ) for various grid sizes. As Figure 5a and 5b shows, the deviation of D increases with grid size. In particular, Figure 5b shows that the deviation for grid sizes larger than 60 m is significantly large. This means that D can be affected dramatically by choosing a grid that is larger than 60 m. This feature can be understood as the balance between grid size and the existence of events within the chosen grid. Therefore, a grid size of 20 or 40 m would be considered preferable; however, in the case of the 20 m grid size, the D value is significantly higher compared with the 40 m grid size and the average value. This could also be attributable to the imbalance between grid size and events. Therefore, in this case, a grid size of 40 m should be appropriate for calculating D regardless of the concept used to choose the grid for the D estimation. We also investigated the temporal change in the D estimation for the 40 m grid. We searched for the grid with the maximum number of events occurring at each 2 h time step. Then, we estimated the D values shown in Figure 5c. It should be noted that because of the assumption of a point source for model simplicity, we did not consider the case in which the rupture area was bigger than the grid size.

Stimulated rock volume Vstim

An appropriate grid size for the Vstim estimation for the Basel case was also considered. As mentioned in the previous section, a grid size that can detect single events is desirable. First, we calculated the minimum distance to the closest event for all combinations of every single event from their hypocenter locations, which provided an average value of 11.95 m. Therefore, we considered that the grid size should be at least twice this value and a grid size of 25 m was deemed reasonable. The detectability of the variously sized grids was evaluated. We discretized the study area with various sizes of grid (10–50 m). Then, we counted the number of events in the grids with respect to all the seismic events that occurred within the study area during the in study period. The results are summarized as the histogram in Figure 6a. Smaller grids of 15–25 m successfully detected single seismic events, but considering the interevent distance, the 25 m grid size was deemed optimum for the Basel case. In addition, Häring et al. (2008) also use a 25 m grid for the estimation for the entire stimulated rock volume; therefore, we decided to use the 25 m grid for the Vstim estimation.

The Vstim estimated using the 25 m grid is shown in Figure 4b. The color of the grids represents the number of events in each grid. Those grids located adjacent to the feed point are shown in light colors, and those in the periphery are shown as blue. This suggests that event density was higher in the near-field of the injection well than in the far-field. This is quite reasonable because pore pressure near the injection well could rise sufficiently high to initiate shear slip not only on well-orientated fractures but also on other fractures. As the grid distribution of Figure 4b shows, this method can evaluate Vstim precisely and realistically. Our Vstim can take into account the shape of the reservoir, fracture system, and even the seismically silent zone, because this method uses all the microseismic information. The distribution of the grids corresponds to the primary model of the permeable fracture system within the reservoir, in which the permeability was enhanced by shear slip. However, several grids are shown spatially isolated from the main body of Vstim in Figure 4b. These independent grids suggest that this method requires improvement to model the flow system precisely, and this will be addressed in our future studies.

Because this method is based on microseismic hypocenters, we must take into account the spatial uncertainty of their locations. The procedure to estimate the stimulated volume can use the grids in place of the hypocenter locations. Therefore, when there is large uncertainty regarding the location of a single seismic event, even though this event may not be detected by the best grid, it will be detected by other grids as long as its hypocenter is within the study area. However, there is the possibility that this event might be detected in a grid that already has one seismic event. Consequently, spatial uncertainty does not affect the results of Vstim significantly because most of the 25 m grids detected single seismic events. Therefore, instead of spatial uncertainty regarding every single event, we evaluated the uncertainty in discretizing the study area. The center of the study area was the injection point, and that point was considered the center point of the central grid. We conducted a bootstrap test that randomly changed the center point of the central grid 300 times within the initial central grid as shown in Figure 6b. It can be seen that the uncertainty associated with discretizing the study area increases with time, but it remains sufficiently small to be neglected (Figure 6c).

RESULTS AND DISCUSSION

Evaluation of seismic moment density

We calculated D, Vstim, and Mopossible using a 40 m grid for D and a 25 m grid for Vstim to check whether our model could reproduce the cumulative seismic moment derived from the observational data. Using all the microseismic data from the study area, the computations of Vstim and D were implemented every 2 h (Figure 7a and 7b, respectively). The results for Vstim show a smooth increase during and after the stimulation. The gradient of Vstim changed several times at approximately 3.0, 3.8, 5.1, 5.8, and 6.2 days, and it became steeper with the time. Clear correlation between the magnitude history shown in Figure 7a and the hydraulic records shown in Figure 6c was not observed; however, the gradient of Vstim appeared to change several hours after the increase of the flow rate. The gradient of Vstim was steepest at day 6 of the stimulation when the wellhead pressure was 30 MPa. Then, at the shut-in, the Vstim gradient decreased almost simultaneously, which was the only case of synchronization between the gradient change and hydraulic operation. Subsequently, Vstim increased logarithmically with a smaller gradient. This is also an interesting result showing that our proposed method was able to detect a slight change in Vstim successfully; however, the investigation of Vstim is beyond the scope of this study.

For the D estimation, we prepared two different types of scenario. The first scenario was to choose the grid with the maximum number of events as the optimum grid in every time step (referred to as Max num). The value of DMax num decreased at around day 3.5 and it kept decreasing until around day 6.0. This observation can be interpreted as many smaller events being included in the grid chosen in that term. The other scenario for D was to choose the same grid that maximized the number of events on day 2 (hereafter, Dsecond day). Hence, DMax num and Dsecond day should have the same values on day 2. As day 2 was an early stage of the stimulation, the near-field of the injection was the area primarily stimulated. Therefore, the grid used to estimate Dsecond day should be located within that region. This grid might show the best response to the change in the pore pressure because the pressure in the well can directly reach the near-field. As Figure 7b shows, Dsecond day increased continually during the stimulation period because the same grid was chosen at each 2 h interval.

To evaluate the model output of DMax num and Dsecond day, we also calculated the observational seismic moment density Dobs, which is obtained from the summation of the observed seismic moment of all the seismic events at each time step divided by Vstim, as in equation 5:  
D(observed)=ΣMo/Vstim.
(5)

It was revealed that Dobs maintained a lower value from the start of the stimulation until the occurrence of the largest event, whereupon it increased drastically. Then, Dobs decreased slightly because of the increase of Vstim and the values nearly converged. Considering our assumption regarding D, the convergent Dobs can be considered as a representative parameter describing the average seismic moment within this study area. Therefore, our assumption regarding D was correctly proven. The convergent Dobs is the parameter that should be estimated as early as possible during the stimulation process. However, practically, we are able to obtain a representative value of D only after hydraulic stimulation and the occurrence of a series of seismic events. Therefore, it is necessary to estimate several values of D based on different scenarios. In this case, Dsecond day could be used to estimate the representative D at around day 4 of the stimulation. However, DMax num was underestimated throughout the stimulation period, and it did not reach the representative value of D.

The difference between Dobs and D based on some scenario can indicate the possible level of seismic hazard. If the difference was large, the probability of occurrence of larger events would increase because D represents the PoSeMo per rock volume, even though D cannot provide specific information on the seismic hazard. Therefore, this value should be used in the form of the PoSeMo model, after multiplication by the stimulated volume Vstim, because the seismic moment density itself is not an easily interpreted quantity and it does not constitute an informative index.

Evaluation of PoSeMo

The value of Mopossible based on DMax num and Dsecond day was calculated and compared with the observed value (ΣMoobs) in Figure 7c. Here, the value of ΣMoobs is the observed value that is independent of our model. Because Vstim increased smoothly and simply, Mopossible shows behavior similar to that of the D value. The functions MopossibleMax num and Mopossiblesecond day were overestimated significantly with respect to ΣMoobs, which is used for the evaluation of seismic activity. The value of MopossibleMax num was overestimated during the stimulation phase, but it showed a good fit overall, whereas the value of Mopossiblesecond day was always higher than MopossibleMax num and ΣMoobs. Just after the occurrence of the largest event, ΣMoobs increased drastically and it overtook both model outputs. After the occurrence of the largest event, Mopossiblesecond day and ΣMoobs converged to almost the same value, but MopossibleMax num was lower than ΣMo.

To compare the relevant model of McGarr (2014), we estimated the cumulative seismic moment based on his model using equation 6:  
ΣMo=2μ(3λ+2G)3ΔV,
(6)
where μ is the friction coefficient, λ and G are Lame’s elastic parameters (G is the modulus of rigidity), and ΔV is the volume of the injected fluid. Equation 6 can be simplified to  
ΣMo=2GΔV,
(7)
when it is assumed that μ=0.6 (a typical result from laboratory tests) and that λ=G. To use equation 6, assuming a friction coefficient of 0.85, λ=31.561  GPa, which is estimated from the velocity of the P-/S-waves (VP=5940  m/s and VS=3450  m/s), the density of granite (2.75g/cm3) from Mukuhira et al. (2013), and G=26.639  MPa from E=65  GPa and ν=0.22 for Poisson’s ratio (Valley and Evans, 2015), we calculated the cumulative seismic moment based on McGarr’s model. Here, McGarr’s cumulative seismic moment specified for each field using equation 6 is referred to as ΣMoMcGarr Basel, and for the general case based on equation 7, it is referred to as ΣMoMcGarr. The estimates are shown in Figure 7c. The functions ΣMoMcGarr Basel and ΣMoMcGarr increased as a function of the injected volume (maybe time). Then, they started to decrease with the flow back of the injected fluid. The value of ΣMoMcGarr Basel was larger than ΣMoMcGarr because λ was larger. The value of ΣMo estimated by McGarr’s model showed a significantly higher value than our model and twice that of ΣMoobserved. This can be interpreted as McGarr’s model not considering the missing injected fluid related to the highly permeable part of the reservoir, as well as to aseismic deformation, as mentioned by Kwiatek et al. (2015).

Evaluation of possible maximum magnitude

As previously described, our model evaluates the potential maximum magnitude based on the discrepancy between Mopossible and ΣMoobs, indicating the cumulative seismic moment that could be released in the future. First, the discrepancy between Mopossible and ΣMoobs is calculated, and this is then translated into the maximum Mw under the assumption that this differential seismic moment will be released in a single shear slip. It should also be noted that when ΣMoobs is larger than Mopossible, our model is unable to provide a magnitude.

Figure 7d shows the possible maximum magnitude estimated by our model. On day 2 of the stimulation, the values of MwmaxMax num and Mwmaxsecond day were greater than 2.0, which were successfully confirmed by the occurrence of an event of Mw>2.0 at around day 2.3. Except for days 4–5, MwmaxMax num was able to provide estimates that gradually increased with time. The value of MwmaxMax num increased significantly at around day 6, which was about one day before the occurrence of the largest event. Just before the occurrence of the largest event, MwmaxMax num showed a value almost identical to the Mw of the largest event. On the other hand, Mwmaxsecond day provided estimates after day 2, which also increased gradually with time. Furthermore, the value of Mwmaxsecond day was nearly the same as the Mw of the largest event before the occurrence of the largest event. After the occurrence of the largest event, neither Mwmaxsecond day nor MwmaxMax num could provide estimates, whereas ΣMoobs increased drastically. This occurred because of the underestimation of Vstim. In the estimation of Vstim, the effect of the source size is not considered; however, in reality, the stimulated volume of larger events would be bigger than for smaller events. Therefore, after the occurrence of the largest event, Vstim would be underestimated, and this caused the inability to estimate Mwmax.

After day 8.5, Mwmaxsecond day provided estimates once again, which suggested that an event of approximately Mw 2.5 could occur. This estimate could be consistent with the occurrence of other larger events after one month of stimulation (Mukuhira et al., 2013). However, considering the fluid flow hydraulic model behind our model, pore pressure increase would not be responsible for the occurrence of induced seismicity at that time. Therefore, it is difficult to consider whether the estimates from this study would be valid after one month of stimulation. The validity period of our model will be studied further in future work.

Overall, the output of our model (Mwmaxsecond day and MwmaxMax num) was able to estimate the maximum Mw successfully, showing good agreement (almost identical) with the maximum Mw before the occurrence of the largest event. Therefore, our model demonstrated reasonable performance in the evaluation of the seismic activity and the estimation of the maximum magnitude for the Basel case.

We compared the results of our estimates of maximum Mw with those from other related studies. By introducing the b-value from the Gutenberg-Richter relationship to equation 6, McGarr (2014) derives  
Mo(max)=1BB2μ(3λ+2G)3ΔV,
(8)
where B=b/1.5. To obtain a simplified and generalized description, when b=1 (i.e., B=2/3), equation 7 becomes  
Mo(max)=GΔV
(9)
under the assumptions for equation 7. Here, we estimated the maximum Mw based on equation 8, which is referred to as MwmaxMcGarr Basel, with b=1.39, which we estimated from our seismic catalog. Although the b-value should be estimated at every time step, we used a constant value for simplicity because this was not a principal objective of this study. Equation 9 can also be translated to MwmaxMcGarr as well.
Kwiatek et al. (2015) extend McGarr’s theory and proposed equation 10:  
Mo(max)=1BB2μVΔP¯,
(10)
where V is the formation weakened by the injection of fluids and ΔP is the average increase in pore pressure. Kwiatek et al. (2015) use V and ΔP instead of elastic parameters and the volume of the stimulated fluid. This theory can introduce microseismic information because V is estimated from the volume of the ellipsoid of the modeled cloud. We conducted principal component analysis to model the seismic cloud with the ellipsoid at every time step, and we defined the length of the axes as 2σ of each direction of the eigenvector, as in Kwiatek et al. (2015). Ideally, we should estimate the pore pressure increase based on the diffusion equation (Shapiro and Dinske, 2009) because the ellipsoid assumption appears based on diffusion pressure propagation. Here, we used 11.01 MPa for ΔP, which is the average value of our estimates of pore pressure increase based on a geomechanical method (see Mukuhira et al., 2013). This value is considered reasonable because microseismic events started to be observed soon after the wellhead pressure reached 10 MPa (see Figure 3 or 5). With these parameters, we estimated the maximum Mo(max) that was then transformed into MwmaxKwiatek.
Shapiro et al. (2011) also propose the following equation to estimate the maximum magnitude with the Y value related to the minimum principal axis of the fluid:  
Mw(max)=2log10Y+(log10Δσlog10C)/1.56.03,
(11)
where Δσ is the average static stress drop and (geometric constant) C=1. We can assume Y is the length of the minimum principal axis of the ellipsoid that models the stimulated volume. For the input stress drop for equation 11, we previously estimated the average value as Δσ=0.026  MPa (Mukuhira et al., 2013). However, this stress drop appears slightly smaller than those for induced seismicity from other fields (Kwiatek et al., 2015), and therefore, we decided to use Δσ=3  MPa as an acceptable common value, as in Deichmann et al. (2014). The maximum magnitude estimated from equation 11 is referred to as MwmaxShapiro.

Figure 7d shows a comparison of the values of Mwmax estimated in related studies with our estimates. The functions MwmaxMcGarr and MwmaxMcGarr Basel increased logarithmically with the increase of injection volume (time), and they switched positions after the conversion to Mw because of the term that includes the b-value. Both of McGarr’s values of Mw showed very good fit with the maximum Mw=3.51, especially MwmaxMcGarr. However, considering the overestimation of ΣMoMcGarr and ΣMoMcGarr Basel, it is questionable why MwmaxMcGarr and MwmaxMcGarr Basel produced comparable values. Interestingly, the same features were observed for MwmaxKwiatek. We traced back to equation 10 and estimated the cumulative seismic moment based on the ellipsoidal volume of the formation weakened by the injection with the relation ΣMo=2μVΔP, which is also shown as ΣMoKwiatek in Figure 7c. The curve of ΣMoKwiatek showed a different shape to that of ΣMoMcGarr but a similar shape to our Vstim, because it introduced the volume of ellipsoid modeling of the seismic cloud. However, the seismic moment remained significantly overestimated. Although MwmaxKwiatek showed almost the same value as MwmaxMcGarr Basel and it provided a reasonable estimate to some extent, MwmaxShapiro showed a lower value than the other models and it underestimated the magnitude throughout most of the study period.

To summarize, the models based on McGarr’s model, i.e., MwmaxMcGarr, MwmaxMcGarr Basel, and MwmaxKwiatek, showed acceptable estimates of the maximum Mw, but considering the overestimations of ΣMoMcGarr, ΣMoMcGarr Basel, and ΣMoKwiatek, the reliability of their estimates of Mwmax remains open to question. The MwmaxShapiro model did not perform well in this case. The reason for this could be that the shape of the seismic cloud for the Basel case was not ellipsoidal. Although this model did perform reasonably well for the ellipsoidal seismic cloud of Soultz (Shapiro et al., 2011), as in the Basel case, the reservoir had a branched structure (Häring et al., 2008) that made it difficult to represent the seismic cloud with an ellipsoid. It also remains unclear whether the density of seismic events was sparse at the periphery of the ellipsoid; i.e., it is not possible to interpret the part of the ellipsoid where no seismic events occurred.

Comparison with relevant studies revealed that our model performed reasonably well for the Basel case in the estimation of Mopossible and maximum Mw. Some relevant models that evaluated the ellipsoid volume as the stimulated volume had difficulty in evaluating the seismic moment or maximum magnitude; however, this does not mean that our model always performs better than the models of related studies. Each of these models was based on different concepts and each performed well in their study, just as our model performed well in the case of microseismic data from Basel. Therefore, we should apply our model to different fields to investigate its feasibility.

As shown here, the application of the our model to real microseismic data has demonstrated the feasibility of using this model to evaluate seismic activity and to estimate the possible maximum magnitude in near real time using only microseismic information, especially in the short and middle terms. Therefore, our physics-based model of seismic activity could be used for comparison with other models based on statistical concepts.

CONCLUSIONS

We presented our original PoSeMo evaluation model and introduced the underlying theory. In the PoSeMo model, the PoSeMo is defined as the seismic moment that could be released in the future, which is expressed as the product of the seismic moment density and the stimulated rock volume due to hydraulic stimulation. The PoSeMo is evaluated from microseismic information from past to present and compared with the cumulative seismic moment, i.e., the observed value. Then, the difference from the observed value can be translated into the potential future maximum magnitude that could be used for seismic risk assessment. We applied the model to a microseismic data set from Basel, Switzerland and showed that it was capable of evaluating the PoSeMo and the maximum magnitude before the occurrence of the largest event. The feasibility of our model was demonstrated on a practical level by comparison with other relevant studies, when appropriate values of D and Vstim were available. It was also shown that our model could provide useful physics-based information regarding the hydraulic stimulation operation.

We intend to apply this model to other cases of fluid-induced seismicity from EGS or to shale fields, as well as to test its wider feasibility and universal features. Moreover, the future development of the PoSeMo model will be more physics orientated. In this work, we used only the seismic moment and hypocenter distribution, but recent progress in microseismic monitoring might be able to provide additional physical information, such as focal mechanisms, fault length, and stress state. The other important aspect of this model is how local information can be used to describe the randomness of an entire physical feature. The results present the interesting suggestion that one representative parameter sampled from one part of the reservoir could be used to describe the general physical features. This should be clarified through the application of this model to various other fields.

ACKNOWLEDGMENTS

We thank Geo Explorers Ltd. for providing the microseismic data sets and for permission to publish the results. We are also grateful to A. McGarr, two anonymous reviewers, and editor for their constructive comments. We also thank J. D. Zechar, E. Kiraly, V. Gischig, H. Moriya, and C. Dinske for discussions and valuable feedback on this research. This study was supported by JSPS KAKENHI (grant number 15K20865).