In this study, we evaluated the damage potential of shallow, small-to-moderate-magnitude earthquakes (Mw 3–5) associated with induced seismicity on the seismic performances of earth canal dykes, which are common geotechnical systems in regions affected by induced seismicity. Increasing earthquake rates driven by hydraulic fracturing and wastewater injection have enhanced the seismic demands on infrastructure in areas with historically low natural seismicity (Petersen et al., 2018). Although induced earthquakes rarely reach magnitudes capable of causing catastrophic damage, they can still cause significant damage to engineering systems, such as earth canal dykes, along with substantial economic consequences. For example, Gupta and Baker (2019) estimated that the average annual expected loss would be approximately $125 million in Oklahoma owing to induced earthquakes. In 2016, the United States Geological Survey (USGS) published the first 1-year ground-motion hazard forecast that explicitly accounted for induced seismicity (Petersen et al., 2016). In regions affected by induced seismicity, the seismic hazard exceeds that associated with natural tectonic earthquakes alone. Figure 1 illustrates this effect by comparing the peak ground acceleration (PGA) hazard curves for Oklahoma City (35.50°, −97.55°) derived from the USGS 1-year induced seismicity forecast with the 2018 National Seismic Hazard Model (NSHM) that includes only tectonic sources (Petersen et al., 2018). Accounting for induced seismicity can substantially increase the estimated PGA hazard; consistent with these findings, Liu et al. (2019) projected a significant increase in seismic risk to buildings in central United States due to induced seismicity.

In engineering practice, seismic hazard and risk analyses commonly include damage thresholds for different types of infrastructure. Engineers typically define these thresholds using either earthquake magnitude (M) or ground-motion intensity measures (IMs), such as PGA and peak ground velocity (PGV). For example, the seismic performance assessments of buildings often entail a minimum magnitude cutoff in probabilistic seismic hazard analyses by assuming that smaller earthquakes generate insufficient energy to damage engineered structures. In seismic risk assessments, the fragility functions need to be dependent on magnitude to capture such behaviors. To avoid this added complexity, low-magnitude events are often excluded from hazard calculations by prescribing a minimum damaging magnitude. These magnitude thresholds generally rely on empirical observations or analyses using natural earthquake ground motions. However, the thresholds derived for tectonic earthquakes may not be applicable to induced seismicity owing to the greater amplitudes of short-period ground motions from smaller earthquakes at short distances observed in induced earthquakes. Consequently, small-magnitude induced earthquakes may pose a greater damage potential than their natural counterparts. This distinction underscores the need to explicitly evaluate the magnitude thresholds for damage potential in regions affected by induced seismicity.

Only a limited number of studies have addressed magnitude thresholds for damage from induced seismicity. For example, Baird et al. (2020) evaluated the damage thresholds for modern light-frame buildings using ground motions recorded from induced seismicity in central United States and western Canada; they concluded that earthquakes with magnitudes between 4.0 and 4.25 cause minimal damage, whereas those exceeding 4.5 may have high potential for damage. Similarly, Green et al. (2019) investigated the minimum magnitude of an induced earthquake required to trigger liquefaction and utilized field observations as well as parametric analyses to determine that earthquakes as small as magnitude 4.5 could induce liquefaction in susceptible soils. Although these studies highlight the importance of the magnitude threshold for induced seismicity, they focus exclusively on either light-frame building damage or liquefaction hazard. Because the minimum magnitude cutoffs used for hazard calculations serve as proxies for the magnitude dependence of fragility functions, the appropriate thresholds are expected to vary across engineering applications. For instance, a magnitude of 4.5 may be suitable for defining the seismic demands on light-frame buildings but inappropriate for assessing the seismic risk to earth canal dykes, which are prevalent in regions affected by induced seismicity and are the primary focus of this study.

