A seismologically consistent surface rupture length model for unbounded and width-limited events

Abstract

A new surface-rupture-length (SRL) relationship as a function of magnitude $(M)$ , fault thickness, and fault dip angle is presented in this article. The objective of this study is to model the change in scaling between unbounded and width-limited ruptures. This is achieved through the use of seismological-theory-based relationships for the average displacement scaling and the aid of dynamic fault rupture simulations to constrain the rupture width scaling. The empirical data set used in the development of this relationship is composed of $122$ events ranging from $M$ 5 to $8.1$ and SRL 1.1 to 432 km. The dynamic rupture simulations data set includes $554$ events ranging from $M$ 4.9 to 8.2 and $SRL$ 1 to 655 km. For the average displacement $(\bar{D})$ scaling, three simple models and two composite models were evaluated. The simple average displacement models are a square root of the rupture area $(\sqrt{A})$ , a down-dip width $(W)$ , and a rupture length $(L)$ proportional model. The two composite models have a $\sqrt{A}$ scaling for unbounded ruptures and transition to $W$ and $L$ scaling for width-limited events, respectively. The empirical data favors a $\bar{D} ~ \sqrt{A}$ scaling for the entire regime (unbounded and width-limited ruptures) followed by a $\sqrt{A}$ scaling for unbounded that changes to $L$ scaling for width-limited ruptures. The selected models exhibit better predictive performance compared to linear $\log (SRL) ~ M$ type models, especially in the large magnitude range, which is dominated by width-limited events. A comparison with published SRL models shows consistent scaling for different fault types and tectonic environments.

Keywords

Probabilistic fault displacement hazard analysis surface fault rupture seismic hazard displacement scaling relations seismology empirical model

Introduction

Surface-rupture fault displacement hazard analyses, both probabilistic and deterministic, require an estimate of the surface rupture length (SRL) to compute the shape of the slip profile (Lavrentiadis and Abrahamson, 2019; Moss and Ross, 2011; Petersen et al., 2011; Youngs et al., 2003). Studies such as Wells and Coppersmith (1994), Wesnousky (2008), and Wells and Youngs (2015) proposed empirical models for the SRL scaling using observations from past earthquakes. These models are straightforward to develop and have good predictive performances within the range of data, but they may exhibit poor extrapolation for large events as the mechanisms that control the SRL scaling in large earthquakes, which are not well understood due to the scant empirical data. Another approach for developing scaling relationships is based on the use of theoretical considerations and constraints. Examples of such models are Leonard (2010, 2014), who developed sets of equations that describe the scaling between seismic moment, rupture area, length, width, and average displacement by imposing a self-consistent constraint. This type of theory-based model may have slightly larger aleatory variability at the center of the data compared to the purely empirical models, due to the incorporation of seismological constraints, but exhibit better extrapolation behavior.

The impact of the thickness of the seismogenic zone in the geometry of the rupture has been observed by previous studies, for instance, Hanks and Bakun (2002) proposed a magnitude break in their $\log (A) ~ M$ relationship at magnitude $M = 6.71$ , while Leonard (2014) proposed bi-linear and tri-linear relationships for $M ~ \log (L)$ and $M ~ \log (SRL)$ relationships, depending on the style of faulting and tectonic environment. However, to the authors’ knowledge, a relationship that explicitly considers the effect of the finite seismogenic zone thickness in the $\log (SRL) ~ M$ scaling does not exist. The recent development of extensively documented community-based fault displacement databases and advances in computer-based earthquake rupture simulation methods (both described below) motivated the development of the model presented herein, which explores this issue. The objective of this study is to develop a new SRL relationship that captures the changes in scaling between unbounded and width-limited ruptures. The impact of the change in scaling is expected to be more pronounced in active crustal regions (ACRs) with thin seismogenic zones, such as California, as opposed to stable continental regions with thicker seismogenic zones, such as Australia or the northeastern part of the North American continent.

A sketch illustrating the two conditions is presented in Figure 1. Assuming a similar static stress drop between the two regions, events of similar magnitude are expected to have similar rupture areas $(A)$ (Equation 1 from Brune, 1971). In addition, making the assumption that ruptures grow with similar aspect ratio scaling, small events that produce surface rupture, which do not reach the bottom of the seismogenic zone, are expected to have similar rupture geometries (rupture length and width), between the two regions. However; for moderate-to-large events, the thinner seismogenic zone is expected to impede the growth of the rupture width at shallower depths, requiring the rupture length to grow faster in order to accommodate the same rupture area, which will result in a magnitude break and steeper scaling for SRL.

\log (M_{o}) = 3 / 2 \log (A) + \log (Δ σ) - 0.387

(1)

In the following sections, the term “unbounded rupture” is used to classify the ruptures of the events not limited by the width of the seismogenic zone, and the term “width-limited ruptures” is used to classify the ruptures of the events that are constrained by the finite width of the seismogenic zone.

Data

Two sets of data were used in developing the proposed model: an empirical data set used for the SRL scaling (Subsection: Empirical Data) and a numerical physics-based simulation data set used to constrain the rupture width scaling (Subsection: SCEC Simulations). The empirical and numerical simulation data sets can be accessed from the repositories provided in the Subsection: Data and Resources.

Empirical data

The empirical data sets used in the model development include: (i) the Fault Displacement Hazard Initiative (FDHI) data set (Sarmiento et al., 2022), (ii) additional events from Wells and Coppersmith (1994), which they classified as reliable, and (iii) the surface rupture events from Baize et al. (2020) which were not part of the FDHI data set. Events less than $M 5$ were excluded from the final data set, as well as the 1892 Laguna Salada, 1978 IzuOshima, and 1993 Killari earthquakes which were considered outliers. Potential reasons for considering the previous events as outliers are the Laguna Saldana earthquake ruptured in the 1800s but was mapped in the 2010s, so large parts of the rupture may have been irrecoverable; the 1978 Izu Oshima earthquake likely included an offshore segment that was not mapped in the FDHI data set; and the 1993 Killari occurred in a stable continental region with a lot of deformation probably accommodated via folding/warping that was unmapped.

