Magnitude and nucleation time of the 2017 Pohang Earthquake point to its predictable artificial triggering

A damaging Mw5.5 earthquake occurred at Pohang, South Korea, in 2017, after stimulating an enhanced geothermal system by borehole fluid injections. The earthquake was likely triggered by these operations. Current approaches for predicting maximum induced earthquake magnitude (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{\max }$$\end{document}Mmax) consider the volume of the injected fluid as the main controlling factor. However, these approaches are unsuccessful in predicting earthquakes, such as the Pohang one. Here we analyse the case histories of induced earthquakes, and find that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{\max }$$\end{document}Mmax scales with the logarithm of the elapsed time from the beginning of the fluid injection to the earthquake occurrence. This is also the case for the Pohang Earthquake. Its significant probability was predictable. These results validate an alternative to predicting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{\max }$$\end{document}Mmax. It is to monitor the exceedance probability of an assumed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{\max }$$\end{document}Mmax in real time by monitoring the seismogenic index, a quantity that characterizes the intensity of the fluid-induced seismicity per unit injected volume.

In November 2017 near the city of Pohang, South Korea, shortly after borehole fluid-injection operations targeted to a stimulation of an EGS, an unexpectedly strong earthquake occurred. The earthquake was likely triggered by the EGSstimulation operations 2,[9][10][11][18][19][20][21] , which raised questions regarding the ability to forecast the maximum magnitude of an anthropogenic earthquake.
There is considerable active research [22][23][24][25][26][27] in developing strategies to forecast M max . Popular approaches consider the injected fluid volume, ΔV f , as the controlling parameter 22,24,25 . Some approaches 22,28,29 implicitly assume that fluid injection is the sole source of strain in the surrounding rocks, and conclude that the seismic moment of the largest induced earthquake (along with the cumulative seismic moment of the induced seismicity), M max 0 , is limited by a quantity that is proportional to the fluid volume, whereby M max 0 / ΔV f . Especially in seismically active regions, this singular role of fluid injection on the local strain may be problematic, because this simple assumption fails to adequately account for tectonic sources of strain and effective stress changes 27 .
Another approach is to compute the moment of the largest arrested rupture 25 that was induced by pore-fluid pressure perturbation in the layer around an area of the earthquake fault plane. A comparison of the stress intensity factor (a function of the pore-pressure perturbation) and fracture toughness yields a different power law 25 , M max 0 / ΔV 3=2 f . In this power law, the proportionality coefficient 25 depends on the background stress drop, friction, elastic properties and thickness of the pressureperturbed layer, with the seismic moment of the maximum arrested rupture eventually being controlled by the injectioninduced pressure perturbation. Such a perturbation is independent of the existing tectonic stresses at the injection site. Thus, such an approach has a potential for failing in tectonically active regions and areas situated near critically stressed faults. However, one can adjust the proportionality coefficient of the power law M max 0 / ΔV 3=2 f to the tectonic conditions 25 by relating it to the seismogenic index, Σ, which quantifies the induced seismicity that may be produced at a given injection site in response to a unit injected volume 30,31 (see 'Methods' section).
A direct application of the seismogenic index model 23,30,31 to the largest earthquake in a time series of fluid-induced seismicity has also provided an estimate 24 of M max . However, such an estimate is derived from a statistical formulation of the magnitude-frequency distribution and yields the magnitude of an event whose occurrence probability is nearly 63% (see 'Methods' section). This approach provides no constraints on the potential to trigger larger earthquakes.
Recent studies on the 2017 Pohang Earthquake have emphasised that none of the abovementioned approaches are adequate in describing this earthquake 2,9,10,19 , whereas case histories 22,25 have shown that M max is significantly correlated with ΔV f . However, magnitude scaling with ΔV f is a function of the geometry of the processes under consideration, as three 22 -and two 25 -dimensional pore-pressure-diffusion-based approaches have yielded various power-law dependencies of M max 0 ðΔV f Þ. The duration of the earthquake triggering can be more indicative of the underlying physical process. For example, the time scaling of normal diffusion is the same in both of the abovementioned cases. Moreover, ΔV f is closely related to the elapsed time ΔT from the beginning of fluid injection to the M max earthquake.
Here, we analyse published case histories and show that M max scales with log ðΔTÞ. Furthermore, this magnitude scaling approach makes the artificial triggering nature of the 2017 Pohang Earthquake more noticeable than magnitude scaling with ΔV f . However, ΔT is unknown a priori; therefore, the seismogenic response of the surrounding rocks due to fluid injection should be monitored in real-time to address the M max problem. We propose to monitor the worst-case probability of a hypothetical M max earthquake and show that such a strategy could be useful in the cases of Pohang and other damaging induced earthquakes (e.g. the Denver seismicity of 1965-67). We show that a significant probability of the 2017 Pohang Earthquake was predictable.

