Interpreting the cause of bound earthquakes at underground injection experiments

Constraining the maximum possible magnitude (M_MAX) of an induced earthquake sequence is a challenging process with important implications for managing risks. CAP-tests are a suite of statistical tests that can infer, quantify, and select best-fitting M_MAX models via an earthquake catalogue's magnitudes. We use CAP-tests to discern between bound/unbound earthquake sequences at underground laboratories, where high-resolution and near-field geophysical observations are abundant. There, we find clear evidence for bound sequences, where magnitude growth was restricted during stimulation. Furthermore, bound sequences tend to be associated with stimulations that occurred within intact rock. On the other hand, unbound sequences tended to be associated with stimulations where hydraulic fractures interacted with relatively large pre-existing faults/fractures. We further examine bound sequences by fitting magnitude growth to a generalized family of M_MAX functions. This process appears to be able to aggregate bound sequences into categories consistent with theoretical considerations (e.g., tectonic, tensile-crack, or shear-crack). These results provide a basis for validating and interpreting bound sequences in controlled experiments, which is important for extrapolating to larger-scale observations. Overall, CAP-tests appear to be a promising avenue for constraining M_MAX from earthquake catalogue data.

1 Introduction

Any process that has the potential to alter stresses in the subsurface also has the potential to trigger earthquakes (Moein et al., 2023). To date, many types of anthropogenic operations have induced earthquakes including wastewater disposal, hydraulic fracturing, enhanced geothermal systems, geological carbon sequestration, reservoir impoundment, and mining (Majer et al., 2007; Foulger et al., 2018; Schultz et al., 2020). In some cases, these events have grown large enough to be felt, damaging, or even harmful (Atkinson et al., 2016; Grigoli et al., 2018). In reaction to some of these cases, social concerns have resulted in subsurface development moratoriums (Kettlety et al., 2021; Muntendam-Bos et al., 2022).

The need to manage the risks of induced earthquakes has been recognized (Bommer, 2022; Zhou et al., 2024). The de facto approach to risk management centralizes around the traffic light protocol (Bommer et al., 2006), which delineates when an operation must stop (i.e., the red-light). Typically, red-lights are designed as magnitude thresholds (Schultz et al., 2021a, b, 2023a). Despite the widespread adoption of traffic light protocols, relatively little is understood about the growth of induced earthquake magnitudes. Some recent papers have begun examining this topic; for example, by attempting to forecast the next largest event (Mendecki, 2016; Cao et al., 2020, 2024; Schultz et al., 2023b; Verdon et al., 2023; Verdon and Eisner, 2024; Yin et al., 2024). However, these approaches are unable to discern the adequacy of input models. These topics are important for understanding trailing earthquakes (Schultz et al., 2022) and magnitude jumps (Verdon and Bommer, 2021), which are among the most important factors for designing risk-informed red-lights (Schultz et al., 2021a).

A common theme underlying these efforts is constraining the maximum possible magnitude (M_MAX), given its prominence for induced seismicity hazard quantification (Bommer and Verdon, 2024). While there are numerous theoretical models that could possibly limit M_MAX (McGarr, 2014; Hallo et al., 2014; Galis et al., 2017; Elsworth, et al., 2025; Im and Avouac, 2025; Sáez et al., 2025), there are relatively few methods available to reliably validate them with empirical data (Holschneider et al., 2014; Pisarenko and Rodkin, 2022; Kijko, 2025). In fact, empirically constraining M_MAX from a catalogue is known to be a challenging task that is usually only possible in special circumstances (Holschneider et al., 2011; Zöller et al., 2016). Given the conflict between needing M_MAX constraints against the lack of robust approaches, expert solicitation is often used as a substitute (DeDontney et al., 2016). Thus, there is a need to replace subjective/opinion-based methods with more quantitative measures for M_MAX.

Based on this need, CAP-tests were designed as a suite of statistical tests to infer and quantify the presence of M_MAX within a catalogue of data. Specifically, this approach compares the distribution of earthquake magnitudes (M) against the distribution of jumps in the largest magnitude events (ΔM_LRG): when they are the same distribution, the catalogue is unbound and there is no M_MAX (and vice versa) (Schultz, 2024). Synthetic testing suggests that this problem reformulation can be more sensitive to M_MAX than traditional approaches (Schultz, 2024, 2026a). CAP-tests are broken into three component parts: the KS-test (Kolmogorov-Smirnov test) as a hypothesis test to infer the presence of M_MAX, the MLE-test (Maximum Likelihood Estimation test) to quantify the value of M_MAX, and the EW-test (Ensemble Weighting test) to select the theoretical M_MAX model that explains the catalogue data best. We provide the full methodological details of CAP-tests in Sect. 2.

The application of CAP-tests to relatively larger magnitude induced seismicity cases found that these cases behaved in an unbound manner (Schultz, 2024) – suggesting that the use of theoretical M_MAX models (McGarr, 2014; Hallo et al., 2014; Galis et al., 2017; Elsworth, et al., 2025; Im and Avouac, 2025; Sáez et al., 2025) are not appropriate during these fault reactivations. On the other hand, smaller-scale cases of induced seismicity (e.g., Utah FORGE, Preston New Road, Helsinki St1) from hydraulic stimulation indicated strong evidence for bound magnitude growth via an M_MAX (Schultz et al., 2025). Furthermore, some stimulation stages transitioned from a bound process into an unbound one; this was interpreted as these particular stages reactivating pre-existing faults instead of stimulating new fractures. That said, subsurface unknowns/uncertainties hamper an unambiguous interpretation for these field-scale cases.

To better address subsurface unknowns/uncertainties, several underground laboratories (UGLs) have performed hydraulic fracturing experiments with varying scales/scopes. These intermediate-scale experiments aim to fill a gap between laboratory studies and field-scale observations: accessing near-field geophysical observations (that are infeasible at the field-scale) in a semi-realistic in situ setting. Said another way, UGL data was collected with the intention of delineating and interpreting fracture stimulation. Thus, the complementary data from UGLs is ideal for better interpreting bound sequences (found via CAP-tests). This comparison has the potential to elucidate an underlying mechanism for bound induced seismicity sequences. This is significant, since it starts to link statistical inferences of M_MAX to physical processes.

In this study, we collect data from three UGLs: the Äspö Hard Rock Laboratory (Zang et al., 2019, 2021), the Sanford Underground Research Facility (Dobson et al., 2018; Kneafsey et al., 2018; Morris et al., 2018), and the Grimsel Test Site (Gischig et al., 2016; Amann et al., 2018). We perform CAP-tests on these UGL datasets and observe both bound and unbound stimulation stages. We cross-reference these observations against the wealth of geophysical observations and interpretations from prior multi-disciplinary studies. Bound cases tend to occur when stimulation is restricted to the growth of new fractures, while unbound cases tend to reactivate (relatively large) pre-existing faults/fractures. Furthermore, we empirically search for the functional form of M_MAX models; we delineate four categories that are predominantly consistent with theoretical expectations. Finally, we discuss how these results create a consistent picture with field-scale results.

2 Methods

2.1 Methodological Overview

The methods of this paper concern the robust inference and estimation of M_MAX from a catalogue of earthquakes. To do so, we have composed a suite of statistical tests. First, simple-tests can be used to pre-screen cases and provide qualitative or semi-quantitative inferences of M_MAX bounds. Candidate catalogues can then be examined more deeply using CAP-tests, which use a series of statistical tests to robustly discern the presence of M_MAX and quantify the best model hypothesis. Last, we outline criteria for assessing the absence of M_MAX. These tests are used to logically guide interpretations amongst the three possibilities: that M_MAX bounds a catalogue, that M_MAX is not bounding a catalogue, or that the catalogue is not well-resolved enough to make an interpretation either way. An overview of this workflow's logical sequencing has been provided (Fig. 1).

Figure 1 Workflow to discern M_MAX bounds. Simple-tests are used to pre-screen for potentially bound catalogues. Afterwards, CAP-tests can more rigorously identify bound cases.

2.2 Simple-tests