The SRL values for the events from Wells and Coppersmith (1994) and Baize et al. (2020) were obtained directly from the corresponding studies, while for the events in the FDHI data set, SRL was estimated based on the length of the event coordinate system (ECS) defined in Lavrentiadis and Abrahamson (this issue). ECS’s goal is to define the along-strike and perpendicular-to-strike distance metrics for the FDHI events whose points are georeferenced on the latitude and longitude coordinates. Comparisons of the SLR reported by Wells and Coppersmith (1994) and Baize et al. (2020) with those obtained from ECS showed unbiased consistent estimates. Figure 2 shows the SRL versus $M$ distribution of the final empirical data set. It is composed of 64 events from FDHI, 52 events from Wells and Coppersmith (1994), and 6 events from Baize et al. (2020). The earthquake $M$ ranges from 5.0 to 8.1, while SRL ranges from 1.1 to 432 km. It contains 64 strike-slip, 26 normal, and 32 reverse earthquakes.

The fault width $(W_{\lim})$ required to define the transition to a width-limited rupture was estimated with Equation 2 based on the thickness of the seismogenic zone $(D_{seis})$ and the fault dip angle $(θ)$ . $D_{seis}$ was based on the fault location, as summarized in Table 1, whereas $θ$ was based on the style of faulting, as summarized in Table 2. In the absence of fault-specific information, representative $D_{seis}$ and $θ$ values were obtained from Huang et al. (2023), while the number of events belonging to each region is also summarized in Table 1.

W_{\lim} = \frac{D_{seis}}{\sin (θ)}

(2)

Table 1.

Seismogenic thickness and number of events of different regions

Region	Seismogenic thickness (km)	Number of events
California	15	29
Guatemala	15	1
New Zealand	15	3
Indonesia	15	1
Japan	17	7
Himalayan region	30	12
Australia	40	11
Other regions	22	57

Table 2.

Dip angle values for different styles of faulting

Style of faulting	Dip angle (deg)
Strike-slip	90
Normal	60
Reverse	45

Dynamic rupture simulations

The data set used for constraining the rupture width scaling is based on physics-based fault displacement simulations conducted by the Southern California Earthquake Center (SCEC) as part of the FDHI efforts. The physics-based approach used in this study is referred to as the dynamic rupture model (Harris et al., 2018); it constructs spontaneously evolving earthquake ruptures under mechanical causative conditions (e.g., fault geometry, friction laws, stress conditions, and surrounding rock properties). Figure 3 illustrates a simplified workflow of the dynamic rupture model. To simulate realistic fault displacement, we follow the methodology (Wang and Goulet, 2021) originated from the 1992 Landers earthquake. We adopt a depth-dependent lithostatic load for vertical stress and self-similarly heterogeneous shear stress that permits unprescribed variability in rupture style and size. With regard to the friction law, we similarly select the slip-weakening law (Ida, 1972). In contrast to the general theoretical considerations such as those used by Leonard (2010), the dynamic rupture model employs a physically plausible parametric uncertainty (i.e., heterogeneous initial stress) to capture the displacement variability across the rupture plane due to the elasto-plasto-dynamic response to imposed stresses. The SRL is then computed using numerical criteria. For the purpose of this article, the length is represented by summing the length of fault segments with surface displacement larger than 1 cm.

The fault displacements were simulated by numerically solving the 3D elastoplastic spontaneous rupture propagation with the Support Operator Rupture Dynamics (SORD) code (Wang and Day, 2017, 2020; Wang and Goulet, 2021). This application is highly optimized and scalable on current cutting-edge supercomputers. This is a requirement since to capture a $M$ range of 5–8 associated with rupture lengths of several kilometers to hundreds of kilometers, the computational demands vary from seconds on tens of CPUs to hours on hundreds of thousands of CPUs. The simulations described here were performed on Theta at the Argonne Leadership Computing Facility and Frontera at the Texas Advanced Computing Center.

A series of dynamic fault rupture simulations under various model setups were therefore performed to understand and model the effect of fault width in the SRL scaling. The general parameters were set based on the work of Wang and Goulet (2021), with additional variations, as described below, all for vertical strike-slip faults. One set of simulations was performed in which the fault width was limited to 19 km, with two variations in the normal loading stresses pattern: fixed with depth in one case, and linearly increasing with depth in the other. The other set of simulations involved limiting the fault width to 15 km while maintaining the other model parameters fixed. It is noted that although a width limit is provided for the dynamic rupture models (15 and 19 km), the resultant ruptures, which are generated based on the initial conditions and physical laws, are not prescribed, which results in the scattered characteristics shown in Figures 4 and 5. The simulations output includes determining the moment mangitude ( $M$ ), which is not a priori imposed, the rupture geometry, and the displacement field across the fault plane. The simulations were continued to generate a wavefield and to provide surface ground motions, which were verified against existing ground motion models (GMMs), to ensure that the simulations were reasonable and technically defensible. A detailed description of the dynamic rupture model can be found in Wang and Goulet (2022).

In total, 554 simulations were performed with magnitude ranging from 4.9 to 8.2 and SRL ranging from 1 to 655 km. Figure 4 compares the $M - SRL$ scaling of the empirical and SCEC data sets. Overall, the empirical and simulation data sets appear consistent, supporting the use of the simulation data sets to inform the mean SRL scaling in the empirical data. The aleatory variability of the SCEC data set is expected to be less than that of the empirical data, as the goal of the simulations was to capture average scaling effects and not to fully represent the randomness and uncertainty in earthquake ruptures.

Model development

The derivation of the SRL model is divided into two parts. In the first part, an average model for the rupture width $(W)$ scaling of unbounded and width-limited ruptures based on the SCEC simulations is proposed. In the second part, the SRL model based on the empirical data is developed with the aid of the $W$ model to capture the transition between unbounded and width-limited ruptures.