In this study, we evaluated the minimum earthquake magnitude capable of damaging earth canal dykes in areas with induced seismicity. Accordingly, we considered dyke properties representative of the systems commonly encountered in practice and examined three earthquake scenarios corresponding to potential magnitude thresholds. To this end, we adopted a modified version of the ground motion model (ASK14) developed by Abrahamson et al. (2014) to select the motions representing induced earthquakes. We used selected time histories of the ground motions as inputs to analytical and numerical models to evaluate the dyke performances. We present the results in terms of the probability of exceeding the specified deformation threshold for each earthquake scenario. Herein, we first describe the modifications to the ASK14 model and the procedure used to select ground motions; then, we outline different methods used to estimate seismically induced deformations, including advanced dynamic analyses. Finally, we present the results of the dynamic analyses and propose minimum-magnitude damage thresholds for earth canal dykes.

Estimating IM (e.g., PGA) amplitudes using ground-motion models (GMMs) is a fundamental step in seismic hazard analysis. Most existing GMMs were developed for tectonic (natural) earthquakes and have been primarily calibrated using large-magnitude events at moderate-to-large distances. Consequently, these models do not extrapolate well to small-magnitude, short-distance, and shallow earthquakes typical of induced seismicity. In particular, such models systematically underestimate the short-period ground motions at distances less than 10 km for magnitudes below 5.0. This bias is attributable to the limited availability of near-source recordings from small-magnitude tectonic earthquakes as well as the historical assumption that such events pose negligible damage potential. However, recent studies have demonstrated that small-to-moderate induced earthquakes at shallow depths and short distances can contribute significantly to seismic hazard (Petersen et al., 2016; Petersen et al., 2018; Liu et al., 2019). These findings motivate modification of tectonic-based GMMs for application to induced seismicity. The guiding principle for such modifications is ensuring that the magnitude and distance scaling of a modified GMM for small-to-moderate earthquakes is consistent with the observed ground motions from shallow events while preserving the scaling at large magnitudes. In this study, we modified the ASK14 GMM reported by Abrahamson et al. (2014) and extended its applicability to earthquakes with magnitudes below 5.6 at rupture distances of less than 15 km. The modification of the GMM is discussed subsequently, and additional details are available in an earlier technical report (Abrahamson et al., 2020).

When applied to small-magnitude events, the primary source of inappropriate extrapolation in existing GMMs formulated for natural tectonic settings is the treatment of the finite-fault term; this term accounts for deviations from the point-source assumption by incorporating the rupture distance R r u p into the geometric spreading formulation. For large earthquakes, high-frequency ground motions originate from the slip of distributed regions of the rupture plane rather than from the point closest to the site alone. To capture this effect, the ASK14 model introduces a finite-fault term H(M) in the geometric spreading term ln R r u p 2 + H M 2 . Because the near-source recordings used to formulate most GMMs are predominantly from large-magnitude tectonic earthquakes, this term is well constrained for such events but poorly represents small-to-moderate earthquakes at short distances that are more representative of induced seismicity. Small-magnitude earthquakes have compact rupture areas, and the seismic slip controlling the ground motions is concentrated over a limited region; hence, the finite-fault effects should be weaker and explicitly dependent on magnitude. Although the ASK14 model includes a magnitude-dependent reduction of the finite-fault term for M < 5, this adjustment is insufficient to reproduce the short-distance scaling of small-magnitude events observed in regions with induced seismicity, such as Oklahoma (Rennolet et al., 2018). We modified the finite-fault term of the ASK14 model to explicitly include magnitude dependence based on rupture-dimension scaling, thereby reducing the term for M < 5.6 while preserving the original formulation for larger magnitudes. The revised finite-fault term is defined as shown in Equation 1: H M = 4.5 , M ≥ 5.6 e 1.5 + 1.15 M − 5.57 , M < 5.6 . (1)

This modification improves the physical basis of the finite-fault representation relative to the original ASK14 formulation. Applying the revised finite-fault term improves the model fit at short distances for earthquakes with magnitudes of 4.5–5.5. However, residual analyses reveal underprediction of the short-period ground motions for smaller events at short distances. Figure 2 illustrates this behavior through the within-event residuals of the pseudospectral acceleration (PSA) at T = 0.1 s that remain positive for the M3 and M4 bins at distances less than 10 km. To address this remaining bias, a second modification is introduced to steepen the short-distance attenuation slope for small magnitudes at rupture distances less than 15 km. The revised distance-scaling formulation is given by Equations 2a–c: f d i s t M , R r u p = c 1 + f R R r u p , M T m M , (2a)