Here, we briefly define some simple-tests to provide an initial assessment of M_MAX. This is accomplished through an examination of the Gutenberg-Richter magnitude frequency distribution (GR-MFD) (Ishimoto and Iida, 1939; Gutenberg and Richter, 1944) and discrepancies between observed largest event minus the expected largest events. This sub-section constitutes the first step in our overall workflow (Fig. 1).

2.2.1 Fitting the Gutenberg-Richter Magnitude Frequency Distribution

The GR-MFD is a relationship that describes the amount and frequency of earthquake magnitudes: $N={\mathrm{10}}^a{\mathrm{10}}^{-bM}$. The a-value describes the scaling to the total amount of earthquakes N, while the b-value is the proportionality of big-to-small events. Given a catalogue of earthquake magnitudes M, this relationship can be fit through Maximum Likelihood Estimators (Marzocchi and Sandri, 2003). For real datasets, a lower bound truncation called the magnitude-of-completeness (M_c) is introduced to account for detection incompleteness. Many methods exist to evaluate M_c robustly, to account for the incompleteness in the detection of low magnitude events, thusly to avoid biases during the GR-MFD fitting process. In this study, M_c is selected by searching for the value that maximizes the goodness-of-fit metrics like R² (Schultz et al., 2018) or minimizes the negative log-likelihood of the GR-MFD. In cases with a goodness-of-fit plateaus/valleys, we select the value of M_c closest to plateau start. Sometimes we are conservative in our M_c choice by selecting a value slightly larger than optimal. Note that this conservativism will have a detrimental effect on finding bound cases. The M_c value selected by this process is typically near (but skewed right-ward of) the peak bin of the non-cumulative GR-MFD.

As part of the simple indicators of M_MAX, a visual comparison of the GR-MFD fit against the observed data is used to qualitatively examine for the presence of M_MAX. Catalogues that are relatively deficient in large magnitudes (compared to their fits) are possible candidates for an M_MAX. We note that this deficiency in large magnitude events is a key metric for identifying if/when M_MAX can be constrained (Holschneider et al., 2011; Schultz, 2024, 2026a).

2.2.2 Appraising the deficiency of large events

To provide additional semi-quantitative assessments of an M_MAX, we examine the empirical degree-of-truncation (δM_LRG), which is the discrepancy between the largest observed event minus the largest expected event. The expected M_LRG can be estimated: if a GR-MFD is assumed, then order statistics suggests a modal value of $M_{\mathrm{LRG}}=M_{\mathrm{c}}+{\log }_{\mathrm{10}}\left(N\right)/b$ (van der Elst et al., 2016). Where N is the total number of events larger than the magnitude-of-completeness M_c and b is the previously described b-value. We can use the prior GR-MFD fits to determine the expected value of M_LRG and then compare that against the observed value (i.e., δM_LRG). We also use the inverted cumulative distribution function (van der Elst et al., 2016) to determine the percentile of the δM_LRG discrepancy – or the likelihood of this degree-of-truncation occurring, assuming an unbound catalogue.

As part of the simple indicators of M_MAX, catalogues that exhibit negative δM_LRG suggest the presence of an M_MAX. We note that the δM_LRG discrepancy is a proxy metric for the resolvability of CAP-tests. Usually values of −0.5M (or less) for M_LRG−M_MAX are required to confidently assert the presence of an M_MAX (Schultz, 2024, 2026a). Note that this M_LRG−M_MAX difference is the expected value where quantitative inferences of M_MAX can start being made with 95 % confidence, from theoretical considerations (Eq. 15; Holschneider et al., 2011).

2.3 CAP-tests

The CAP-tests are a suite of statistical methods aimed at discerning the influence of M_MAX on a catalogue; each test is rooted in fundamentally different statistical frameworks, to ensure the robust validation of results. In this sense, when all the CAP/simple-tests suggest a similar bound/unbound result, we can be (more) confident that we have reached the right interpretation – even if there might be data/method issues.

In this sub-section, we explicitly define each of the CAP-tests. We also refer readers to a past study that defines the CAP-tests in detail, provides comprehensive sensitivity tests, and highlight real-data applications (Schultz, 2024; Schultz et al., 2025; Verdon and Schultz, 2026). These tests build upon the simple pre-screening assessments (Sect. 2.2), logically answering a sequence of questions to discern the influence of M_MAX more rigorously (Fig. 1).

2.3.1 The Kolmogorov-Smirnov test (KS-test)

The first test is rooted in the statistical framework of hypothesis testing. Here, we take advantage of the fact that the distribution of magnitudes (M) and the distribution of jumps in the sequence of large events (ΔM_LRG) is the same if unbound, but differ when there is an M_MAX upper bound (Schultz, 2024, 2026a). This fact ideally lends itself to hypothesis testing via the KS-test (Berger and Zhou, 2014). Given a catalogue, both M and ΔM_LRG can be observed. We can then compare these two observations against each other to test if they are drawn from the same distribution (or not), via the KS-test. This approach is advantageous in that it is non-parametric – i.e., it is completely data-driven and imposes no assumptions about the kind of distributions M or ΔM_LRG were drawn from. Because of this, we do not need to fit the data to a GR-MFD or have any knowledge/estimates of the b-value to perform our KS-test. Confidence in the KS-test is reported as compliments of standard p-values, where 95 % is a common threshold used to declare statistical significance.

Since this test is only interested in discerning the existence of an M_MAX, additional catalogue realizations can be drawn through reshuffling the order of events. In this sense, bootstrapping can be employed to generate numerous catalogue realizations in which the KS-test is repeated. This provides more robust p-value estimates. Testing on both synthetic and real datasets suggests that the KS-test is significantly more sensitive to discerning M_MAX than approaches that attempt to appraise GR-MFD fits (Schultz, 2024, 2026a). As well, synthetic testing on unbound cases shows that this KS-test produces false-positives at the rate expected for a p-value distribution (Schultz, 2024). To be able to discern the influence of M_MAX, M_LRG−M_MAX discrepancies of −0.5M or better are usually required, consistent with theoretical expectation (Holschneider et al., 2011).

While this formulation of the KS-test is powerful, we also provide a word of caution towards a potential interpretation pitfall: this method is testing for differences between the distributions M and ΔM_LRG. The presence of an M_MAX is one possible reason for this difference, but others may also confound a clear interpretation (e.g., temporal changes in b-value, kinked distributions, tapered distributions). Thus, the KS-test should be suitably pre-processed or complemented with other tests to increase the certainty of an M_MAX interpretation.

2.3.2 The Maximum Likelihood Estimator (MLE-test)

The next test is rooted in the statistical framework of Maximum Likelihood Estimation. If there is some suggestive evidence for the existence of an M_MAX, the next natural step is to quantify this value. The log-likelihood function is defined as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\ln \left(L\left(M;\mathit{\theta }\right)\right)& =\sum _i\ln \left(f_M\left(M;b,M_{\mathrm{c}},M_{\mathrm{MAX}}\right)\right)\\ & -\sum _j\ln \left(f_M\right(\mathrm{\Delta }{M}_{\mathrm{LRG}};b,\mathrm{0},M_{\mathrm{MAX}}\\ & -M_{\mathrm{LRG}}\left)\right)\end{array}\end{array}$$

where the probability density function of the GR-MFD is given by f_M(M), with a set of model parameters θ (Schultz, 2024). This function essentially entails two parts: the log-likelihood for the catalogue magnitudes M (right-hand side of equation, first term) and the log-likelihood of the jumps in largest events ΔM_LRG (right-hand side of equation, second term). The optimal set of model parameters θ are then solved for via numerical methods to maximize the log-likelihood, given the observed catalogue data M and ΔM_LRG. We perform this optimization in two-steps: the first using the standard log-likelihood constrain the b-value (Marzocchi and Sandri, 2003) (with the optimal M_c estimates) and then using the composite log-likelihood to constrain M_MAX. In this study, we consider the simplest M_MAX variant for f_M(M), which is a truncated GR-MFD (Schultz, 2024).