Rupture width modeling

Although the physics-based simulations have complete slip profiles on fault planes, the rupture widths are sensitive to where rupture arrest may be subjectively measured. To avoid this, we adopt the approach originally developed for extracting rupture widths from finite fault models by Mai and Beroza (2000) to define the effective rupture width (the autocorrelation width of depth-variable slip profile) and apply this method to all simulation-based data sets for rupture width measurements. Then, the effective rupture widths, estimated from the SCEC simulations, were used to determine the width scaling relationship. A linear functional form with a plateau for the upper limit is proposed:

\log_{10} (W) = min (β_{1} + β_{2} M + ϵ_{W}, \log_{10} (W_{\lim}))

(3)

in which the log of the rupture width scales linearly with magnitude until it reaches the fault width, where it remains constant. The coefficient $β_{1}$ is the model intercept, $β_{2}$ controls the magnitude scaling, and $ϵ_{W}$ is the aleatory term that is modeled with a normal distribution with a zero mean and $σ_{W}$ standard deviation $(ϵ_{W} ~ N (0, σ_{W}))$ The model coefficients were estimated with a linear regression; the best estimates and standard error values of the model coefficients are provided in Table 3. In the SCEC simulations, events were classified as width-limited if the effective width was within ±2.5 km of the prescribed maximum width. Figure 5 compares the scaling of the proposed model with the SCEC simulations, and Figure 6 presents the regression residuals, where it can be seen that there is a small misfit in the different magnitude bins, but overall, the model captures the fundamental scaling of the SCEC simulations and the transition between unbounded and width-limited ruptures. A tri-linear relationship is expected to improve the fit; however, it will also introduce additional complexity, which cannot be adequately characterized with the limited number of simulations; for this reason, the simpler bi-linear model was selected. In addition, any simplifications or differences in the width scaling between the SCEC simulations and empirical data will be mapped into the aleatory variability of the SRL model described in the next section.

Table 3.

Rupture width model coefficients

Coefficient	$β_{1}$	$β_{2}$	$σ_{W}$
Estimate	−1.1602	0.3395	0.0544
Standard error	0.0325	0.0053	−

SRL modeling

In developing the SRL model, the candidate functional forms are derived in Section: Functional Form Derivation, and the model coefficients are estimated in Section: Functional Form Derivation along with a discussion on the selection of the preferred and alternative models.

Functional form derivation

Starting with the definition of the seismic moment $(M_{0} = μ \bar{D} L W)$ , and substituting it into the definition of the moment magnitude ( $M = 2 / 3 \log_{10} (M_{0}) - 10.7$ ), (Kanamori, 1977), a relationship between $M$ , average rupture slip $(\bar{D})$ , subsurface rupture length $(L)$ , subsurface rupture width $(W)$ , and crust stiffness $(μ)$ is obtained:

M = \frac{2}{3} \log_{10} (μ) + \frac{2}{3} \log_{10} (\bar{D}) + \frac{2}{3} \log_{10} (L) + \frac{2}{3} \log_{10} (W) - 10.7

(4)

Treating $μ$ as model constant and combining it with the $- 10.7$ factor, the previous equation reduces to:

M = \frac{2}{3} \log_{10} (\bar{D}) + \frac{2}{3} \log_{10} (L) + \frac{2}{3} \log_{10} (W) + c

(5)

Here and in the remainder of this article, $c$ with no subscript is used to indicate any arbitrary constant, not a specific model coefficient.

Two end-member scaling relationships for average displacement often discussed in the literature (Bodin and Brune, 1996; Hanks and Bakun, 2002; Leonard, 2010; Pegler and Das, 1996; Romanowicz, 1992, 1994; Romanowicz and Rundle, 1993; Scholz, 1982, 1994, 1997, 1998) are the L model, in which average displacement is proportional to rupture length $(\bar{D} ~ L)$ , and the W model, in which average displacement is proportional to rupture width $(\bar{D} ~ W)$ . More recent studies suggest that slip scaling is between L and W models (Bodin and Brune, 1996; Mai and Beroza, 2000). Leonard (2010) proposed the $\bar{D} ~ \sqrt{A}$ model, in which the average displacement scales proportionally to the square root of the rupture area. Here, we investigated all three scaling relationships in building the $SRL ~ M$ model, as well as two composite scaling laws that follow $\bar{D} ~ \sqrt{A}$ for the unbounded ruptures and transition to $\bar{D} ~ L$ and $\bar{D} ~ W$ for width-limited events, respectively. In particular, the $SRL$ relationship based on $\bar{D} ~ \sqrt{A}$ scaling is denoted model 1, the relationship based on $\bar{D} ~ W$ scaling is denoted model 2, and the relationship based on $\bar{D} ~ L$ scaling is denoted model 3. The scaling relationship that starts with $\bar{D} ~ \sqrt{A}$ and transitions to $\bar{D} ~ W$ is denoted model 2′ while the scaling relationship that transitions to $\bar{D} ~ L$ is denoted model 3′.

By substituting the $\bar{D} ~ \sqrt{A}$ scaling in Equation 5, the following scaling relationship between $M$ , $L$ , and $W$ is obtained for model 1:

M = \log_{10} (L) + \log_{10} (W) + c

(6)

Combining Equation 5 with the $\bar{D} ~ W$ scaling, the model 2 scaling relationship becomes:

M = \frac{2}{3} \log_{10} (L) + \frac{4}{3} \log_{10} (W) + c

(7)

While, combining Equation 5 with the $\bar{D} ~ L$ model, the model 3 scaling relationship becomes:

M = \frac{4}{3} \log_{10} (L) + \frac{2}{3} \log_{10} (W) + c

(8)