with the magnitude taper given by T m M = 1 , M < 4 5.5 − M 1.5 , 4 ≤ M ≤ 5.5 0 , M > 5.5   , (2b) and f R R r u p , M = c 2 ⁡ ln max R e f f , 5 + 0.1 R c + 0.1 , R _ e f f < R c 0 , R _ e f f < R c . (2c) Here, c 1 and c 2 are period-dependent coefficients estimated from the data used to develop the ASK14 model. The taper function confines the distance-scaling adjustment to small-to-moderate magnitudes, with full effect for M ≤ 5.0 and no effect at M = 5.5, consistent with the magnitude-dependent geometric spreading observed in the NGA-West2 GMMs as well as the Mogul (Nevada) earthquake (Anderson et al., 2009). The cutoff distance R_c was set to 15 km based on inspection of the residual trends, while the term max R e f f , 5 models the flattening trend of the short-distance residuals at distances less than 5 km. Together, these two modifications enable the ASK14 model to reproduce the short-distance scaling observed for earthquakes with magnitudes between 3 and 5, which are representative of induced seismicity. This modified GMM predicts substantially higher short-period ground motions at distances less than 5 km (Figure 3). Figure 3 shows the comparison of the response spectra for M3, M4, and M5 earthquakes at a rupture distance of 3 km and a site condition with V_s30 = 400 m/s as computed using the original ASK14 model, modified ASK14 model, and the A15 GMM developed by Atkinson (2015) to explicitly account for the induced-seismicity effects. The modified ASK14 spectra closely match those from the A15 model and exceed the original ASK14 predictions, particularly at high frequencies. This agreement indicates that the modified ASK14 GMM captures the ground-motion amplitudes for small-magnitude, short-distance, and high-frequency conditions representative of induced-seismicity earthquakes more accurately.

Selecting input ground motion time histories is a critical step in dynamic analysis. We defined the median target response spectra for three scenarios using the modified ASK14 model with magnitudes of 3, 4, and 5 at a rupture distance of 3 km. These target spectra represent ground motions characteristic of induced seismicity. For each scenario, we selected ten pairs of two-component horizontal time histories from the NGA-West2 database (Bozorgnia et al., 2014). The selection criteria included event magnitudes within ±0.5 of the target magnitude, rupture distances less than 15 km, and RotD50 spectral shapes consistent with the target spectrum. We applied light spectral matching to the selected records by preserving the variability between the two horizontal components while matching their geometric mean to the target spectrum. This approach retains the peak-to-trough variability of the individual components, which is not captured by the RotD50-based standard deviation of the ASK14 GMM. For the M3 and M4 scenarios, we performed spectral matching over the range of 0.02–1.0 s, whereas the M5 scenario entailed an extended range of 0.02–2.0 s. To verify the suitability of the matched records for deformation analysis, we compared their Arias intensity (AI) values with the expected AI value for each scenario. We estimated the target AI using the conditional GMM for AI proposed by Abrahamson et al. (2016), which relates the AI to the PGA, PSA at 1.0 s, magnitude, and distance defined for the target spectrum. The AI values of the matched time histories agree well with the target AI estimates (Supplementary Appendix A), confirming their appropriateness for the deformation analysis. Figure 4 presents the responses of the spectrally matched M4 records as an example.

We scaled the input time histories using constant factors derived from sampled epsilon values (i.e., the number of standard deviations above or below the median) ranging from −3 to +3 in increments of 0.5 to represent the aleatory variability of the ground-motion amplitudes. We based the scaling standard deviation on the ASK14 GMM over the range of 0.075–0.3 s, which is representative of short periods that are of interest to induced-seismicity events. In this context, the epsilon values can also be considered as averaged values over the same period range. Table 1 summarizes the resulting logarithmic standard deviations of the spectral acceleration values ln S a from the ASK14 model for the M3, M4, and M5 scenarios. Because the present study focuses on evaluating the minimum-magnitude damage thresholds rather than the magnitude-dependent variability, we adopted the average standard deviation for the M4 scenario (0.87) to scale the time histories across all magnitudes; using this value, we computed the scale factors corresponding to the 13 sampled epsilon values listed in Table 2.

