## 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.

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

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.

#### Evaluation of seismic activity

To evaluate seismic activity, $Mopossible$ is compared with the cumulative seismic moment, which is referred to as $\Sigma Mo$ (observed), calculated from all the seismic events that occurred by the given time. When $Mopossible$ is larger than $\Sigma 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 $\Sigma 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.

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.

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.

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).

### 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, $\u0192$ 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\sigma $) 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.

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 ($\Sigma Moobs$) in Figure 7c. Here, the value of $\Sigma 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 $\Sigma 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 $\Sigma Moobs$. Just after the occurrence of the largest event, $\Sigma Moobs$ increased drastically and it overtook both model outputs. After the occurrence of the largest event, $Mopossiblesecond day$ and $\Sigma Moobs$ converged to almost the same value, but $MopossibleMax num$ was lower than $\Sigma Mo$.

### Evaluation of possible maximum magnitude

As previously described, our model evaluates the potential maximum magnitude based on the discrepancy between $Mopossible$ and $\Sigma Moobs$, indicating the cumulative seismic moment that could be released in the future. First, the discrepancy between $Mopossible$ and $\Sigma 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 $\Sigma 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 $\Sigma 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.

_{o}(max) that was then transformed into $MwmaxKwiatek$.

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 $\Sigma MoMcGarr$ and $\Sigma 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 $\Sigma Mo=2\mu V\Delta P$, which is also shown as $\Sigma MoKwiatek$ in Figure 7c. The curve of $\Sigma MoKwiatek$ showed a different shape to that of $\Sigma 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 $\Sigma MoMcGarr$, $\Sigma MoMcGarr Basel$, and $\Sigma 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).