Adopting the rupture width scaling relationship from Section: Rupture Width Modeling assuming the surface to subsurface rupture length relationship is of the form: $lo g_{10} (SRL) = c_{1} \log_{10} (L) + c$ (Leonard, 2010), and solving for $lo g_{10} (SRL)$ , model 1 transforms to:

\log_{10} (SRL) = c_{1} (M - min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c

(9)

model 2 transforms to:

\log_{10} (SRL) = c_{1} (\frac{3}{2} M - 2 min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c

(10)

and model 3 transforms to:

\log_{10} (SRL) = c_{1} (\frac{3}{4} M - \frac{1}{2} min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c

(11)

Considering Equations 9 and 10, the functional form the composite model 2′ is:

\begin{matrix} \log_{10} (SRL) = {\begin{matrix} c_{1} (M - min (β_{1} + β_{2} M, \log_{10} (W_{\lim})) + d_{1}) + c & M \leq M_{\lim} \\ c_{1} (\frac{3}{2} M - 2 min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c & M > M_{\lim} \end{matrix} \end{matrix}

(12)

where $M_{\lim}$ is the magnitude at which the transition from unbounded to width-limited events occurs. The term $d_{1}$ is introduced to ensure a continuous transition between the two branches. Using Equation 3, $M_{\lim}$ is calculated as: $M_{\lim} = 1 / β_{2} (\log_{10} (W_{\lim}) - β_{1})$ . The transition from unbounded to width-limited events occurs at the $(M_{\lim}, W_{\lim})$ coordinate pair. Setting the two branches equal to each other at the transition point and solving for $d_{1}$ , we obtain:

d_{1} = - \frac{1}{2} M_{\lim} + \log_{10} (W_{\lim})

(13)

With the previous equation, the two branches of model 2′ are combined into one equation:

\log_{10} (SRL) = c_{1} (\frac{3}{2} M - \frac{1}{2} min (M, M_{\lim}) - min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c

(14)

Similarly, considering Equations 9 and 11, the functional form of the composite model 3′ is:

\begin{matrix} \log_{10} (SRL) = {\begin{matrix} c_{1} (M - min (β_{1} + β_{2} M, \log_{10} (W_{\lim})) + d_{2}) + c & M \leq M_{\lim} \\ c_{1} (\frac{3}{4} M - \frac{1}{2} min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c & M > M_{\lim} \end{matrix} \end{matrix}

(15)

Following a similar derivation to model 2′, the functional form of model 3′ can be expressed as:

\log_{10} (SRL) = c_{1} (\frac{3}{4} M + \frac{1}{4} min (M, M_{\lim}) - min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c

(16)

The style of faulting is implicitly included in the previous scaling relationships through $W_{\lim}$ ; normal and reverse faults have shallower dip angles compared to strike-slip faults leading to larger fault widths for the same seismological zone thickness. However, an additive term for dip-slip faults $(F_{D})$ is also added in the final functional forms to account for differences in the static stress drop between the different styles of faulting. In preliminary regressions, separate additive terms for normal and reverse faults were also evaluated but were rejected to favor model simplicity as both reverse and normal faults had similar additive term values. The static stress drop is defined as:

Δ σ = C μ \frac{\bar{D}}{L_{C}}

(17)

where $L_{C}$ is the characteristic length, and $C$ is a constant which is a function of the rupture geometry. In the log scale, a shift in the static stress drop will be accompanied by a shift in the average displacement $(\log (\bar{D}) = \log (Δ σ) + c)$ , which, in turn, will shift the intercept of the $SRL ~ M$ relationship. Thus, the final functional form for model 1 becomes:

\begin{matrix} \log_{10} (SRL) = c_{0} + c_{1} (M - min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c_{2} F_{D} + ϵ_{SRL} \\ = c_{0} + c_{1} X_{1} + c_{2} F_{D} + ϵ_{SRL} \end{matrix}

(18)

which can be expressed as a linear model with $X_{1} = (M - min (β_{1} + β_{2} M, \log_{10} (W_{\lim})))$ . The functional form for model model 2′ is:

\begin{matrix} \log_{10} (SRL) = c_{0} + c_{1} (\frac{3}{2} M - 2 min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c_{2} F_{D} + ϵ_{SRL} \\ = c_{0} + c_{1} X_{2} + c_{2} F_{D} + c_{3} + ϵ_{SRL} \end{matrix}

(19)

with $X_{2} = (3 / 2 M - 2 min (β_{1} + β_{2} M, \log_{10} (W_{\lim})))$ to be expressed as a linear model. The functional form for model 3 is:

\begin{matrix} \log_{10} (SRL) = c_{0} + c_{1} (\frac{3}{4} M - \frac{1}{2} min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) + c_{2} F_{D} + ϵ_{SRL} \\ = c_{0} + c_{1} X_{3} + c_{2} F_{D} + c_{3} + ϵ_{SRL} \end{matrix}

(20)

with $X_{3} = (3 / 4 M - 1 / 2 min (β_{1} + β_{2} M, \log_{10} (W_{\lim})))$ . Model 2′ final functional form is:

\begin{matrix} \log_{10} (SRL) = c_{0} + c_{1} (\frac{3}{2} M - \frac{1}{2} min (M, M_{\lim}) - min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) \\ + c_{2} F_{D} + ϵ_{SRL} \\ = c_{0} + c_{1} {X'}_{2} + c_{2} F_{D} + c_{3} + ϵ_{SRL} \end{matrix}

(21)

with $X'_{2} = (\frac{3}{2} M - \frac{1}{2} min (M, M_{\lim}) - min (β_{1} + β_{2} M, \log_{10} (W_{\lim})))$ , and model 3′ final functional form is:

\begin{matrix} \log_{10} (SRL) = c_{0} + c_{1} (\frac{3}{4} M + \frac{1}{4} min (M, M_{\lim}) - min (β_{1} + β_{2} M, \log_{10} (W_{\lim}))) \\ + c_{2} F_{D} + ϵ_{SRL} \\ c_{0} + c_{1} {X'}_{3} + c_{2} F_{D} + c_{3} + ϵ_{SRL} \end{matrix}

(22)

with $X'_{3} = (\frac{3}{4} M + \frac{1}{4} min (M, M_{\lim}) - min (β_{1} + β_{2} M, \log_{10} (W_{\lim})))$ .

Finally, a Wells and Coppersmith type model was also evaluated, hereafter referred to as model 0, to investigate the impact of width-limited ruptures on the current state of practice models. The functional form for model 0 is:

\log_{10} (SRL) = c_{0} + c_{1} M + c_{2} F_{D} + ϵ_{SRL}

(23)

In all previous equations, $c_{0}$ is the model intercept, $c_{1}$ controls the magnitude scaling, $c_{2}$ captures the median shift between strike-slip and dip-slip events ( $F_{D}$ is zero for strike-slip and one for reverse and normal faults), and $ϵ_{SRL}$ is the aleatory term. More information on the modeling of the aleatory variability is provided at the end of Subsection: Model Regression.

Model regression

All candidate models were estimated using a maximum likelihood linear regression and the empirical data set. Table 4 provides the log-likelihood $(L)$ and Akaike information criterion (AIC, Akaike, 1998) for the different models; a higher $L$ and lower AIC indicate a better model fit to the data. Table 5 summarizes the best estimates and standard errors of the model coefficients.

Table 4.

Log-likelihood $(L)$ and Akaike information criterion (AIC) of the candidate rupture length models

Model	$L$	AIC
Model 0	6.06	−4.12
Model 1	15.42	−22.85
Model 2	14.15	−20.30
Model 3	14.30	−20.60
Model 2′	8.91	−9.81
Model 3′	12.22	−16.44

Table 5.

Rupture length candidate model coefficients

Model	Coefficient	$c_{0}$	$c_{1}$	$c_{2}$	$σ_{SRL - unb}$
Model 0	Estimate	2.9680	0.6629	−0.1836	0.2331
Model 0	Standard error	0.2126	0.0305	0.0484	−
Model 1	Estimate	−3.8783	0.9509	−0.1517	0.2494
Model 1	Standard error	0.2191	0.0363	0.0457	−
Model 2	Estimate	−4.1798	0.7152	−0.1390	0.2520
Model 2	Standard error	0.2333	0.0276	0.0463	−
Model 3	Estimate	−3.6049	1.1316	−0.1645	0.2517
Model 3	Standard error	0.2110	0.0436	0.0459	−
Model 2′	Estimate	−3.0977	0.8044	−0.1291	0.2631
Model 2′	Standard error	0.2019	0.0327	0.0486	−
Model 3′	Estimate	−4.2507	1.0235	−0.1759	0.2560
Model 3′	Standard error	0.2406	0.0403	0.0466	−

From a statistical perspective, candidate models 1 to 3 have similar good performance, with model 1 having the best performance. Models 3′ and model 2′ which have poorer performance, while model 0 has the worst performance showcasing the limitations of a purely empirical model. Models 1 to 3′ have a break in the magnitude scaling when the width of the rupture reaches the width of the seismogenic zone, allowing them to fit the data better.

However, based on seismological theory, models 1 to 3′ should also have unit slopes with respect to the magnitude scaling terms (e.g., $X_{1}$ , $X_{2}$ , $X_{3}$ ), if the average displacement in the empirical data follows the assumed scaling law and the subsurface rupture length $(L)$ was the response variable. With SRL being the response variable, a magnitude scaling term greater than unity implies that the surface-to-subsurface rupture length ratio increases with magnitude, while a magnitude scaling term less than unity implies the opposite. Intuitively, a larger than one magnitude scaling term is expected, as for large events, a bigger part of the rupture is expected to reach the surface. Considering that model 1 exhibited the best predictive performance, but the magnitude scaling term was slightly less than one, the final preferred model uses model 1 functional form (Equation 18) with a fixed unit slope. Similarly, model 3′ has a slope consistent with seismological theory and decent predictive performance; thus, model 3′ with a fixed unit slope (Equation 21) is selected as an alternative physically plausible model. The statistical measures for the goodness of fit and coefficients of the preferred and alternative models are summarized in Tables 6 and 7, respectively.

Table 6.

Log-likelihood $(L)$ and Akaike information criterion (AIC) of the preferred and alternative models

Model	$L$	AIC
Model 1, fixed slope (preferred model)	14.49	−22.99
Model 3′, fixed slope (alternative model)	12.05	−18.10

Table 7.

Rupture length preferred and alternative model coefficients

Model	Coefficient	$c_{0}$	$c_{1}$	$c_{2}$	$σ_{SRL - unb}$
Model 1, fixed slope	Estimate	−4.1728	1	−0.1287	0.2503
(preferred model)	Standard error	0.0242	−	0.0425	−
Model 3′, fixed slope	Estimate	−4.1113	1	−0.1855	0.2553
(alternative model)	Standard error	0.0248	−	0.0434	−

The findings of the regression analyses favor a $~ \sqrt{A}$ scaling for the entire regime or a $~ \sqrt{A}$ scaling that transitions to a $~ \sqrt{L}$ scaling for width-limited events for the average displacement compared to a pure $~ W$ and $~ L$ scaling or a $~ \sqrt{A}$ scaling that transitions to a $~ \sqrt{W}$ scaling for width-limited events. The SCEC simulations show some extra support for $~ \sqrt{A}$ scaling. While a limited width (15 and 19 km) is set for all simulations, the smooth transition near the limited width and deeper scattering width for larger earthquakes indicate that a fixed limited width is not suitable for the general scaling. As shown in simulations of Wang and Goulet (2022), the deep penetration of large earthquakes keeps increasing average displacements with the rupture length approaching a $\sqrt{A}$ scaling. In addition, Thingbaijam et al. (2017) observed $D ~ A^{0.429 \pm 0.134}$ for reverse events, $D ~ A^{0.858 \pm 0.214}$ for normal events, and $D ~ A^{0.597 \pm 0.112}$ for strike-slip events, which also support a $\sqrt{A}$ scaling for reverse and strike-slip events, while for normal events a relationship that is in between linear and square-root scaling. In hazard calculations, the square-root dependence has been applied as one of the logic stress dealing with epistemic uncertainties in the Uniform California Earthquake Rupture Forecast version 3 (Field et al., 2015).

Figure 7 presents the median magnitude scaling for the preferred and alternative models for strike-slip events for a seismological zone thickness equal to 15 and 20 km against the empirical data for the same fault type and seismological zone thickness up to 25 km. For the 15 km thick zone, the magnitude break occurs at $M 6.9$ . For the 20 km crust, a larger rupture area is needed for the event to become width-limited; thus, the magnitude break occurs closer to $M 7.3$ . The change in magnitude scaling is more significant in the preferred model than in the alternative model. The empirical data show consistent trends; at small magnitudes, there is no systematic difference between events in thin and thick parts of the seismogenic zones, while at larger magnitudes, there is a positive shift for events in thin parts of the seismogenic zone compared to events on thicker seismogenic zones.

Figure 8 shows the preferred and alternative model scaling for strike-slip and reverse-slip events for a seismological zone thickness equal to 15 km. The vertical offset between the two models is introduced by the $F_{D}$ term, and the difference in the magnitude break is caused by the different dip angles between strike-slip and reverse faults. Although the same seismological thickness was assumed for both fault styles, the shallower dip angle of reverse faults allows for a wider fault width which can accommodate larger events before they become width-limited. The preferred model captures the style of faulting effects primarily through the break in magnitude scaling suggesting that they are more noticeable at large magnitudes, while the alternative model captures the style of faulting effects primarily through the additive term suggesting they are similar for small-to-moderate and large events. The magnitude break for the strike-slip events occurs at $M 6.9$ and for reverse events at $M 7.3$ . Consistent trends are also observed in the empirical data.

Figure 9 compares the regression residuals from the W&C type model (model 0) with the residuals from the preferred and alternative models (models 1 and 3′ with fixed slopes). All models have similar zero-mean-centered residuals for unbounded ruptures. However, Figure 9a shows a positive bias in the W&C type model residuals for width-limited ruptures, whereas in Figure 9b and c the residuals for width-limited ruptures of the selected models are centered closer to zero. Furthermore, the width-limited residuals exhibit less bias in the preferred compared to the alternative model. This comparison illustrates the advantage of seismological-theory-based models; a quadratic functional form would address the positive bias for the residuals but it would provide no basis for its existence. Relating the change in magnitude scaling to the finite thickness of the seismogenic zone increases the confidence in an extrapolation behavior that is scientifically defensible.

Preliminary regressions showed a magnitude-dependent aleatory variability that is by a factor of two smaller for width-limited ruptures as compared to unbounded ruptures. A rupture width transition parameter $(δ_{W})$ is introduced, and defined as:

δ_{W} = \log_{10} (W_{unb}) - \log_{10} (W_{\lim}) = β_{1} M + β_{2} - \log_{10} (W_{\lim})

(24)

where $W_{unb}$ is the theoretical rupture width for an infinite thickness seismological zone estimated from the first leg of rupture width model from Equation 3 $(\log_{10} (W_{unb}) = β_{1} M + β_{2})$ with the estimated coefficients in Table 3. A negative $δ_{W}$ indicates an unbounded rupture, while a positive $δ_{W}$ indicates a width-limited rupture, with the transition occurring around zero. Using the parameter $δ_{W}$ , a heteroscedastic standard deviation model for the aleatory variability is proposed $(ϵ_{SRL} ~ N (0, σ_{SRL} (δ_{W})))$ that gradually shifts from unbounded to width-limited ruptures based on a sigmoid function:

σ_{SRL} (δ_{W}) = \frac{σ_{SRL - unb}}{1 + S (10 \ln (9) δ_{W})}

(25)

where $σ_{SRL}$ is the magnitude-dependent standard deviation, and $σ_{SRL - unb}$ is the standard deviation for unbounded ruptures reported in Table 7. The sigmoid function is defined as $(S (x) = 1 / (1 + \exp (- x)))$ . Figure 10 presents the regression residuals and proposed standard deviation versus the parameter $δ_{W}$ for the preferred and alternative models. The factor $10 \times \ln (9)$ is used so that $80 %$ of the standard deviation change occurs between $δ_{W} = - 0.1$ and $0.1$ , which is consistent with the empirical observations.

Comparison with existing relationships

A series of comparisons with available models was performed to ensure the reasonableness of the proposed relationships. Figure 11 compares the preferred and alternative models with existing SRL relationships. “WC94 All,”“WC94 SS,” and “WC94 N” correspond to the Wells and Coppersmith (1994) $SRL ~ M$ scaling relationships for all, strike-slip, and normal events; the scaling relationship for reverse events is omitted due to low significance levels. “W08 All,”“W08 SS,” and W08 R” are the Wesnousky (2008) SRL relationships for all, strike-slip, and reverse shallow-crustal events. The Wesnousky (2008) SRL relationship for normal events is not included in this comparison due to poor constraints. “L10 SS” and “L10 DS” correspond to the Leonard (2010) $SRL ~ M$ relationships for dip-slip and strike-slip shallow-crustal events, while “L10 SCR” corresponds to stable continental regions. “L14 SCR DS” is the updated Leonard (2014) relationship for dip-slip events in stable continental regions. Finally, “WY15 OLS” and “WY15 Deming” correspond to the Wells and Youngs (2015) $SRL ~ M$ relationships estimated with ordinary least squares and errors in variables regression methods.

The developed models are in overall agreement with the existing relationships over different magnitude ranges (Figure 11). In Figure 11a, the existing models for all and strike-slip events cover the range of both the preferred and alternative scaling relationships. The relationships from Wells and Coppersmith (1994) and Leonard (2010) are in better agreement with the proposed model at the small-to-moderate magnitude range. In contrast, Wesnousky (2008) and Wells and Youngs (2015) are in better agreement with the proposed relationships at the large magnitude range. For normal and reverse events in shallow-crustal environments (Figure 11b and c), the existing relationships overestimate the size of SRL for small-to-moderate magnitudes, while the proposed relationships predict larger SRL for large events. Leonard (2010) has a similar slope to the developed models, whereas Wells and Coppersmith (1994) and Wesnousky (2008) have a lower slope compared to the developed models. Finally, Figure 11d compares existing models for stable continental regions with the predictions of the proposed models, which were produced for dip-slip events, at a 30 km thick crust and a 60° fault angle. The proposed models have similar slopes to the existing relationships. Due to the increased thickness of the crust and shallow dip angle, the proposed relationships are predominately linear, with the magnitude break only occurring at very large magnitudes $(M > 8)$ .

Range of applicability

A $M 5$ lower limit is recommended as it is the smallest magnitude in the regression data set. For the upper limit, the largest event in the data set was $M 8.1$ . However, the authors believe that seismological constraints will help the extrapolation up to $M 8.5$ , albeit with increased epistemic uncertainty.

Discussion and conclusions

The proposed SRL models capture the change in magnitude scaling between unbounded and width-limited ruptures through the use of seismological constraints and dynamic fault rupture simulations. Seismological theory was used to derive candidate scaling relationships between the moment magnitude $(M)$ , subsurface length $(L)$ , and rupture width $(W)$ . The dynamic rupture simulations constrained the $W ~ M$ scaling. The empirical data were used to decide between the alternative scaling relationships for average displacement, to model the difference between $L$ and SRL scaling, and to capture the aleatory variability in SRL scaling. The empirical data supports a $\sqrt{A}$ -type scaling for the entire regime followed by the alternative model, which has a $\sqrt{A}$ -type scaling for unbounded ruptures and transitions to $L$ -type scaling for width-limited events. Other scaling relationships that were evaluated were pure $W$ -type and $L$ -type scaling, and composite relationship with $\sqrt{A}$ scaling for unbounded events that transition to $W$ scaling for width-limited ruptures.

Compared to a simple linear regression between SRL and $M$ , the proposed functional form provides a better fit to empirical data sets. Both the linear and the proposed functional forms have similar fits for unbounded ruptures, but the proposed relationship provides a better fit for width-limited events. A quadratic functional form could have been used to resolve the positive bias at large magnitudes, but it would lack a seismological basis limiting the confidence for its extrapolation. A width transition parameter, $δ_{W}$ , is proposed to capture the difference in aleatory variability between unbounded and width-limited ruptures. The parameter $δ_{W}$ corresponds to the log of the ratio of the theoretical rupture width for an infinitely thick seismogenic and the actual fault width. A comparison with existing SRL relationships showed consistent scaling with the proposed models for a range of fault types and tectonic environments.

Future studies should evaluate the effect of the thickness of the seismogenic zone using fault-specific information. In addition, a fully consistent scaling model for magnitude, rupture area, subsurface rupture length and width, and SRL for unbounded and width-limited ruptures using a more comprehensive data set following the presented framework is encouraged.

Footnotes

Acknowledgements

The authors thank Stephane Baize for providing the Baize et al. data set and Alexandra Sarmiento for the insightful discussions during the development of this model. We would like to also thank the two anonymous reviewers and associate editor for the review and constructive comments that helped to improve the final article.

Correction (April 2025):

Article updated to correct Article type.

Authors’s note

Yongfei Wang is now affiliated to Verisk Extreme Event Solution, Boston, MA, USA.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the California Energy Commission, California Department of Transportation, Pacific Gas and Electric Company, and Southern California Earthquake Center (SCEC), SCEC Contribution #12712. SCEC is funded by the National Science Foundation (NSF) and the U.S. Geological Survey (USGS) through cooperative agreements with the University of Southern California (USC). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the sponsoring agencies. An award of computer time was provided by the INCITE program and this research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. This research was also partially computed on Frontera computing project at the Texas Advanced Computing Center through an allocation made possible by the National Science Foundation award OAC-1818253.

ORCID iDs

Grigorios Lavrentiadis

Yousef Bozorgnia

Christine A Goulet

Data and resources

The regression code and regression data sets are provided in: https://github.com/NHR3-UCLA/LWABC23_SRL_model. The statistical regressions were performed using the computer software R and stats package (R Core Team, 2022). The open-source software package Support Operator Rupture Dynamics (SORD) can be downloaded from .

References

Akaike

(1998) Information theory and an extension of the maximum likelihood principle. In: Parzen

Tanabe

Kitagawa

(eds) Selected Papers of Hirotugu Akaike. New York: Springer, pp. 199–213.

Baize

Scotti

Ritz

J-F

Ferry

Larroque

Mathot

Audin

Nurminen

Boncio

Delouis

(2020) Surface ruptures during moderate earthquakes: Is that so rare? In: Conference: GSA 2020 connects online, Montreal, QC, Canada, 26–30 October.

Bodin

Brune

(1996) On the scaling of slip with rupture length for shallow strike-slip earthquakes: Quasi-static models and dynamic rupture propagation. Bulletin of the Seismological Society of America 86: 1292–1299.

Brune

(1971) Correction to tectonic stress and the spectra of seismic shear waves from earthquakes. Journal of Geophysical Research 76: 1971.

Field

Biasi

Bird

Dawson

Felzer

Jackson

Johnson

Jordan

Madden

Michael

Milner

Page

Parsons

Powers

Shaw

Thatcher

Weldon

Zeng

(2015) Long-term time-dependent probabilities for the Third Uniform California Earthquake Rupture Forecast (UCERF3). Bulletin of the Seismological Society of America 105: 511–543.

Hanks

Bakun

(2002) A bilinear source-scaling model for M-log A observations of continental earthquakes. Bulletin of the Seismological Society of America 92: 1841–1846.

Harris

Barall

Aagaard

Roten

Olsen

Duan

Liu

Luo

Bai

Ampuero

Kaneko

Gabriel

Duru

Ulrich

Wollherr

Shi

Dunham

Bydlon

Zhang

Chen

Somala

Pelties

Tago

Cruz-Atienza

Kozdon

Daub

Aslam

Kase

Withers

Dalguer

(2018) A suite of exercises for verifying dynamic earthquake rupture codes. Seismological Research Letters 89: 1146–1162.

Huang

J-Y

Abrahamson

Chao

S-H

Sung

C-H

(2023) New empirical source scaling laws for crustal earthquakes incorporating the fault dip angle and seismogenic thickness effects. [Unpublished Manuscript].

Ida

(1972) Cohesive force across tip of a longitudinal-shear crack and Griffiths specific surface-energy. Journal of Geophysical Research 77: 3796–3805.

10.

Kanamori

(1977) The energy release in great earthquakes. Journal of Geophysical Research 82: 2981–2987.

11.

Lavrentiadis

Abrahamson

(2019) Generation of surface-slip profiles in the Wavenumber domain. Bulletin of the Seismological Society of America 109: 888–907.

12.

Lavrentiadis

Abrahamson

(this issue) A local event coordinate system for surface fault rupture. Earthquake Spectra.

13.

Leonard

(2010) Earthquake fault scaling: Self-consistent relating of rupture length, width, average displacement, and moment release. Bulletin of the Seismological Society of America 100: 1971–1988.

14.

Leonard

(2014) Self-consistent earthquake fault-scaling relations: Update and extension to stable continental strike-slip faults. Bulletin of the Seismological Society of America 104: 2953–2965.

15.

Mai

Beroza

(2000) Source scaling properties from finite-fault-rupture models. Bulletin of the Seismological Society of America 90: 604–615.

16.

Moss

RES

Ross

(2011) Probabilistic fault displacement hazard analysis for reverse faults. Bulletin of the Seismological Society of America 101: 1542–1553.

17.

Pegler

Das

(1996) Analysis of the relationship between seismic moment and fault length for large crustal strike-slip earthquakes between 1977–92. Geophysical Research Letters 23: 905–908.

18.

Petersen

Dawson

Chen

Cao

Wills

Schwartz

Frankel

(2011) Fault displacement hazard for strike-slip faults. Bulletin of the Seismological Society of America 101: 805–825.

19.

R Core Team (2022) R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

20.

Romanowicz

(1992) Strike-slip earthquakes on quasi-vertical transcurrent faults: Inferences for general scaling relations. Geophysical Research Letters 19: 481–484.

21.

Romanowicz

(1994) Comment on “A reappraisal of large earthquake scaling” by C. Scholz. Bulletin of the Seismological Society of America 84: 1675–1676.

22.

Romanowicz

Rundle

(1993) On scaling relations for large earthquakes. Bulletin of the Seismological Society of America 83: 1294–1297.

23.

Sarmiento

Madugo

Bozorgnia

Shen

Mazzoni

Lavrentiadis

Dawson

Madugo

Kottke

Thompson

Baize

Milliner

Nurminen

Boncio

Visini

(2022) Fault displacement hazard initiative database. Technical Report, Report GIRS-2021-08. Los Angeles, CA: UCLA Natural Hazards Risk and Resiliency Research Center.

24.

Scholz

(1982) Scaling laws for large earthquakes: Consequences for physical models. Bulletin of the Seismological Society of America 72: 1–14.

25.

Scholz

(1994) A reappraisal of large earthquake scaling. Bulletin of the Seismological Society of America 84: 215–218.

26.

Scholz

(1997) Size distributions for large and small earthquakes. Bulletin of the Seismological Society of America 87: 1074–1077.

27.

Scholz

(1998) A further note on earthquake size distributions. Bulletin of the Seismological Society of America 88: 1325–1326.

28.

Thingbaijam

KKS

Mai

Goda

(2017) New empirical earthquake source-scaling laws. Bulletin of the Seismological Society of America 107: 2225–2246.

29.

Wang

Day

(2017) Seismic source spectral properties of crack-like and pulse-like modes of dynamic rupture. Journal of Geophysical Research-solid Earth 122: 6657–6684.

30.

Wang

Day

(2020) Effects of off-fault inelasticity on near-fault directivity pulses. Journal of Geophysical Research: Solid Earth 125: e2019JB019074.

31.

Wang

Goulet

(2021) Validation of fault displacements from dynamic rupture simulations against the observations from the 1992 Landers earthquake. Bulletin of the Seismological Society of America 111: 2574–2594.

32.

Wang

Goulet

(2022) Constraining fault displacements for strike-slip events using physics-based simulations. In: The 12th national conference on earthquake engineering, Salt Lake City, UT, 27 June–1 July.

33.

Wells

Coppersmith

(1994) New empirical relationships among magnitude, rupture length, rupture. Bulletin of the Seismological Society of America 84: 974–1002.

34.

Wells

Youngs

(2015) Improved regression relations for earthquake source parameters. In: SSA 2015 annual meeting, Pasadena, CA, 21–23 April.

35.

Wesnousky

(2008) Displacement and geometrical characteristics of earthquake surface ruptures: Issues and implications for seismic-hazard analysis and the process of earthquake rupture. Bulletin of the Seismological Society of America 98: 1609–1632.

36.

Youngs

Arabasz

Anderson

Ramelli

Ake

Slemmons

McCalpin

Doser

Fridrich

Swan

Rogers

Yount

Anderson

Smith

Bruhn

Knuepfer

PLK

Smith

DePolo

O’Leary

Coppersmith

Pezzopane

Schwartz

Whitney

Olig

Toro

(2003) A methodology for probabilistic fault displacement hazard analysis (PFDHA). Earthquake Spectra 19: 191–219.