Results and discussion
Maximum magnitudes scale with logarithm of elapsed times. We used existing compilations 22,24,25,28,29 of published M max data for fluid-injection-induced earthquakes, and included ΔT data from corresponding and some additional literature sources [5][6][7]12,18,24, . In some publications times, ΔT can be directly found. However, in the other cases, we estimated ΔT from published plots. We provide ΔT values up to the second digit in the leading order. This precision is sufficient for our consideration below. The data and corresponding references are given in the Supplementary Data 1 file.
We consider M max as functions of the injected fluid volumes, ΔV f , and the elapsed times, ΔT, (Fig. 1). The complete data set includes field observations of fluid injections in deep boreholes (data points with Mw > −3), and observations obtained from mine-like (underground laboratories) and laboratory experiments (data points with Mw < −3). The laboratory and mine-like observations are characterised by the presence of free-surfacetype boundaries in the vicinity of the injection sources: sample surfaces or mining tunnel walls. The results from these experiments must not be described by the same scalings and/or factors as the field observations in deep boreholes since the latter are better approximated by an infinite continuum. We, therefore, consider the field observations of fluid injections in deep boreholes separately to the mine-like and laboratory observations (Fig. 1d, e). We find that M max is well correlated with both ΔV f and ΔT (Fig. 1). However, magnitude scaling with ΔT describes M max better by providing a more compact set of data points (Fig. 1d, e). For example, a linear regression of the ΔT -M max scaling provides significantly better statistics than the equivalent fitting of the ΔV f À M max scaling (e.g. a significantly higher coefficient of determination and significantly lower standard deviations of fitting parameters (Supplementary Information file, Eqs. 10 and 11)). Moreover, ΔT -M max scaling is also significantly more successful in describing the 2017 Pohang Earthquake.
Our findings suggest the following. The growth of the characteristic size, L, of an impacted rock domain (stimulated volume) due to fluid injections can be governed by various processes, such as (non-)linear pore-pressure relaxation 23 , which sometimes includes a significant portion of the poroelastic coupling 59 , aseismic deformation 60,61 , and/or a diffusion-driven event interaction 11 . The growth of L can be approximated as a power-law function, ΔT ξ , where ξ is indicative of the physical process(es) controlling this growth. For example, ξ = 1 in the case of a ballistic propagating perturbation, ξ = 1/2 in the case of a linear pore-pressure relaxation (a normal pore-pressure diffusion), and 1/ 3 ≤ ξ < 1/2 in the case of a non-linear diffusion-like process (an anomalous diffusion process) of opening additional 3-D pore space 23 . Asymptotically, such power laws apply to both the triggering front and backfront 62 (approximate outer and inner envelopes of the seismically activated stimulated volume). L also continues to grow after the termination of injection 62 . The size L influences the frequency-magnitude statistics of the induced seismicity 23 , because the length statistics of triggered rupture surfaces that are related to the stimulated volume (i.e. contacting, intersecting or belonging to it) are controlled by this size 23 . The probability of a triggered earthquake with a large rupture length, R max , which corresponds to a given M max , increases with increasing L. This is also the case for so-called runaway ruptures because increasing the stimulated volume raises the probability that it will contact a large fault 23,25 . This relationship can be approximated as R max / L because the contacting probability 23 is a function of R max =L. The well-known relationships between the rupture length R of an earthquake and its seismic moment, M 0 ∝ R 3 , and between the seismic moment and the moment magnitude 22,23,24,25 , Mw¼ 2ðlog M 0 À 9:1Þ=3, then yield M max 0 / L 3 and M max % 2log 10 L þ const. The substitution of L ∝ ΔT ξ into these relationships yields the following scalings: M max % 2ξlog 10 ΔT þ const and M max 0 / ΔT 3ξ . We find that the scalings M max 0 / ΔT 3=2 and, respectively, M max % log 10 ΔT þ const generally describe the field-scale observations (Fig. 1b, e and the linear regression results, Supplementary Information file, Eqs. 10 and 11). These relations correspond to the power law L ∝ ΔT 1/2 , which is typical for the size-time dynamics of rock volumes stimulated by normal pressure diffusion processes. Thus, our observations indicate a statistically dominant role of linear pore-pressure diffusion in triggering M max earthquakes. The following three scenarios 23 can be relevant for such earthquakes: (i) the rupture surface of the earthquake is contained within the stimulated volume (sometimes strictly such events are called induced earthquakes 18,19,23 ); or (ii) the nucleation domain of the rupture surface and its part are contained within the stimulated volume; or finally, (iii) the rupture surface is touched or slightly intersected by the stimulated volume (sometimes events, addressed in the scenarios (ii) and (iii) are called triggered ones 18,19,23 ). The cases of ruptures arrested inside or outside of the stimulated volume as well as runaway ruptures are parts of these scenarios. At least in the scenarios (ii) and (iii), magnitudes of triggered earthquakes are very likely controlled by tectonic features of corresponding geological sites rather than by the scale L of the stimulated volume. However, the occurrence probability of such earthquakes depends on this scale, as it follows from our consideration above. The Pohang Earthquake is an example of such an event. Figure 2a-e indicates that it corresponds to the scenario ii. This is also in agreement with the conclusion of the ORAC Committee 18 that the Pohang Earthquake was triggered by an EGS-stimulation impact on the earthquake hypocenter domain. This impact could have various forms, for example, a direct injection-produced pore-pressure perturbation or a combined pore-pressure and stress perturbations 11 of induced earthquakes in the stimulated volume.
The abovementioned magnitude scalings naturally include situations of earthquake triggering after the termination of fluid injection since the law, L ∝ ΔT ξ , also describes the continuing growth of the stimulated rock volume after the termination of injection 62 . These magnitude scalings are therefore more adequate than those with ΔV f . Large earthquakes have an enhanced occurrence probability at or shortly after the In respect to the post-injection growth of the size L of the stimulated volume, we should note that this volume will be relevant for earthquake triggering if corresponding perturbations of the pore-pressure and/or of the (poro)elastic stresses are still significant (e.g. above their fluctuations of tidal and seasonal nature). Thus, the scaling M max À ΔT has its natural limits. For example, in the case of a pressure-diffusion earthquake triggering, the seismogenic post-injection time period is of the order of the total duration time of the preceding injection operations, t 0 , (see 'Methods' section). Then, the corresponding realistic scale of the spatial domain where the earthquake triggering will be probable is limited by the size of the order of several characteristic lengths ffiffiffiffiffiffiffiffiffi ffi 2t 0 D p , where D is a representative hydraulic diffusivity of hydraulic paths to the critically stressed faults.
Monitoring the probability of a hypothetical M max . Maximum induced magnitude and ΔT are unknown a priori parameters for a particular fluid-injection experiment. However, one could estimate the exceedance probability of an assumed critical M max during a given injection. The magnitudes of large runawayinduced ruptures are determined via the surrounding tectonic fault networks [23][24][25]64 , with their frequency-magnitude distribution given by the Gutenberg-Richter statistic of the tectonic environment. However, their probability is enhanced due to the fluid injection. We take this into account by combining the occurrence probability of an earthquake with a magnitude ≥ M max with the seismogenic index model 23,30 (see 'Methods' section). This model contains the seismicity parameters, Σ and b, which are functions of the seismotectonic features of the surrounding fault systems. A more or less seismically active fault system will potentially be activated if these parameters exhibit non-stationary behaviour during rock stimulation. Therefore, Σ and b should be monitored in real-time. We propose to monitor the worst-case probability based on the real-time evolution of these parameters: where we introduce two functions of the observation time t: sup ΣðtÞ ¼ maxfΣðτÞj0 < τ < tg and inf bðtÞ ¼ minfbðτÞj0 < τ < tg. potential alternative is to obtain a realistic estimate of the b-value from (a priori) available regional seismicity data. However, real-time monitoring of Σ is feasible, such that the observation of significant events then becomes sufficient (see 'Methods' section and Supplementary Information file).
2017 Pohang Earthquake and other case histories. The epicentre of the 2017 Pohang Earthquake 18 was 510 m from the site of the Pohang-EGS project and~15 km away to the northnortheast from the centre of Pohang, a city with a population of 500,000. The earthquake occurred on 15 November 2017, after a series of water injections into two deep boreholes, PX-1 and PX-2, that were drilled into granodioritic rocks to 4215 and 4340 m depth, respectively (Fig. 2a-e). The hypocenters of the aftershocks of the Pohang Earthquake 65 are distributed approximately in the depth range between 2 and 7 km and in a lateral SSW-NNE zone of 11km long (Fig. 2b, c). The injection operations stimulated a one order of magnitude smaller rock volume (approximately of 1 km size) indicated by the induced seismicity 20 occurred before the Pohang Earthquake (Fig. 2d, e). Apparently, the hypocenter of the Pohang Earthquake was within the stimulated volume (corresponding to scenario ii of earthquake triggering described above).
We obtain the approximate lower and upper bounds of Σ by treating the Pohang-EGS fluid injections in two different ways (see 'Methods' section). Σ is approximately between −2 and −1 (Fig. 3). We observe the tendency for Σ to increase with t. This tendency may be indicative of a gradual involvement of more seismically active domains in the stimulated volume. For example, there was probably an expansion of the stimulated volume to a more seismogenic part of the major fault system that was initially intersected by borehole PX-2 and indicated by a massive mud loss during the drilling operations 18 in late 2015. Moreover, typical fault rocks were contained in the cuttings from borehole PX-2 next to the mud-loss depth interval 18,19 . Thus very likely, a fluid injection occurred nearly directly into a pre-existing large-scale fault.
The Mw3.3 event on 15 April 2017 was particularly alarming, as it indicated a possible sudden increase in Σ to −1. Our estimates of Σ are completely based on the seismicity data coming from the stimulated zone. Thus they are not representative for the total focal area of the Pohang Earthquake. Neither are they indicative for the final size of the Pohang Earthquake. The reasons of rupture stopping or rupture arrest can be illuminated by investigations of the rupture segmentation during the Pohang Earthquake sequence 66 (early aftershocks). It has been observed that the initial propagation of the main rupture segment and its subsidiary segment was likely arrested by two other fault segments, one to the northeast and another one to the southwest from the hypocenter domain. However, our estimates of Σ indicate a dangerous tendency of non-stationarity in the temporal evolution of the seismogenic index, and thus an enhanced probability of triggering a runaway large-scale earthquake rupture. We use Eq. (1) to estimate such a probability of an assumed 5.5 magnitude earthquake (i.e. the probability of the Pohang Earthquake in this particular case history), which we denote below as M max . Please note that our estimates are not a prediction of a maximum possible earthquake magnitude but rather an estimate of the worst-case exceedance probability of an earthquake with the magnitude assumed.   Fig. 3d, green line). The worst-case probability of the 2017 Pohang Earthquake was ultimately very significant, 15-17% ( Fig. 3b-d, red and dark-red lines).
This worst-case probability was calculated using Σ, which characterises the stimulated zone, which is an order of magnitude smaller than the aftershock domain of the Pohang Earthquake. In the case of heterogeneous and/or unsteady tectonic stresses (the possibility of such a situation is indicated above by a nonstationarity of Σ; for example, the Mw5.4 Gyeongju earthquake 2,10,19 occurred on 12 September 2016 about 40 km south of the Pohang EGS could contribute to such a heterogeneous and/or unsteady stresses situation), the real probability of the Pohang Earthquake could be even higher. A strong indication of a possibility of a runaway rupture can also be seen from a very high worst-case probability of an event of the stimulation-zone size (this is a largest event for which Eq. (1) can be rigorously applied based on seismicity parameters observed in the stimulation zone). For the Pohang EGS, such an event has a magnitude of the order of Mw4-Mw4.5 (the rupture size of the order 10 3 m and the stress drop in the range 1-5 MPa). The temporal behaviour of the probability of a stimulation-zone size event during the Pohang-EGS-stimulation activity is similar to that calculated for the Pohang Earthquake (Fig. 3e). However, the probability of such an event has reached 50-80% (Fig. 3c, d the dark-red line). Such a high probability of a stimulation-zone size event was already achieved immediately after the induced Mw3.3 event (Fig. 3e). However, there were no earthquakes of the strength between Mw3.3 and Mw5.5. This is an indication that a similar event may have become the 2017 Pohang Earthquake with an unstable runaway rupture. Thus, during the stimulation operations, the Mw2.3 event should be considered critical, and a safe stimulation strategy would be to keep the induced seismicity approximately below Mw2.0.
The low b-value and high Σ increased the probability of the 2017 Pohang Earthquake. A comparison of the 2006 Basel EGS stimulation 12 , which injected more fluid (~11,500 m 3 ) than the Pohang one, yielded a significantly higher corresponding Σ (~0.25). However, the b-value (~1.44) was also higher, such that the probability of a Mw5.5 event at Basel was highly unlikely. However, the probability of a Mw3.4 event (such an earthquake led to the termination of the Basel EGS project 12 ) was~15% (Fig. 3c, green line). This level of probability could be predicted during the earlier stages of the Basel EGS project because an extremely high Σ could be observed shortly after the onset of the stimulation (a very high Σ remained nearly stationary throughout the injection period 31 ).
The 1962-1968 Denver earthquakes 4,43,44,68 are another example of how the induced seismicity could be predicted, with the large ΔV f being the main driver in this instance. Approximately 630,000 m 3 of wastewater was injected down a single 3671-m-deep well that was drilled into the fractured Precambrian crystalline basement below the Rocky Mountain Arsenal 44 (Fig. 4a). We obtain an estimate of Σ(t) for the Denver site by assuming a continuous injection experiment. This provides sup ΣðtÞ in the -2.5 to -1.75 range (Fig. 4b). However, we obtain a sudden increase in the seismogenic index to -0.7 when we assume that the second injection period that began in September 1964 was an independent one (Fig. 4b). We compute the probability of a hypothetical Mw5.5 event (Fig. 4c). We use this magnitude for a direct comparison with the Pohang case study; corresponding to various sources 4,22,43,68 none of the Denver earthquakes were larger than Mw5.3 or even Mw4.85. We assume b = 0.85, based on the published data 4 , and observe that even for the lower estimate of the seismogenic index, the probability of such a high magnitude event (>20%) was higher than that in Pohang (Fig. 4c). Furthermore, the probability was already high (~5%) before the beginning of the second phase of pressure-driven injections (Fig. 4c), and significant earthquakes (Mw4.6-4.7) had already begun to occur in early 1965. The largest earthquake 43 (≤Mw5.3) occurred on 9 August 1967.
The largest earthquakes in these three case histories occurred after the injection periods. The most critical parameters were a very high Σ (with a moderate ΔV f and high b-value) in Basel, a very low b-value (with a moderate to high Σ and very moderate ΔV f ) in Pohang, and an extremely high ΔV f (with a moderate to high Σ and nearly normal b-value) in Denver. The interplay between these parameters takes the form of Eq. (1). This equation forecasts a worst-case probability of an assumed maximum- magnitude earthquake. It can be applied for a real-time monitoring of underground fluid-injection operations. This equation (along with our Eqs. (8) and (9)) can become a useful ingredient of traffic-light-type 3 or value-at-induced-risk-type 21 management approaches of injection-induced seismicity (for example, contributing to them by monitoring the worst-case probability of a critically large stimulation-zone size event as a proxy of the runaway rupture probability). The elapsed injection time ΔT is also an important implicit parameter of Eq. (1). The occurrence time of the Pohang Earthquake and its magnitude are in agreement with our proposed ΔT -M max scaling for fluidinjection triggered seismicity (Fig. 1b, e). The nearly direct injection into the earthquake fault was likely a reason that a rather small injected fluid volume triggered a very large earthquake, in the upper right-hand domain of this scaling.
Presumably and in agreement with this scaling, the Pohang Earthquake nucleation required a rather long time of porepressure-related perturbation propagation along the fault zone. This perturbation propagation was forced by multiple borehole fluid injections into the underground.