If this test is only interested in discerning a stationary value of M_MAX, additional catalogue realizations can be drawn through reshuffling the order of events. Similar to the KS-test, bootstrapping can be employed to generate numerous catalogue realizations in which the MLE-test is repeated. This provides more robust M_MAX estimates. Testing on both synthetic and real datasets suggests that the MLE-test is sensitive to quantifying M_MAX within a hundredth of a magnitude unit when M_LRG-M_MAX discrepancies are −0.5M or better. In cases where the MLE-test is applied to unbound catalogues, bootstrapped estimates of M_MAX will be much larger than M_LRG and standard deviations can be on the order of 1 magnitude unit or greater.

2.3.3 The Ensemble Weighting test (EW-test)

The third and final CAP-test is rooted in the statistical framework of likelihood inference. Together, the two prior tests provide suitable evidence for the existence of M_MAX. However, M_MAX may be a function of time or injected volume, becoming some non-stationary value throughout the catalogue duration. Certainly, M_MAX processes relevant for induced seismicity and hydraulic fracturing have been proposed in the past (McGarr, 2014; Hallo et al., 2014; Galis et al., 2017; Elsworth et al., 2025; Im and Avouac, 2025; Sáez et al., 2025). In this sense, having an approach that can distinguish the best proposed M_MAX model (given the data) would be insightful.

The EW-test starts by considering an ensemble of proposed M_MAX models to explain the catalogue data (Schultz, 2024). Using the previously defined log-likelihood function, both Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) can be defined for each M_MAX model (Schwarz, 1978; Akaike, 1998); note that we have used the small sample size corrections for AIC/BIC (Sugiura, 1978; McQuarrie, 1999). Next, we compute the differences in AIC/BIC scores, by subtracting the score of the best M_MAX model. These score differences can be translated into relative model weights by an exponential function (Wagenmakers and Farrell, 2004). We combine AIC/BIC weights into a single weight by taking an average between the two. In this study we will consider four standard M_MAX models: McGarr-like (i.e., log-proportional to injected volume; M_MAX∝log ₁₀(V¹)) (McGarr, 2014; Hallo et al., 2014; Elsworth et al., 2025), Galis-like (i.e., $M_{\mathrm{MAX}}\propto {\log }_{\mathrm{10}}\left(V^{\mathrm{3}/\mathrm{2}}\right))$ (Galis et al., 2017), a constant tectonic upper bound (Kanamori and Anderson, 1975), and the unbound null hypothesis. The unbound null hypothesis consists of three model parameters (M_c, b-value, GR-MFD variance), while all bound models have one additional model parameter (i.e., K=4) to account for the slope of the volume-based relationship.

The interpretation of the M_MAX model weights is straightforward: larger model weights indicate a better explanation of the data. The model with the largest weight is the best explanation of the data (within the ensemble). To quantify the statistical significance of weight differences between two models, the relative odds ratio can be computed as the ratio of the two model weights (larger/smaller). Ratios of 1+ are insignificant, 3+ are substantial/positive, 10+ are strong, and 100+ is decisive (Kass and Raftery, 1995). We note that the best model in an ensemble does not necessarily imply the veracity of the model; there could be another (unknown) model that explains the data better than all of those yet considered.

In synthetic testing, the EW-test can accurately and confidently discern the true M_MAX model (Schultz, 2024). Usually, only a handful of ΔM_LRG observations are required to confidently identify the true model (i.e., with odds ratios of 3–10 or better). However, like the prior CAP-tests, M_MAX must be influencing the catalogue for meaningful inferences to be made. Said another way, if the M_MAX is much larger than M_LRG, then the EW-test will not be able to distinguish between bound/unbound hypotheses. A schematic diagram explaining the EW-test has been provided (Fig. 2).

Figure 2 Schematic demonstration of the EW-test. The upper panels show hypothetical catalogues of earthquakes. Earthquake magnitudes (M, pink circles), the sequence of largest events (M_LRG, red line), and jumps in the largest event magnitudes (ΔM_LRG) can be directly observed. Potentially, some unobservable physical process could be bounding these catalogues (M_MAX, blue line). Left and right panels detail contrasting concepts for a bound and unbound catalogue, respectively. If there is an M_MAX, then it can be statistically inferred from the observables. Middle panels show the differences between distributions of M and M_LRG, when bounded by M_MAX (or not). The distribution of M follows a GR-MFD: analytical cumulative (thick black line) and non-cumulative (thin black line) agree with numerical cumulative (pink line) and non-cumulative (pink bars) distributions. The distribution of ΔM_LRG will differ from the GR-MFD when bounded by M_MAX. Bottom panels show the results of EW-tests using these concepts. Weights of two hypothesized models, an unbound M_MAX (pink area) and a bound M_MAX (blue area), change as new values of M_LRG are observed. EW-tests can quickly infer the presence of the true M_MAX model, from the equivocal a priori assumption (white horizontal line).

2.4 Truly unbound or just lacking data?

In an ideal experiment, injection would continue indefinitely, providing excellent resolving power for M_MAX by sampling infinitely many large events near M_MAX. Obviously, real-data cases are sample/time-limited, however. CAP-tests have the potential to discern influence from M_MAX in a (sample-limited but still) well-resolved catalogue. However, an issue arises when attempting to perform these tests on poorly-resolved catalogues, since CAP-tests cannot distinguish an M_MAX in this case. For example, if CAP-tests fail to indicate a bound catalogue, then this could be because the catalogue is either truly unbound or simply is not well-resolved enough (i.e., it has a low degree-of-truncation). Distinguishing between these two possibilities can be nebulous. We outline some guiding metrics to assist in making this assessment, which is our last workflow step in discerning the (apparent) absence of M_MAX more rigorously (Fig. 1).

The first consideration is the consistency among simple/CAP-test results. Cases that are truly bound tend to have all tests unanimously indicate a bound process, and vice versa. For reference, sensitivity testing on bound real-data that was sequentially decimated indicated a loss of resolving power, consistent with the degree of decimation (Schultz et al., 2025). Specifically, the EW-test tends to lose resolving power first, followed by the KS-test, with the MLE-test generally being the most sensitive. In fact, the use of multiple sub-tests rooted in disparate statistical methods was an intentional design choice to cover strength/deficiencies of each individual approach (Schultz, 2024). Note that this consistency indicator is only relevant for semi-well-resolved cases that are truly bound; very-poorly-resolved cases would still appear as unbound.

The second metric considered is the size of the catalogue N, within 1–2 magnitude units of the current M_LRG. Synthetic tests indicate that catalogues between 10¹–10² (above the M_c) are usually required to confidently assert M_MAX (Schultz, 2024). Similarly, real-data cases also noted that EW-tests can confidently assert an M_MAX model after observing 30–200 events – although most cases required ∼ 50–100 events (Schultz et al., 2025). Thus, we can use catalogue size as another rough indicator for how well-resolved M_MAX is.

The third metric we will consider is the currently observed M_LRG. Specifically, we will consider M_LRG in a contextual comparison against relevant cases that are already known to be bound. For example, if two independent neighbouring stages (with similarly sized catalogues) give diverging bound/unbound responses, then the observed M_LRG can be used to make inferences. Scenarios where the unbound stage has an observed M_LRG greater than the bound stage is more likely to be truly unbound. On the other hand, scenarios where the unbound stage has an observed M_LRG significantly less than the bound stage is potentially just a poorly-resolved case.

3 Data, Results, and Interpretations

We apply this workflow to data collected at a few UGLs that were aimed at better understanding the hydraulic fracturing process (Fig. 3). UGLs include the Äspö Hard Rock Laboratory (Äspö HRL), Sanford Underground Research Facility (SURF), and the Grimsel Test Site (GTS). Each of these UGLs had differing scopes, scales, and aims. Here, we will cover each UGL in a parallel style: we introduce each UGL, examine CAP-test results at the UGL, and then briefly interprets those results. We will start with the simpler UGLs and then transition into more complex cases, to pedagogically build upon the complexity of interpretations.

Figure 3 Global locations of datasets. Locations of the test sites considered: the Äspö Hard Rock Laboratory (Sweden), Sanford Underground Research Facility (South Dakota, USA), and the Grimsel Test Site (Switzerland). Supporting locations from past studies (Schultz, 2024; Schultz et al., 2025), are also shown.

3.1 The Äspö Hard Rock Laboratory

3.1.1 Overview of the Äspö HRL