The probability of each epsilon value was computed according to Equation 3 by assuming a normal P ε i = Φ ε i + 0.25 − Φ ε i − 0.25 , (3) where Φ x is the cumulative standard normal distribution. The probability for each epsilon value is presented in Table 2.

We used the selected time histories as inputs to dynamic analyses that compute the deformations using various approaches, as described below. We performed these analyses for the M3, M4, and M5 scenarios (Supplementary Appendix A) by applying the scale factors listed in Table 2.

Figure 8 presents some typical results from the numerical analyses in terms of the horizontal and vertical displacements corresponding to the ground motions producing the maximum deformation for the K y = 0.10 g scenario. The numerical simulation results indicate that the horizontal displacements are primarily concentrated along the canal slopes, whereas the vertical displacements are localized near the crest. Similar displacement patterns are observed for the other ground motions.

We present the results of the dynamic analyses using both analytical and numerical models, expressed in terms of exceedance probability curves for specific deformation thresholds. We computed the conditional probability of exceeding a deformation level z given an earthquake scenario and magnitude at a rupture distance R as provided in Equation 6: P Def > z | M , R = 1 2 N T H ∑ j = 1 N S F ∑ i = 1 N T H H Def p o s M , S F j , T H i − z + H Def n e g M , S F j , T H i − z P S F j , (6) where H(x) denotes the Heaviside function, N S F = 13 is the number of scale factors, and N T H = 20 is the number of time histories. For the analytical models, we averaged the conditional exceedance probabilities over the four dyke configurations listed in Table 3.

Figures 9a–c show the exceedance probability curves obtained from the analytical models (TFM and stick–slip) for yield accelerations of K y = 0.05 g, 0.1 g, and 0.15 g. Figure 9d shows a comparison of these curves with those derived from the numerical dynamic analyses. The numerical analyses consistently predict larger displacements than the analytical models. This difference is because the analytical approaches primarily account for shear-induced deformation mechanisms, whereas the numerical simulations capture both deviatoric and volumetric deformation components (to a certain extent). Consequently, the numerical models provide a more complete representation of the physical processes governing seismically induced displacements. These trends are consistent with the findings of Macedo et al. (2022), who reported similar differences between analytical and numerical predictions for rockfill dams.

We next discuss the exceedance probabilities associated with a 10-cm displacement threshold that we adopted as the onset of potential damage for the dyke systems evaluated in this study. These discussions are focused on K y = 0.1   g as this value is representative of many dyke systems. Table 5 summarizes the exceedance probabilities for all K y values for both the analytical and numerical approaches. For the M3 scenario, the analytical models predict that the probability of exceeding 10 cm is below 0.001 for all K y values; the numerical simulations yield a slightly higher probability of approximately 1/800 for K y = 0.1   g , but this value remains small and indicates that M3 earthquakes are unlikely to contribute significantly to seismic risk. In contrast, the M4 scenario exhibits substantially higher exceedance probabilities. For K y = 0.1   g , the probability of exceeding 10 cm is approximately 0.07 based on the analytical models and 0.30 based on the numerical simulation. The probabilities increase further for lower K y values (Table 5). Overall, the numerical analyses predict exceedance probabilities that are four to eight times higher than those obtained with the analytical methods. This ratio can guide the approximation of results at other K y values. For example, for K y = 0.05   g , the analytical models yield an average exceedance probability of approximately 0.15; applying a representative ratio of 6, the corresponding numerical exceedance probability could approach 0.9. For the M5 scenario, the exceedance probabilities are high for all yield accelerations. For K y = 0.1   g , both the analytical and numerical models predict probabilities of approximately 0.8; the analytical exceedance probabilities for K y = 0.05   g and 0.15 g also exceed 0.65, as shown in Figure 9. In this case, the analytical and numerical predictions converge because all models estimate large displacements, resulting in similar values for the summation in Equation 6.

In summary, the M3 scenario is unlikely to impose significant seismic risk, whereas the M4 and M5 scenarios are both associated with substantial risk. Among the scenarios evaluated herein, the M4 case represents the lowest magnitude at which significant risk emerges. However, because the transition from negligible to significant risk likely occurs between the M3 and M4 scenarios, we recommend a practical lower bound of approximately 3.5, above which induced earthquakes may pose meaningful seismic risk to earth canal dykes.