Methods
Seismogenic index and event probability. We use the seismogenic index model 23,30,31 to analyse the case histories. Σ quantifies the potential for an earthquake to occur at a particular geologic site due to fluid injection. This model has been used to characterise induced seismicity in various studies, such as waterdisposal-related seismicity 8 , seismicity triggered by hydraulic fracturing operations 16,69,70 , and seismicity related to the Pohang-EGS activity 20 . The seismogenic index model 23,30 follows from the mass conservation of an approximately incompressible pore fluid (i.e. approximately pressure independent density). In the approximation of a monotonic non-decreasing injection pressure, this model describes the cumulative number of induced earthquakes, N ≥M (t), with moment magnitudes ≥M, that occur during the injection period t as: The model assumes the Gutenberg-Richter statistic for the frequency of magnitude ≥M induced events, log 10 N ≥ M ¼ a À bM, where a and b are the parameters of the Gutenberg-Richter distribution. Σ is the seismogenic index of the rock volume stimulated by the injection. Conventionally, V f (t) is given in m 3 . Σ can be estimated using Eq. (2) and a frequency-magnitude statistic of induced seismicity: Σ ¼ log 10 N ≥ M ðtÞ À log 10 ΔV f ðtÞ þ bM: We use this equation to compute Σ in our analysis (Figs. 3 and 4b) (details of these computations are shown in the Supplementary Information file and the Supplementary Data 2 file). We calculate two estimates of Σ in the Pohang case study to obtain approximate lower and upper bounds of this quantity as follows.
We first neglect the non-monotonic character of the injection operations and estimate Σ(t) directly by assuming that all of the injection phases are a single continuous injection experiment. This results in a somewhat underestimated Σ (an approximate lower bound) because the fluid volume was also accounted for during the intervals where the pore pressure was lower than those reached in the previous injection phases. We then assume that all of the injection phases are independent experiments and compute the corresponding values of Σ(t) ( Supplementary  Information file). This calculation results in a somewhat overestimated Σ (an approximate upper bound) because some events that were induced by previous injection phases can influence the number of induced events in the later injection phases. During long-term bleed-off phases of the injected fluid, the relation between these two Σ estimates can turn round because of possible late triggered events and resulting lower injected volumes at occurrence times of such events (Fig. 3a). Equation (2) has been proposed 24 to compute an estimate of M max . Substituting N ≥M = 1 (this must be valid for the largest event in the observed earthquake series) into Eq. (2) yields 24 M max ¼ 1 b ½Σ þ log 10 ΔV f ðtÞ. The b-value is usually close to 1, which means this relationship is very close to the result for the maximum arrested rupture 25 . Furthermore, this approach takes into account both the injected volume and seismogenic index that characterises the seismotectonic features of the injection site. However, the problem due to such an approach is inherited from the statistical nature of Eq. (2). If the observed events follow a Gutenberg-Richter magnitude distribution and the seismicity process can be approximated as a Poisson process (both are common realistic approximations), then the probability of the occurrence of a magnitude M or larger event will be 30 : We obtain W ev ≥ M max ¼ 1 À e À1 % 0:632 for the maximum expected event with N ≥M = 1. This M max estimate, therefore, predicts a magnitude of an event, which occurrence probability is 63% or even higher. It provides no constraints for larger magnitudes. If we use this estimate and express the term ½Σ þ log 10 ΔV f ðtÞ from Eq. (2), we will obtain an equivalent most-probable type 24 of M max estimate: where M is the magnitude used in Eq. (3) for calculating Σ. This estimate of the expected M max is shown in Figs. 3a and 4b by the yellow and dark-blue lines, respectively.
The probability of an arbitrary magnitude event can be estimated by combining Eq. (4) with Eq. (2) to obtain 30 : This equation provides estimates of the probability of triggered events with magnitudes ≥M occurring in the time period from the beginning of fluid injection at the time moment 0 until the given time moment t.
Event probability after injection termination. A specific model of the seismicity triggering process is required to predict the seismicity rate after the termination of a fluid injection. One can calculate the seismicity rate after the termination of injection 63 by assuming both a triggering mechanism that is governed by porepressure diffusion and a constant injection rate into a homogeneous porous medium. The decay rate of the induced seismicity is similar to Omori's law, which describes the decay rate of aftershock activity after tectonic earthquakes, whereby the decay exponent, p d , is larger than one 63 . The following approximation 63 yields the seismicity rate, R ev (t), after the termination of an injection at a constant rate (at least for τ 0 = O(1)): where t 0 is the duration of the fluid injection and τ 0 = t/t 0 denotes the normalised time after the termination of the injection (τ 0 ≥ 1).
The following approximation of R ev (t 0 ) for events with magnitudes ≥M can be directly obtained from Eq. (2): Numerical computations 63 have shown that p d approaches 2 for long durations (τ 0 on the order of 2 and larger). This is also supported by studies of the seismicity induced by large-scale massive fluid disposals in Oklahoma 8 . The estimated seismicity rates for the Fenton Hill and Soultz case studies 63 (τ 0 on the order of 1) have provided high values of 7.5 and 9.5, respectively, for the p d exponent.
We obtain the following relationship when we use Eq. (6) for the events with magnitudes ≥M that were triggered during the time range [t 0 , t] (fluid injection starts at the time moment 0): log 10 N ≥ M ðt 0 ; tÞ % ½Σ þ log 10 V f ðt 0 Þ À bM þ log 10 1 À ðt 0 =tÞ p d À1 p d À 1 : Equation (8) simplifies to log 10 N ≥ M ðt 0 ; tÞ % ½Σ þ log 10 V f ðt 0 Þ À bM À log 10 ðp d À 1Þ in the limit of long time periods after the termination of fluid injection. We obtain log 10 N ≥ M ðt 0 ; tÞ ≤ ½Σ þ log 10 V f ðt 0 Þ À bM for p d ≥ 2. Therefore, the event probability after the termination of fluid injection is roughly limited by its probability during the period [t 0 , t]. The short elapsed time of Δt after an injection termination causes the event number to increase as log 10 N ≥ M ðt 0 ; tÞ % ½Σ þ log 10 V f ðt 0 Þ À Mb þ log 10 ðΔt=t 0 Þ. Let us assume that the injection was terminated immediately after a magnitude M cr event (assumed to be critical) that occurred at t 0 . We can then use right-hand side of Eq. (3) with N ≥M = 1 to compute the seismogenic index, which will provide the following estimate of the event number: log 10 N ≥ M ðM cr ; tÞ % bðM cr À MÞ þ log 10 1 À ðt 0 =tÞ p d À1 when Σ is either stationary or maximal at t 0 (the latter should usually be the case). The event probabilities corresponding to Eqs. (8) and (9) are determined by substituting these equations into the left-hand side of Eq. (4), respectively (Fig. 3c, d).

Data availability
The data supporting the findings of this study are provided in the manuscript, in the Supplementary Information and in the Supplementary Data 1 and Data 2 files. The data and the computation details underlying Fig. 1 are provided in the Supplementary Data 1 file and in the Supplementary Information file. The map outline in the inset of Fig. 2a was obtained using https://www.generic-mapping-tools.org/. The earthquake location data underlying Fig. 2 are available as supplementary information of Kim et al. 65 and Woo et al. 20 . The hydraulic and seismic data underlying Fig. 2f are available as the supplementary information and the source data of Yeo et al. 11 . The seismic and hydraulic data underlying Fig. 3 and corresponding computation details are provided in the Supplementary Information file and in the Supplementary Data 2 file. The corresponding hydraulic data are based on the source data of Yeo et al. 11 . The corresponding seismicity catalogue is based on the open-access Report of the Korean Government Commission 18 . The data underlying Fig. 4 are collected from plots and tables of the publications on the Denver case study 4,43,44,68 . These data and the computation details are provided in the Supplementary Information file.