The Äspö HRL is located on the Baltic east coast of Sweden, near the Simpevarp Peninsula (Fig. 3). This UGL was constructed as a testbed to study the potential for nuclear waste storage by the Swedish Nuclear Fuel and Waste Management Company starting in 1986. By 1995, construction efforts reached the 450 m target depth into the granitoids of the Trans-Scandinavian Igneous Belt (Stanfors et al., 1999).

The intention of recent (June 2015) injection experiments at the Äspö HRL was to test the feasibility of cyclic stimulation as an alternative means to enhance subsurface permeability, while also reducing the severity of induced seis micity (Zang et al., 2019, 2021). Six injection stages (HF1–HF6) were conducted in a single borehole (28 m long at 410 m depth), with an average stage spacing of ∼ 3 m, and encountering either Ävrö granodiorite or fine-grained diorite-gabbro. Stage locations were chosen to avoid natural fractures (Zimmerman et al., 2019). Each stage used between 4.1-27.2 L of injected fluid, spread between 4–6 injection sub-cycles (Zang et al., 2019). Some stages (HF1, HF2, HF4, and HF6) used traditional hydraulic fracturing techniques, while other stages (HF3 and HF5) employed cyclic stimulation (Zang et al., 2019). Resulting microseismicity was predominantly recorded in stages HF1-HF2 (López-Comino et al., 2017; Niemz et al., 2020, 2021), with events reaching up to M_w−3.5 (Kwiatek et al., 2018) and slipping with reverse or strike-slip motions (López-Comino et al., 2021). Cyclic stimulations tended to produce less seismicity (Zang et al., 2019) and more complex fracture networks (Stephansson et al., 2019; Zhuang et al., 2019), albeit with less permeability enhancement (Zimmerman et al., 2019). Data for the Äspö HRL is publicly available (Zang et al., 2024) and a spatiotemporal summary of stimulation events is plotted (Fig. 4).

Figure 4 Hydraulic stimulation and earthquake response at the Äspö HRL. (a) Timings and injection rates of six stages (blue line) are compared against induced earthquake magnitudes (circles). (b) Locations of six stages (diamonds; HF1-HF6), along the well bore (black line), and the resultant earthquakes (circles). Injection stages and corresponding earthquakes are colour coordinated in all panels.

3.1.2 CAP-tests results at the Äspö HRL

To begin assessing if some process might be restricting magnitude growth at the Äspö HRL, we fit the GR-MFD to the catalogue data from each injection stage (i.e., HF1-HF6). To account for magnitude errors, we employ a 50-trial bootstrap process in which the catalogue magnitudes are dithered by ±0.1 (this includes a dithered M_c). Only four of the stages (HF1-HF3, and HF6) have enough data to examine. Overall, the earthquakes here appear to be deficient in large-magnitude events (Fig. 5). Correspondingly, the δM_LRG discrepancies range between −0.2 to −0.4 for all stages. Assuming an unbound catalogue, these degrees-of-truncation (δM_LRG) would be <  1st, 3rd, <  1st, and 6th percentile events, respectively. These initial assessments are suggestive of some M_MAX upper bound restricting catalogue growth at the Äspö HRL.

Figure 5 Simple magnitude statistics for the Äspö HRL. GR-MFD of cumulative events (circles), alongside best fit to the data (solid lines) and the magnitude-of-completeness (dashed line). Data is colour coordinated according to injection stages (i.e., HF1–HF6).

Next, we use the CAP-tests to detect and assess the potential for M_MAX more rigorously. The KS-test is performed 50 times, in which the catalogue magnitudes are dithered (including a dithered truncation magnitude). The KS-test also performs 100 reshuffles within each trial. Average KS-test confidences are > 99.99 %, >  99.99 %, 99.87 %, and 68.85 % for stages HF1-HF3 and HF6, respectively. We perform 50 MLE-tests using dithered catalogues and 100 reshuffles within each trial. Similarly, MLE-tests suggest the standard error in fitted M_MAX values is 0.02, < 0.01, 0.02, and 0.18 for stages HF1-HF3 and HF6, respectively. Next, we seek test hypotheses for the functional form of M_MAX: our EW-tests focus on McGarr-like (McGarr, 2014; Hallo et al., 2014; Elsworth, et al., 2025), the Galis-like model (Galis et al., 2017), a constant tectonic upper bound (Kanamori and Anderson, 1975), and the unbound null hypothesis. Similarly, EW-tests also show evidence for an M_MAX bound process, most strongly evidenced with stage HF2 (Figs. 6 and S1–S3). By the end of HF2 stimulation, the Galis-like model is > 100 times more likely than the unbound model (Fig. 6).

Figure 6 Using the EW-test to discern between M_MAX models for HF2 at the Äspö HRL. (a) The chronological sequence of earthquake magnitudes (red circles), the observed M_LRG (red line), and the expected M_LRG at the 10/50/90 percentiles (red area) are compared against M_MAX hypotheses (blue lines). (b) The corresponding ensemble weights update as new data is encountered.

Last, we organize our results for all injection stages at the Äspö HRL – for convenience to the reader. These results are summarized below (Table 1).

Table 1 Summary of results at the Äspö HRL. All the prior results of our simple-tests and CAP-tests are compiled here for convenience. Additionally, we have coordinated individual entries according to their interpretation: blue/bold for bound, pink/italic for unbound, and uncoloured/regular for indeterminate.

3.1.3 Interpretations for the Äspö HRL

All the stages at the Äspö HRL appear to have some evidence of a bound process restricting the growth of earthquake magnitudes. The most convincing of the cases is also the most well-resolved one: for HF2, all the simple/CAP-tests unanimously indicate a bound process with strong statistical significance. As the cases become less well-resolved, the statistical confidence also diminishes. For example, the HF1 and HF3 stages have most CAP-tests indicating a bound process, only the EW-test of HF3 narrowly falls short of statistical significance. The HF6 stage produces the most mixed results: simple-tests, MLE-test, and EW-test are indeterminate, while the KS-test is unbound. It is worthwhile to mention that HF6 differed in stimulation approach in two ways: the first sub-stage of HF6 was the largest injected volume (with reduced volume for later sub-stages) and the last three sub-stages were performed on a different day.

Overall, we interpret stage HF2 as certainly being bound. Stages HF1 and HF3 as likely to be bound, but likely needing more well-resolved catalogues to better discern M_MAX. We hesitate to make a clear statement for HF6, which appears indeterminate due to data limitations.

The interpretation that all the Äspö HRL stages were (likely) bound corresponds with the geophysical interpretations from relevant studies. The intention of the Äspö HRL was to test various stimulation programs against the complexity/growth of hydraulic fractures (Zang et al., 2019, 2021). Stage intervals were chosen to avoid natural fractures, and impression packers noted the generation of new hydraulic fractures (Zimmerman et al., 2019). Furthermore, the progressive growth of the HF2 hydraulic fracture plane was inferred jointly from the microseismic and deformation constraints (Niemz et al., 2020, 2021). In this sense, a progressively growing fracture aligns well with the bound interpretation: the finite extent of a fracture limits M_MAX via geometric considerations (Kanamori and Anderson, 1975). As stimulation continues, the fracture continues to grow; thus, the value of M_MAX would increase alongside the injected volume. Correspondingly, each of the Äspö stages indicated bound growth (albeit with varying degrees of confidence).

These interpretations at the Äspö HRL constitute the simplest interpretation. We describe a scenario as to how hydraulic fracturing would be linked to a bound M_MAX interpretation. Regardless of the fracture network's complexity, the finite spatial extent of stimulated fractures ultimately restricts magnitude growth.

3.2 The Sanford Underground Research Facility

3.2.1 Overview of the EGS Collab Experiment 1

The SURF is located near Lead, South Dakota, and has repurposed the Homestake gold mine (Fig. 3). The SURF is a research facility operated by the South Dakota Science and Technology Authority to study rare-process physics (Heise, 2015). The EGS Collab takes advantage of the SURF, by using this facility to host hydraulic stimulation experiments at field-scale depths of ∼ 1.5 km (Dobson et a., 2018; Kneafsey et al., 2018; Morris et al., 2018).