The marginal hazard for a given scenario is defined as the product of the occurrence rate of the scenario and the conditional probability of exceedance given that the scenario occurs, as depicted in Equation 7: R a t e   M , R . P P S A > z | M , R . (7)

In regions of induced seismicity, the earthquake occurrence rates can be high. Consequently, relatively large epsilon values (e.g., 2–3) can control the marginal hazard even for small-to-moderate-magnitude events. Under these conditions, while the median ground motions (epsilon ≈ 0) from M3–M4 earthquakes are generally modest and often non-damaging, ground motions associated with large epsilon values could pose substantial hazard. For example, using the modified ASK14 GMM, the 2.5-epsilon PSA at T = 0.2 s for an M3.5 earthquake at a rupture distance of 3 km is approximately 1.0 g compared to a median PSA of approximately 0.1 g. Thus, although most M3.5 induced earthquakes are unlikely to cause damage, rare large-epsilon realizations associated with high activity rates could generate damaging ground motions.

In engineering practice, earthquakes with magnitudes below structure-dependent thresholds are often assumed to generate insufficient energy for damaging engineered systems. These damage thresholds are typically derived from observational data or dynamic analyses using ground motions consistent with GMMs developed for tectonic (natural) earthquakes. However, the thresholds established for tectonic events should not be applied directly to induced seismicity, which is characterized by relatively large short-period ground-motion amplitudes from small-magnitude earthquakes at short distances. This discrepancy is because tectonic-based GMMs do not adequately represent magnitude and distance scaling under conditions typical of induced seismicity, which limit their applicability in such settings.

As demonstrated in this study, modifying both the finite-fault term and short-distance scaling of tectonic-based GMMs provides a practical pathway for improving ground-motion predictions for small-magnitude induced earthquakes. Using a modified version of the ASK14 model, we defined the target response spectra for three scenarios corresponding to earthquakes of magnitudes 3, 4, and 5 at a rupture distance of 3 km. We used these spectra to select specific ground motions and conducted dynamic deformation analyses of dyke structures using both analytical and numerical approaches. The numerical simulations consistently yielded larger displacements than the analytical models, reflecting the limitations of analytical formulations in capturing the deformation mechanisms. Although the numerical models employed in this study capture the physics involved in the generation of seismic displacements better, a more rigorous representation requires advanced constitutive models, such as those formulated within the framework of critical-state soil mechanics.

The dynamic analysis results indicate that ground motions consistent with the M3 scenario are unlikely to produce significant damage in the evaluated dyke systems regardless of the strength level ( K y = 0.05 g, 0.10 g, and 0.15 g). In contrast, the M4 scenario shows a non-negligible probability of damaging deformation, particularly for dykes at low-to-moderate strengths ( K y = 0.05 g and 0.10 g). The ground motions associated with the M5 scenario are expected to cause damage across all strength levels considered in this study. Collectively, these results indicate that for a rupture distance of 3 km, the minimum earthquake magnitude capable of producing significant deformation (>10 cm) in low-to-typical-strength dyke systems is between the magnitudes of 3.0 and 4.0. Consequently, the commonly adopted minimum magnitude threshold of 5.0 used in seismic hazard applications cannot be applied universally in regions affected by induced seismicity. Based on the scenarios evaluated in this study, a practical lower threshold bound of approximately M3.5 is therefore suggested. However, we acknowledge that these findings are based on a limited set of scenarios, including the use of a single representative dyke configuration in the numerical analyses. Additional case history evidence and broader parametric investigations are required to characterize the seismic risk to earth canal dykes more fully and to refine the minimum magnitude thresholds identified herein. Future efforts in this direction would therefore include evaluating other dyke geometries and more extensive ground motions. Finally, it is important to recognize that minimum magnitude thresholds are conceptually linked to the use of magnitude-independent fragility functions in seismic risk assessments. Such thresholds are commonly introduced as simplifications to exclude small-magnitude events from hazard calculations, thereby avoiding explicit magnitude dependence in fragility formulations. Future studies should also critically examine the validity of this assumption and its consistency with fragility-based risk frameworks for small-magnitude and short-distance earthquakes.