The EGS Collab Experiment #1 intended to connect injection and production boreholes via stimulated fractures in a controlled environment (Kneafsey et al., 2020; Morris et al., 2018). Injection and production boreholes were drilled subparallel to the estimated minimum principal stress direction (Oldenburg et al., 2017) in the phyllites of the Precambrian-aged Poorman Formation (Kneafsey et al., 2020). These metamorphic rocks are strongly foliated and highly anisotropic (Frash et al., 2019; Vigilante et al., 2017). The rock mass is generally thought to be low permeability ($\sim {\mathrm{10}}^{-\mathrm{18}}$ m²), with a thermally altered stress field (Singh et al., 2019). Some natural fractures were noted in borehole cores (Fu et al., 2021), with at least one noteworthy fracture that is naturally/actively flowing (Wu et al., 2021a). Starting in May 2018, stimulation stages took place at three wellbore intervals (N164, N142, N128), covering both hydraulic stimulation and hydraulic characterization programs; stimulation programs injected on the order of 10 s L of water per interval (Morris et al., 2018; White et al., 2019). Resulting microseismicity was recorded for each interval, with ∼ 2000 located events (Schoenball et al., 2019, 2020; Chai et al., 2020; Qin et al., 2024). Studies have covered topics including modeling heat transport (Wu et al., 2021b), strain/deformation (Guglielmi et al., 2021), and hydraulic fracture propagation (Li and Zhang, 2023). Data for the EGS Collab Experiment #1 is publicly available (https://gdr.openei.org/, last access: July 2025) and a spatiotemporal summary of stimulation events is plotted (Fig. 7).

Figure 7 Hydraulic stimulation and earthquake response for SURF Experiment #1. (a) Timings and injection rates of three stages (blue line) are compared against induced earthquake magnitudes (circles). (b) Locations of three stages (diamonds; N164, N142, N128), along the well bore (black line), and the resultant earthquakes (circles). Injection stages and corresponding earthquakes are colour coordinated in all panels.

3.2.2 CAP-tests results from the EGS Collab Experiment 1

To begin assessing if some process might be restricting magnitude growth at SURF Experiment #1, we fit the GR-MFD to the catalogue data from each injection stage (i.e., N164, N142, N128). To account for magnitude errors, we employ a 50-trial bootstrap process in which the catalogue magnitudes are dithered by  ± 0.1 (this includes a dithered M_c). Only one of the stages (N164) appear to be appreciably deficient in large magnitude events (Fig. 8). Correspondingly, the δM_LRG discrepancy is significant for N164 (−1.15M), but fairly small for N142 (−0.05M) and N128 (−0.28M). Assuming an unbound catalogue, the δM_LRG discrepancy for N164 would be a < 1st percentile event; stages N142 (10th percentile) and N128 (33rd percentile) are more commonly expected occurrences. These initial assessments are suggestive of some M_MAX upper bound restricting catalogue growth for one stage at SURF Experiment #1.

Figure 8 Simple magnitude statistics for SURF Experiment #1. GR-MFD of cumulative events (circles), alongside best fit to the data (solid lines) and the magnitude-of-completeness (dashed line). Data is colour coordinated according to injection stage (i.e., N164, N142, N128).

Next, we use the CAP-tests to detect and assess the potential for M_MAX more rigorously. The KS-test is performed 50 times, in which the catalogue magnitudes are dithered (including a dithered truncation magnitude). The KS-test also performs 100 reshuffles within each trial. Again, KS-test results are split, with average confidences of 99.96 % for stage N164, but 34.23 % and 74.39 % for stages N142 and N128, respectively. We perform 50 MLE-tests using dithered catalogues and 100 reshuffles within each trial. MLE-tests are also split: standard error in fitted M_MAX values are 0.10, 2.18, and 1.31 for stages N164, N142, and N128, respectively. Similarly, EW-tests also shows suggestive evidence for an M_MAX bound process with N164 (Fig. 9), but certainly not for stages N142 and N128 (Figs. S4–S5 in the Supplement).

Figure 9 Using the EW-test to discern between M_MAX models for N164 for SURF Experiment #1. (a) The chronological sequence of earthquake magnitudes (red circles), the observed M_LRG (red line), and the expected M_LRG at the 10/50/90 percentiles (red area) are compared against M_MAX hypotheses (blue lines). (b) The corresponding ensemble weights update as new data is encountered.

Last, we organize our results for all injection stages at SURF Experiment #1 – for convenience to the reader. These results are summarized below (Table 2).

Table 2 Summary of results for the SURF EGS Collab Experiment #1. All the prior results of our simple-tests and CAP-tests are compiled here for convenience. Additionally, we have coordinated individual entries according to their interpretation: blue/bold for bound, pink/italic for unbound, and uncoloured/regular for indeterminate.

3.2.3 Interpretations for the EGS Collab Experiment 1

CAP-test results at SURF Experiment #1 have similarities and differences from those at the Äspö HRL. For example, stage N164 is similar to the stages at the Äspö HRL: both simple-tests and CAP-tests unanimously agree on an M_MAX bound – with varying degrees of confidence. On the other hand, stages N142 and N128 differ in that they produce strong and unambiguous unbound inferences. Because of this, we interpret stages N142 and N128 to be truly unbound.

The interest in considering UGL cases is the wealth of complementary geophysical information to cross-examine against the results of CAP-tests. The stimulation at stage N164 is predominantly understood to be the creation of a new fracture network. This interpretation comes from multiple lines of evidence: the orientation of microseismic fault planes with respect to the ambient stress field (Schoenball et al., 2020), deformation constraints on fracture motion (Guglielmi et al., 2021), significant recovery of injected fluid (White et al., 2019), and direct evidence of new fluid jets intersecting the producing well (Fu et al., 2021). The later N164 injections (i.e., June 2018) were performed from the producing well side, into the newly created fracture network. That said, some complexity in the N164 stimulation suggested limited interaction with natural fractures, via either arrests or redirected continuation of fracture growth (Schoenball et al., 2020; Fu et al., 2021). Overall, these inferences are consistent with finite hydraulic fractures bounding the growth of earthquake magnitudes – which are the same interpretations made for the Äspö HRL.

On the other hand, stages N142 and N128 had significant interactions with natural pre-existing fractures, and likely reactivated them in shear rather than creating a new fracture network. These inferences/interpretations come from the orientation of microseismic fault planes with respect to the ambient stress field, corroboration with well image logs, and observations of shear deformation (Schoenball et al., 2020). The connection to fractures/faults allows for seemingly unbound magnitude growth, as the spatial extent of these pre-existing structures have the potential to host larger events.

These interpretations at SURF Experiment #1 constitute the next level of complexity in interpretation. We describe a scenario where hydraulic stimulation could interact with natural fractures/faults to produce an unbound process. Said another way, the interaction with natural fractures/faults (of sufficient size) can facilitate unbound magnitude growth.

3.3 The Grimsel Test Site

3.3.1 Overview of the GTS

The GTS is located near the Grimsel pass in the central Swiss Alps (Fig. 3), ∼ 450 m below the Juchlistock in the Varsican-aged Aar Massif granites (Schneeberger et al., 2019). The GTS was established in 1984 as an underground research facility to study the safe disposal of nuclear waste. It is owned and operated by the National Cooperation for Radioactive Waste Disposal (Nagra; https://www.grimsel.com/, last access: March 2025).

The intention of recent (February–May 2017) injection experiments at the GTS was to demonstrate the stimulation of fractures at the decameter scale and to better understand how to manage induced seismicity (Gischig et al., 2016; Amann et al., 2018). Twelve injection stages spanning ∼ 1 m intervals were situated in two boreholes (∼ 45 m long) (Gischig et al., 2020): six related to stimulation via hydraulic fracturing (HF1-HF8) (Dutler et al., 2019) and six more via hydroshearing (HS1-HS8) (Krietsch et al., 2020b). Each stage used (on the order of) ∼ 1000 L of injected fluid, spread between 4 injection sub-cycles. The role of natural faults and fractures are a prominent focus of the GTS injection experiments, which injected into either brittle-ductile shear zones, ductile shear zones, or intact rock (Doetsch et al., 2018a). The degree of seismic response for individual stages was strongly heterogeneous in space, with stages HS4, HS5, and HF2 being both the most seismically active and hosting the large events (M−3.0, M−2.8, and M−3.9) (Villiger et al., 2020). Studies at the GTS were diverse, covering topics like geological characterization (Krietsch et al., 2018), stress inversion (Krietsch et al., 2019), tomographic velocity changes (Doetsch et al., 2018b; Schopper et al., 2020), permeability changes (Jalali et al., 2018; Brixel et al., 2020a, b), and inferring fracture propagation from hydromechanical response (Dutler et al., 2019; Krietsch et al., 2020a, b). Data for the GTS is publicly available (https://doi.org/10.3929/ethz-b-000276170, Doetsch et al., 2018c) and a spatiotemporal summary is plotted (Fig. 10).

Figure 10 Hydraulic stimulation and earthquake response at the GTS. (a, b) Timings and injection rates of twelve stages (blue line) are compared against induced earthquake magnitudes (circles). (c) Locations of twelve stages (diamonds; HFX and HSX), along the well bore (black line), and the resultant earthquakes (circles). Injection stages and corresponding earthquakes are colour coordinated in all panels.

3.3.2 CAP-tests results at the GTS

To begin assessing if some process might be restricting magnitude growth at the GTS, we fit the GR-MFD to the catalogue data from each injection stage (i.e., HSX & HFX). To account for magnitude errors, we employ a 50-trial bootstrap process in which the catalogue magnitudes are dithered by ±0.1 (this includes a dithered M_c). Note that many of the stages recorded here have too few events for a meaningful analysis; thus, we predominantly focus our discussions to a subset of stages. Only two of the stages (HS4 and HF2) appear to be appreciably deficient in large magnitude events (Fig. 11). Correspondingly, the δM_LRG discrepancy is large for HS4 and HF2 (−1.31M and −0.61M), but fairly small for the other viable stages (between −0.19M and +0.08M). Assuming an unbound catalogue, the δM_LRG discrepancy for HS4 and HF2 would both be < 1st percentile events; all other stages are commonly expected occurrences. These initial assessments are suggestive of some M_MAX upper bound restricting catalogue growth for just two stages at the GTS.

Figure 11 Simple magnitude statistics at the GTS. GR-MFD of cumulative events (circles), alongside best fit to the data (solid lines) and the magnitude-of-completeness (dashed lines). Data is colour coordinated according to injection stage (i.e., HSX and HFX) and separated into hydroshearing (a; HSX) or hydrofracturing experiments (b; HFX).

Next, we use the CAP-tests to detect and assess the potential for M_MAX more rigorously. The KS-test is performed 50 times, in which the catalogue magnitudes are dithered (including a dithered truncation magnitude). The KS-test also performs 100 reshuffles within each trial. KS-test results are split, with strong confidences of > 99.99 % for stages HS4 and HF2, but unconvincing values for all other stages. We perform 50 MLE-tests using dithered catalogues and 100 reshuffles within each trial. MLE-tests are also split: standard error in fitted M_MAX values are < 0.01 and 0.02 for stages HS4 and HF2, but values are large (0.34–1.92) for all other stages. The EW-tests shows decisive evidence for an M_MAX bound process with HS4 (Fig. 12), and strong evidence for HF2 (Fig. S6). All other stages have substantial-to-strong evidence for an unbound process via EW-tests (Figs. S7–S10).

Figure 12 Using the EW-test to discern between M_MAX models for HS4 at the GTS. (a) The chronological sequence of earthquake magnitudes (red circles), the observed M_LRG (red line), and the expected M_LRG at the 10/50/90 percentiles (red area) are compared against M_MAX hypotheses (blue lines). (b) The corresponding ensemble weights update as new data is encountered.

Last, we organize our results for all injection stages at the GTS – for convenience to the reader. These results are summarized below (Table 3).

Table 3 Summary of results at the GTS. All the prior results of our simple-tests and CAP-tests are compiled here for convenience. Additionally, we have coordinated individual entries according to their interpretation: blue/bold for bound, pink/italic for unbound, and uncoloured/regular for indeterminate.

3.3.3 Interpretations for the GTS

In summary, most of the GTS stages appear to exhibit unbound growth of earthquake magnitudes. Although, this interpretation varies in confidence, depending on the specific stage in question, since many of the stages had relatively few events recorded. The most confidently unbound case is HS5, which fails all of the simple-tests and CAP-tests. On the other hand, there are two exceptions to this general trend: HS4 shows clear and definitive evidence for a bound process, while HF2 shows strong evidence for a bound process. That said, stages with fewer than 50 events are difficult to arrive at a clear interpretation between truly unbound or simply lacking data.

The aims and scope of the GTS were to investigate the response of injection into faults/fractures at the intermediate scale (Gischig et al., 2016; Amann et al., 2018). For example, all HSX stages intentionally injected into previously known fractures of shear damage zones (Krietsch et al., 2018); thus, micro-seismically delineated fault planes were generally consistent with fracture/fault orientations (Villiger et al., 2020), near-field hydromechanical effects were consistent with pre-existing fracture dislocation (Krietsch et al., 2020b), and many stages were hydraulically connected to other boreholes via these stimulated fractures (Brixel et al., 2020b). In particular to HS5, seismicity was well fit to a single fault plane (∼ 16 m diameter) oriented subparallel to the targeted fracture (Villiger et al., 2020), significant pressure increases (∼ 70 %–75 % of injection pressure) were observed in boreholes 7–8 m away (Krietsch et al., 2020a), likely due to transient permeability increases driven by fracture aperture changes (Krietsch et al., 2020b). Furthermore, the diameter of the HS5 fault plane (∼ 16 m) is appreciably larger than the largest event diameter (∼ 0.65 m), assuming a circular crack. HFX stages started hydraulically fracturing intact rock and subsequently propagated/stimulated new fractures (Dutler et al., 2019). However, the newly propagating fractures were inferred to have significant interactions with pre-existing fractures, which served as pressure sinks that arrested further growth (Dutler et al., 2019; Villiger et al., 2020). Note that HF6 is unique, since it injected directly into a pre-existing fracture (by mistake) and thus can be considered a hydroshearing experiment. Furthermore, the propagation of microseismic events along pipe-like geometries was thought to be formed via the intersection with natural fractures (Dutler et al., 2019). Given the prior interpretations of CAP-test results at SURF Experiment #1, we would anticipate that all GTS stimulations should be unbound – because of their prominent connection to pre-existing fractures. Bound cases would require special exceptions to this generalization.

On the other hand, the HS4 stage is a clear exception to this generalization. The HS4 seismicity is not well-fit by a single plane. Instead, it is best fit by four intersecting planes (C1–C4), where C1–C3 are oriented subparallel with pre-existing fractures and C4 is likely a new tensile failure (Villiger et al., 2020, 2021). Clustering in focal mechanism slip style also corresponds to spatial clusters (Villiger et al., 2021). Each cluster has seismicity spatially restricted along discrete linear streaks (Villiger et al., 2021); these streaks grow/propagate alongside injection, although their spatial extent (∼ 1–2 m) is much smaller than other HSX fault planes. In tectonic contexts, streaks are often interpretated as rheological boundaries between seismic and creeping/locked fault segments (Rubin et al., 1999; Waldhauser et al., 2004). In other hydraulic fracturing contexts, streaks have also been observed and analogously interpreted (Rutledge et al., 2004; Evans et al., 2005). Thus, the streaks and clustering at HS4 have been interpreted as fractures channelizing fluid-flow towards highly seismogenic asperities that slip perpendicular to the fluid migration direction (Villiger et al., 2021). In this sense, we argue that HS4 seismicity is spatially bound to these asperities/streaks, giving rise to the bound growth of magnitudes. Correspondingly, the GR-MFD for HS4 (Fig. 11a) starts to roll-off around M−4.3 and M−3.0 was the final M_LRG value; note that M−3.0 roughly corresponds to a circular crack diameter of ∼ 0.5–1.0 m, which is comparable to the spatial extent of the streaks (∼ 1–2 m).

Similarly, HF2 is also an exception to this generalization. The HF2 seismicity was best fit by two intersecting planes, the first which resembles a newly created hydraulic fracture and then a deflection to a more E–W orientation (Dutler et al., 2019; Villiger et al., 2020). It has been suggested that this stage is exceptional in that it potentially has limited interaction with pre-existing structures, being able to propagate fracture growth before leak-off into potential fracture connections (Dutler et al., 2019; Villiger et al., 2020). In fact, the growing hydraulic fracture intersected monitoring boreholes, reaching further than suggested by the microseismic events (Dutler et al., 2019). In this sense, we would argue that HF2 is most similar to the N164 stage at SURF Experiment #1. Correspondingly, both HF2 and N164 express a bound process via CAP-tests.

These interpretations at the GTS constitute the greatest level of interpretive complexity we will consider in this study. We have reiterated scenarios where hydraulic stimulation could interact with natural fractures/faults to produce an unbound process. Regardless of the stimulation program, the interaction with natural fractures/faults (of sufficient size) can facilitate unbound magnitude growth. Exceptions come in the form of hydraulic fracturing with limited interactions to natural fractures (HF2) and shear reactivation with seismicity limited to streaks/asperities (HS4).

4 Discussion

4.1 Consolidating M_MAX interpretations

In this sub-section, we synthesize the results and interpretations of this study in comparison with relevant literature. We begin to jointly interpret the physical rationale behind bound/unbound earthquake sequences. We start by reiterating the interpretations of CAP-tests results at the UGLs and then discuss the correspondence with CAP-tests at the field-scale.

In the simplest case of the Äspö HRL (Sect. 3.1), stages were stimulated among intact rock, creating new hydraulic pathways via stimulated fractures. The finite geometric extent of these stimulated fractures limits the size of the largest events that can occur, at any given instant (Fig. 13). Although, the spatial extent of this fracture-bound M_MAX will change with time/injection. Correspondingly, all Äspö HRL cases indicated a bound M_MAX process via CAP-tests. While newly stimulated fractures were the intention of the Äspö HRL, bound cases were also observed for analogous stages at SURF Experiment #1 (N164; Sect. 3.2) and the GTS (HF2; Sect. 3.3), despite their different experiment objectives.

Figure 13 Interactions between stimulated fractures and reactivated faults during hydraulic fracturing. A series of stages of increasing complexity (text labels) are considered alongside a hypothetical well (black and grey rectangle). In the simple case, the host rock is split in tension via stimulated fractures (grey polygons). In the complex case, small pre-existing faults (black polygons) can also be reactivated in shear slip. The clustered case hydraulically connects a series of stages/perforations into a single fracture/fault network. In the last case, stimulated fractures intersect a large fault system that is critically-stressed for slip. By watching the growth of earthquake magnitudes, CAP-tests can discern between bound cases (blue circles) and unbound cases (pink circles).

On the other hand, stages encountering more complicated interactions with pre-existing fractures/faults tended to exhibit an unbound catalogue. For example, the simplest case of fracture interaction was noted at SURF Experiment #1 (Sect. 3.2): the two stages (N142 and N128) that intersected a natural fracture were also unbound via CAP-tests. Furthermore, injection into faults/fractures was an aim of the GTS (Sect. 3.3), and these stages predominantly exhibited an unbound process. We interpret this correspondence via geometrical considerations again (Fig. 13): an immediate hydraulic connection to a relatively large-scale fracture system provides an opportunity for (apparently) unrestricted magnitude growth. While we have not observed it at these three UGLs (to our knowledge), interactions with smaller-scale pre-existing fractures would be expected to produce a bound sequence.

It is worth emphasizing that stage HS4 at the GTS is exceptional (Sect. 3.3), in that a bound M_MAX process was observed for a case with strong natural fracture interactions. Thus, this stage is an outlier to the general interpretive theme regarding natural faults/fractures. That said, there is also good evidence for an alternative bounding mechanism at HS4. Earthquakes were restricted to smaller streaks along larger fault plane trends (Villiger et al., 2021); these streaks grew in spatial extent with increased injection volume. Streaks are interpreted as fault heterogeneity, where seismic asperities are surrounded by a broader creeping/locked fault segment (Rubin et al., 1999; Waldhauser et al., 2004; Rutledge et al., 2004; Evans et al., 2005). Furthermore, the spatial extent of these streaks is comparable to the fault plane area of the largest events observed. In this sense, HS4 seismicity is likely bound to these asperities/streaks, giving rise to the bound growth of magnitudes. Following this logic, it should be possible to test seismic asperities for a bound M_MAX process (via CAP-tests) in tectonic settings too.

A recent study performed CAP-tests on field-scale experiments (Schultz et al., 2025), including both shale gas hydraulic fracturing at Preston New Road in the UK (Clarke et al., 2019; Kettlety et al., 2021) and the enhanced geothermal system at Utah FORGE (Moore et al., 2019; Niemz et al., 2025). CAP-tests indicated that many clustered stages followed a bound process, while the more hazardous stages were unbound. Interpretations in these field-scale cases followed a similar interpretation to UGL stages in this study (Sect. 4.1.1), albeit with larger uncertainties. This interpretation was partly driven by clustering sensitivity tests that showed omitting/ignoring earlier bound stages tended to diminish statistical confidence (Schultz et al., 2025) – providing an inference to the importance of accurate clustering and hydraulic connectivity for bound sequences. Certainly, natural fractures have demonstrably influenced the propagation of stimulated fractures (Jin et al., 2024). In this sense, the verification of this interpretation at UGLs provides a potential conceptual linkage to the field-scale, given their observational similarities. Prior results for bound processes are consistent with stimulating a fracture system with a restricted extent, while unbound cases likely have reactivated larger fault systems (Schultz, 2024).

4.2 Empirically constraining the functional form of M_MAX vs. V

Up to this point, our study has focused on discerning between bound/unbound sequences – and establishing the importance of natural fault/fracture systems within this dichotomy. However, our discussion has neglected the underlying physical mechanisms responsible for creating this bound M_MAX process. Here, we further discuss the underlying physical mechanisms, for all the known bound sequences with data available (Fig. 3) – both in this study and past studies (Schultz et al., 2025).

Numerous theoretical models of M_MAX and fracture propagation have been proposed in the literature (Eaton and Igonin, 2018). The simplest model considers a physical limitation based on the geometry of the finite fault extent (Kanamori and Anderson, 1975); For a circular fault with radius R and constant stress drop Δσ, the seismic moment release is given via the equation $M_{\mathrm{0}}=\frac{\mathrm{16}}{\mathrm{7}}\mathrm{\Delta }\mathit{\sigma }{R}^{\mathrm{3}}$. Seismic moment can then be translated into an M_MAX magnitude via the moment magnitude relationship (Hanks and Kanamori, 1979). Non-stationary M_MAX models have been proposed for induced seismicity based on the footprint of stress perturbation along a fault (Shapiro et al., 2011), the evolution of seismic moment (McGarr, 2014; Hallo et al., 2014; Elsworth et al., 2025), self-arrested rupture dynamics (Galis et al., 2017), or time-dependent pressure diffusion (Shapiro et al., 2021). Additionally, we can consider the finite extent of fracture propagation via an aseismic shear crack (Danré et al., 2024) or a tensile crack (Davis et al., 2020). These two models suggest time-varying fracture radii based on the equations $R\left(t\right)=\sqrt{s_{d}V\left(t\right)}$ and $R\left(t\right)=\sqrt[\mathrm{3}]{kV\left(t\right)}$, respectively. Here, s_d and k are injection propagation coefficients. M_MAX models have been investigated in numerous prior studies (Kwiatek et al., 2018; Li et al., 2021; Yu et al., 2024; Langenbruch et al., 2024; Lanza et al., 2026).

Examining the theoretical models illustrates a common mathematical theme: seismic moment release (M₀) is proportional to injected volume, raised to some exponent. For example, the McGarr-like model considers M₀∝V¹ (McGarr, 2014), while the Galis-like model considers $M_{\mathrm{0}}\propto V^{\mathrm{3}/\mathrm{2}}$ (Galis et al., 2017). Similarly, a growing tensile crack is equivalent to M₀∝V¹ (Kanamori and Anderson, 1975; Davis et al., 2020), while a growing aseismic shear crack is $M_{\mathrm{0}}\propto V^{\mathrm{3}/\mathrm{2}}$ (Kanamori and Anderson, 1975; Danré et al., 2024; Im and Avouac, 2025; Sáez et al., 2025). Even purely tectonic/geometric limitations can be cast into this form as M₀∝V⁰. Given this commonality, we consider a generic family of M_MAX models that have a functional form of M_MAX∝log ₁₀(Vⁿ). This generalization is advantageous, as the EW-test can consider any arbitrary M_MAX function. Following this logic, we use these arbitrary functions within the EW-test to find the best-fitting n-exponent for each bound case. To ensure a fair comparison, each competing/candidate M_MAX model is first fit to the case data before the EW-test selects the best model.

This investigation indicates that the underlying physical mechanisms are varied among cases (Fig. 14). However, the results tend to aggregate into four distinct categories: geometric/streak-like (n 0.0–0.3), X-like (n 0.5–0.7), McGarr/tensile-like (n 0.9–1.1), and aseismic/shear-like (n 1.4–1.7). Bounding mechanisms can be varied, even between similar stages (e.g., Äspö HRL, HF1-HF3). Higher n-exponents tend to be less confident, likely because aseismic/shear-like exponents emulate/approach unbound-like magnitude growth (van der Elst et al., 2016; Galis et al., 2017). Other geophysical data/studies to corroborate the inferences of physical mechanisms are generally limited to a few cases. For example, the propagation of microseismic events at Utah FORGE stages 1–3 were best-fit by tensile-like models (Clarkson et al., 2025; Lanza et al., 2026; Niemz et al., 2026), which independently corresponds to our best-fit n-exponents (0.97 and 1.06) for clusters 1 (stages 1–2) and 2 (stages 3–6) (Fig. 14), as defined in Schultz et al. (2025). Interestingly, the low n-exponent associated with HS4 at the GTS is consistent with the observation of tectonic streaks/asperities that exhibited limited spatial growth during injection (Villiger et al., 2021). The X-like category is exceptional, since there is (to our knowledge) no previously suggested model for this n-exponent. This could be potentially explained by M_MAX models that do not fit the Vⁿ generalization (Shapiro et al., 2011, 2021), an as-of-yet undiscovered M_MAX model, or observational biases for encountering a relatively large pre-existing fracture network (V⁰) that also starts growing alongside injection (V¹ or $V^{\mathrm{3}/\mathrm{2}}$). Currently, we favour the observational bias interpretation for the X-like category, based on the results of drop-out tests that demonstrate how an intermediate sized pre-existing fracture network that also grows alongside injection can produce a mixed n-exponent with poorly-resolved data (Supplement Sect. S1 and S16). We consider it likely that these considerations were a factor for the SURF N164 and Äspö HF3 cases, which indicated the next closest M_MAX model when restricting choices to our three bound options (Tables 1 and 2). Accounting for some prior growth appears to be an important consideration for discerning the true M_MAX model.

Figure 14 Searching for volume-exponents. The EW-test fit to case data (lines) for both (a) normalized model weights and (b) odds ratios (relative to the unbound null hypothesis). Strength of statistical confidence and n-exponent categories are shown with grey shaded areas (and text callouts). Note that relative odds ratios are capped at 5000 to facilitate easier comparison via plotting.

Additional UGL study alongside refinements to CAP-tests (potentially also considering other candidate M_MAX models) will likely shed light on the physical processes underlying bound sequences. For example, CAP-tests could be used to systematically search for M_MAX models in other datasets. We note that this approach could be straightforwardly adapted to forecast the “next largest event” (Cao et al., 2020), using these weights as part of either a logic-tree (Bommer and Verdon, 2024) or an ensemble of models (Schultz et al., 2023c).

4.3 Exceptions to bound/unbound interpretations

Thus far, we have predominantly interpreted the difference between bound/unbound sequences via distinctions between fracture stimulation and fault reactivation (Fig. 13). While this appears to capture the first-order effects, there likely exists subtleties that can deviate from this conjecture.

4.3.1 Exceptions to the bound-fracture interpretation

We have predominantly interpreted bound sequences arising from the finite (and time/volume dependent) extent of a propagating fracture restricting available slip area. While this interpretation matches both expectation and observation, we can still envision a scenario in which a finite fracture could produce an unbound sequence. For example, in the case that the finite extent of the fracture grows faster than the equivalent area of an unbound M_LRG. In this sense, while the sequence is technically bound, it would never be empirically inferable by magnitude information. Criteria for this condition are given below: for the Vⁿ generalization of M_MAX as well as the tensile-like and shear-like fractures. There, the expected M_LRG formulation is a combination of the population expectation (van der Elst et al., 2016) coupled with a proportionality between injected volume and event counts (Shapiro et al., 2010). Note that Σ is the seismogenic index, a parameter indicating a fault's potential to induce earthquakes per unit of injected fluid volume (Shapiro et al., 2010).

$$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}\left({\displaystyle \frac{\mathrm{2}n}{\mathrm{3}}}-{\displaystyle \frac{\mathrm{1}}{b}}\right)& {\log }_{\mathrm{10}}\left(V\left(t\right)\right)+\left({\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\right){\log }_{\mathrm{10}}\left(c\right)\\ & -\mathrm{9.1}\gg {\displaystyle \frac{\mathrm{\Sigma }}{b}}\left(V^n\mathrm{generalization}\right)\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{b}}\right)& {\log }_{\mathrm{10}}\left(V\left(t\right)\right)+\left({\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\right){\log }_{\mathrm{10}}\left({\displaystyle \frac{\mathrm{16}}{\mathrm{7}}}\mathrm{\Delta }\mathit{\sigma }{s}_d^{\mathrm{3}/\mathrm{2}}\right)\\ & -\mathrm{9.1}\gg {\displaystyle \frac{\mathrm{\Sigma }}{b}}\left(\mathrm{tensile}\mathrm{fracture}\right)\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\left({\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}-{\displaystyle \frac{\mathrm{1}}{b}}\right)& {\log }_{\mathrm{10}}\left(V\left(t\right)\right)+\left({\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\right){\log }_{\mathrm{10}}\left({\displaystyle \frac{\mathrm{16}}{\mathrm{7}}}\mathrm{\Delta }\mathit{\sigma }k\right)\\ & -\mathrm{9.1}\gg {\displaystyle \frac{\mathrm{\Sigma }}{\mathrm{b}}}\left(\mathrm{shear}\mathrm{fracture}\right)\end{array}\end{array}$$

As a reminder of terms: V(t) is the injection volume time-series, n is the volume exponent, b is the b-value, Δσ is the stress drop, while c, s_d and k are injection propagation coefficients for different growth.

4.3.2 Exceptions to the unbound-fault interpretation

We have predominantly interpreted unbound sequences arising from fractures interacting with relatively-large pre-existing faults. While this matches our expectation and observations, we can still envision scenarios where pre-existing faults could produce bound sequences. For example, M_LRG could be limited by the pre-existing faults being too small, the in situ stress resolved on the fault (Gischig, 2015; Norbeck and Horne, 2018), time-dependence of stress perturbations (Segall and Lu, 2015), fault geometry/complexity (Lee et al., 2024), fault roughness (Maurer et al., 2020; Wang et al., 2024), material/stress heterogeneity (Kroll and Cochran, 2021), rheology, or the extent of asperities. In fact, HS4 at the GTS is one such exception: where streaks outlined the rheological boundary between seismic asperities hosted within an aseismic fault. If these factors are limiting M_LRG, CAP-tests should be able to discern their effects.

5 Conclusions

In summary, we have rigorously tested CAP-tests against data from controlled injection experiments at three underground laboratories. Bound sequences are consistent with observations of new hydraulic fracturing growth, while unbound typically reactivate larger pre-existing structures. Furthermore, the EW-test appears to be able aggregate generalized M_MAX functions into categories consistent with theoretical considerations for tectonic/geometric limits, tensile/McGarr-like growth, or aseismic/shear-like growth. This process is potentially able to identify the underlying physical mechanism responsible for M_MAX, specific to an individual earthquake sequence. Overall, the ability to robustly identify and constrain M_MAX will likely be important for both induced seismicity and natural/tectonic earthquake sequences.