Finite source properties of large strike-slip earthquakes

SUMMARY

Earthquake ruptures are complex physical processes that may vary with the structure and tectonics of the region in which they occur. Characterizing the factors controlling this variability would provide fundamental constraints on the physics of earthquakes and faults. We investigate this by determining finite source properties from second moments of the stress glut for a global data set of large strike-slip earthquakes. Our approach uses a Bayesian inverse formulation with teleseismic body and surface waves, which yields a low-dimensional probabilistic description of rupture properties including the spatial deviation, directivity and temporal deviation of the source. This technique is useful for comparing events because it makes only minor geometric constraints, avoids bias due to rupture velocity parametrization and yields a full ensemble of possible solutions given the uncertainties of the data. We apply this framework to all great strike-slip earthquakes of the past three decades, and we use the resultant second moments to compare source quantities like directivity ratio, rectilinearity, average moment density and vertical deviation. We find that most strike-slip earthquakes have a large component of unilateral directivity, and many of these earthquakes show a mixture of unilateral and bilateral behaviour. We notice that oceanic intraplate earthquakes usually rupture a much larger width of the seismogenic zone than other strike-slip earthquakes, suggesting these earthquakes may often breach the expected thermal boundary for oceanic ruptures. We also use these second moments to resolve nodal plane ambiguity for the large oceanic intraplate earthquakes and find that the rupture orientation is usually unaligned with encompassing fossil fracture zones.

1 INTRODUCTION

Large earthquakes involve complex ruptures that can vary strongly between events. The characteristics of these ruptures may be controlled by the structural and tectonic characteristics of the fault zone, and understanding patterns in these ruptures may improve our understanding of the interplay between source phenomenology and the rupture zone. In particular, large strike-slip earthquakes are known to show considerable variability in rupture properties between events (e.g. Hayes 2017; Yin et al. 2021; Bao et al. 2022). Systematically characterizing this variability has the potential to yield insights into the underlying controls on the rupture process. These insights are of societal and scientific interest because these earthquakes present significant global hazard and provide unique windows into the structure and rheology of the lithosphere. Several faults known to host large strike-slip earthquakes are in close proximity to dense population centres. There is also wide speculation that the propagation behaviour and rupture dimensions are dictated by the structural (Ben-Zion & Andrews 1998; Wesnousky 2008) and rheological properties (Abercrombie & Ekström 2001; Boettcher et al. 2007) of the host fault zone. Intraplate oceanic earthquakes are particularly enigmatic, because the explanation for the weakening of the oceanic lithosphere that accommodates these events remains elusive (Lay 2019).

A general quantity for describing the space–time kinematics of earthquake ruptures is the so-called stress glut (Backus & Mulcahy 1976a), which quantifies the breakdown of linear elasticity in space and time (Dahlen & Tromp 1998). Finite-fault slip distributions, which approximate the stress glut as discretized slip on a predefined fault plane, are routinely computed for large events (Mai & Thingbaijam 2014). These solutions provide a high dimensional view of fault ruptures but in practice are challenging to compare between events due to pervasive non-uniqueness in the inverse problem, a priori fault plane parametrization, poor rupture velocity sensitivity and regularization. An alternative technique for characterizing earthquake source properties is the second-moment formulation (Backus & Mulcahy 1976a, b). Instead of approximating the stress glut as a superposition of assigned subevents, this approach involves solving for the second order polynomial moments of the stress glut, yielding a source covariance matrix that approximates the spatiotemporal extent of the source. This technique has been successfully applied in the past (Bukchin 1995; McGuire et al. 2000, 2001, 2002; Clévédé et al. 2004; Chen 2005; Meng et al. 2020), but has received far less attention than slip inversions. This low-dimensional framework makes only minor assumptions about rupture geometries and has the advantages of not requiring an explicitly parametrized rupture velocity and avoiding the discretization challenges that arise when performing slip inversions. The low dimensionality of the solution also facilitates computation of these moments with a Bayesian approach (Atterholt & Ross 2022), which can provide uncertainty estimates crucial for comparing the source processes of different earthquakes.

Our contributions to this study are as follows. We compute second moments for all of the M_w ≥ 7.5 strike-slip earthquakes of the past three decades using a Bayesian inference approach. We use this catalogue to establish baselines for the range of values observed globally and compare values between events, subgroups and other tectonic features. From these analyses we conclude that (i) large strike-slip earthquakes almost always show unilateral or a comparable amount of unilateral and bilateral directivity behaviour, (ii) that large intraplate oceanic earthquakes usually rupture over a much wider depth range and (iii) that oceanic intraplate strike-slip earthquakes are not systematically aligned with fossil fracture zones.

2 PRELIMINARIES

The stress glut is a tensor field representing the expected stress due to the application of Hooke’s law to inelastic strain in a body (Backus & Mulcahy 1976a, b). The stress glut is a useful source characterization quantity, because the stress glut is identically zero outside the source region and can be used to compute displacements anywhere on Earth resulting from an arbitrary source (Dahlen & Tromp 1998). The dimensionality of the stress glut can be significantly reduced by assuming the source mechanism does not change throughout the rupture:

Here, |$\pmb {\Gamma }$| is the stress glut, |$\mathbf {\hat{M}}$| is the normalized mean seismic moment tensor and f is a scalar function of position |$\pmb {\xi }$| and time τ. The second moment formulation is defined by taking the second central moment of the scalar stress glut rate function (⁠|$\dot{f}$|⁠) with terms for the spatial and temporal components. The equation for these moments is given by:

where |$\pmb {\xi }$| and τ are position and time, |$\pmb {\xi }^{c}$| and τ^c are the centroid position and centroid time, which are the first moments of the distribution, and m and n are the spatial order and temporal order of the moment. The central moments of order m + n = 2 correspond to the covariance of the source. Specifically, |$f(\pmb {\xi }^{c},\tau ^{c})^{(2,0)}$| is the spatial covariance, |$f(\pmb {\xi }^{c},\tau ^{c})^{(0,2)}$| is the temporal variance and |$f(\pmb {\xi }^{c},\tau ^{c})^{(1,1)}$| is the spatiotemporal covariance. The distribution |$\dot{f}$| is defined by the superposition of its polynomial moments, and at low-enough frequencies, the contribution of moments of order three or greater may be approximated as zero. Under this approximation, the second moments can be linearly related to displacement:

where d is a vector of the difference between the measured displacements and the theoretical Green’s functions, F is a forward propagation matrix of spatial and temporal integrals and derivatives of the Green’s tensor weighted by the components of the moment tensor M, and p is a vector that contains the independent parameters of the second order stress-glut moments.

Since the standard deviation of a distribution gives a low-dimensional estimate of the width of a distribution, these second moments, which are the covariance of the stress-glut, can be used to compute low-dimensional measures of the volume, duration and directivity of a source distribution. In particular, we define dimensional quantities of the source that describe the shape of the stress-glut distribution about the centroid. These are:

Here, |$r^{c}(\mathbf {\hat{n}})$| is the distance from the centroid in the direction of a unit vector |$\hat{n}$| that defines a spatial deviation ellipsoid in which most of the moment was released. The characteristic spatial deviation of the source is given by |$L^{c}=2 r^{c}(\pmb {\eta })$|⁠, where |$\pmb {\eta }$| is the principal eigenvector of |$\pmb {\dot{f}}^{(2,0)}$|⁠. t^c is the characteristic temporal deviation that captures a time interval in which most of the moment was released. v₀ is the average instantaneous velocity of the rupture centroid. These quantities together provide a physically interpretable, low-dimensional estimate of a rupture’s spatiotemporal behaviour (Backus 1977; Silver & Jordan 1983).

From the aforementioned quantities, we can compute ensembles of parameters that may further illuminate potential differences between ruptures. In particular, we inspect four derived parameters in this study: rectilinearity (R), directivity ratio (α), average moment density (⁠|$\bar{m}$|⁠) and vertical deviation (Z). Rectilinearity is a measure of the degree of elongation along an axis and has been used in seismology to evaluate particle motion (Vidale 1986; Jurkevics 1988). We instead use this measure to evaluate the elongation of ruptures along the principal axis, and define it so the values are bounded between 0 (spherical source) and 1 (linear source):

The variables λ₁, λ₂ and λ₃ are the eigenvalues of the spatial second moment, in order of largest to smallest and yield estimates of the dimensions of the source along its principal axes. The directivity ratio provides an estimate of the degree of directivity of a rupture by comparing the average centroid velocity to the maximum possible average centroid velocity, and has been used in the second moment literature previously (McGuire et al. 2002):

The average moment density comes from the moment tensor density formulation (Aki & Richards 2002) and compares the scalar moment, M₀, to the volume of the ellipsoid representation of the spatial second moments:

As is apparent, the equation for |$\bar{m}$| is very similar to analytic solutions for the stress drop on planar sources (e.g. Eshelby 1957) and has units of stress. Finally, we define a measure of the vertical extent of the source:

Here, |$\pmb {\gamma }$| is the unit vector pointing in the direction with the maximum vertical component. Each of these quantities yield estimates of the source dimensions that have historically been difficult to obtain. In this framework, we can systematically solve for full posterior distributions on these parameters and consequently perform error-informed comparisons.

Because we will subsequently apply this inversion framework to large strike-slip earthquakes, an important point about the assumption of a constant source mechanism should be made here. Several large strike-slip earthquakes, such as the 2012 Wharton Basin and the 2018 Gulf of Alaska Earthquakes, have been shown to rupture multiple conjugate, near-orthogonal faults. This presents no problem in our framework because the faults are nearly orthogonal and have opposite senses of slip, making the source mechanism approximately constant between faults. There is no assumption that the rupture be constrained to a single fault. The orthogonality of the faults in these ruptures is likely not a coincidence, and may be systematic for large fault zones (Thatcher & Hill 1991; Scholz & Choi 2022); thus we may expect that for other ruptures of this type, the assumption of a constant source mechanism is also reasonable. This assumption is less valid for faults with significant curvature (Shimizu et al. 2020), ruptures in which the slip direction changes significantly, or ruptures of multiple faults that are highly non-parallel and non-orthogonal. Though this assumption is a limitation of this methodology, this methodology is strengthened by avoiding constraints of slip to a user-defined fault plane that is typically necessary for fault slip distributions.

3 METHODS

3.1 Data and preprocessing

We use the Global Centroid Moment Tensor (gCMT) catalogue to select large strike-slip earthquakes of the past three decades (Ekström et al. 2012). To find these events, we search the gCMT catalogue for events with M_w ≥ 7.5 and with nodal axis plunges greater than 45°. We then manually evaluate these events to ensure that each event shows predominantly strike-slip behaviour, resulting in the set of events shown in Fig. 1 and Table 1. We choose events of this magnitude because the frequency to which the waveforms for these events need to be fit to resolve stress-glut second moments is inversely related to the duration of the event. Consequently, larger events with correspondingly larger durations may be fit with theoretical Green’s functions at lower frequencies. We directly fit the displacement waveforms from Global Seismographic Network (GSN) stations for data consistency between events. The waveforms used in these inversions are SH and R1 waveforms selected using travel times from the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981) specific to the selected frequency bands in this study. The inclusion of the P waveforms would further constrain the posterior distributions, but we decide not to include these because this would require the computation of theoretical Green’s functions at much higher frequencies. We select 200-s windows around the SH waveforms and 700-s windows around the R1 waveforms. Waveforms at distances less than 30° are excluded to avoid the convolution of the SH waves with the R1 waves, because at short distances these phases have not diverged enough to select SH windows with no R1 waveforms present. An example set of waveforms for the 1999 Izmit earthquake are shown in Fig. 2. We compute the Green’s tensors using the gCMT moment tensor solutions (Ekström et al. 2012). We use the spectral element method software Salvus (Afanasiev et al. 2019) in combination with the 3-D earth model S362ANI+M (Moulik & Ekström 2014) to compute the Green’s tensors. The derivatives of the Green’s functions needed for the forward propagation matrices are computed using a centred finite difference approximation. For the spatial finite difference, we compute the wavefield for a grid of source locations centred on the centroid solution.

As in Atterholt & Ross (2022), we select the frequency bands used in these inversions using duration estimates for each event. In particular, we consider the observation that the contribution of moments of order m + n is approximately proportional to (D/P)^{m + n} (where (D/P) < 1), where D and P are the event duration and period of the signal, respectively (Backus 1977). We select a frequency band such that (D/P)² > 0.05 and (D/P)³ < 0.05 in order to balance maximizing the second moments relative to the zeroth and first moments while minimizing the contribution of moments of order three and greater. For an estimate of duration, we use the empirical relationship used for the gCMT catalogue (Ekström et al. 2012). These frequency bands are given in Table S1. We subsequently show that this relationship provides a fairly reasonable estimate for the source duration for these events. With the aforementioned inequality and duration estimates, we can compute a frequency band for filtering the data, and we filter the waveforms using a Butterworth filter. Once filtered, we visually inspect the fit of the theoretical Green’s functions to the data for quality control.

3.2 Inversion

For the second moment inversion, we follow the procedure outlined in Atterholt & Ross (2022), and we summarize the method here. In essence, our objective is to use the relationship described in eq. (3) to invert for a model vector of stress glut moments that better fits displacement waveforms by accounting for the finiteness of the source. An example of the waveform fits for the 1999 Izmit, 2002 Denali and 2018 Palu earthquakes using the subsequently described inversion are shown in Figs 2, S1 and S2, respectively. The azimuthal distributions for these earthquakes, which are representative of the azimuthal distributions for the earthquakes used in this study, are shown in Fig. S3. We solve this inverse problem using a Bayesian formulation that produces an ensemble of potential solutions given the uncertainty of the inverse problem, outlined in detail by Atterholt & Ross (2022). The posterior distribution for this problem can be written using the relationship:

where σ is a parameter that is updated every sample and quantifies the uncertainty by estimating the inability of the model to fit the data (Gelman et al. 2010). p(p) is the prior distribution on the model parameters; in this study we use uninformed priors. We compute the likelihood using a multivariate normal distribution:

where Σ is a block diagonal covariance matrix that incorporates the temporal correlation structure, which is a function of the parameter σ and the minimum period of the bandpass filter. A limitation of this methodology is that the parameter σ accounts for random, but not systematic, errors in the Green’s functions.

Prior studies on stress-glut second moment inversions (Bukchin 1995; McGuire et al. 2001) have applied the constraint that since second moments constitute a covariance matrix, only the solutions that ensure the second moments are positive definite are valid, or:

To ensure our samples of the posterior are positive definite, we take advantage of the Cholesky Factorization Theorem, which states that every symmetric, positive definite matrix can be represented as the product of a lower triangular matrix and its transpose. For X, there thus exists a matrix L such that X = LL^T. We thus sample over L and construct X from L when computing the likelihood of the sample. This framework allows for the inclusion of prior distributions on the source model; for example, it would be possible to impose a prior on the azimuth of the principal eigenvector or the size of the characteristic values using information from previous source characterizations or empirical scaling relationships. In this study, we use uninformed priors to maximize the independence of this catalogue of solutions from other source studies. We probe the solution manifold using Hamiltonian Monte Carlo (HMC) sampling (Neal 2010), an instance of Metropolis Hastings sampling that uses an analogue from Hamiltonian Dynamics. We choose HMC sampling because our parameter space is large, and HMC sampling applies gradient information to sample large parameter spaces efficiently. To take these gradients, we use reverse-mode automatic differentiation (Innes 2019). In total, we take 10 000 samples, and in each subsequent representation of solution ensembles, we show the final 500 unique samples. In Atterholt & Ross (2022), we demonstrated that this parametrization was appropriate to conservatively achieve convergence for this type of inverse problem using the Gelman–Rubin diagnostic (Gelman & Rubin 1992). The computational time for this inversion varied according to the number of waveforms used, but in general was executed for each event within between 2 and 10 hr on a single CPU core.

4 RESULTS AND DISCUSSION

4.1 Patterns in rupture behaviour

Using the definitions in eqs (4)–(8), we compute ensembles of quantities that capture low dimensional rupture characteristics. We summarize these ensembles in Table 2 and show the full posterior of these quantities in Figs S4, S5 and S6. In Fig. 3, we plot the ensembles of temporal deviations, which show good agreement with the empirical magnitude–duration relationship used to determine the frequency bands used in this study. The initial duration estimates control the quality of the solution by reducing biases from the error of the moments of order m + n ≤ 1 and from the contributions of the moments of order m + n ≥ 3, but this does not suggest a causal relationship between the initial duration estimates and the second temporal moments. These moments thus show that this empirical relationship is a reasonable model for these data.

Also shown in Fig. 3 is a comparison between the characteristic spatial deviation of the rupture, L_c, and the propagation length of the rupture centroid |$||\pmb {v}_{0}||\cdot t^{c}$|⁠. The ratio between these two quantities reflects the degree of directivity of the rupture (McGuire et al. 2002). We separate these quantities into three categories of our design: bilateral (⁠|$||\pmb {v}_{0}||\cdot t^{c}/L_{c}\le \frac{1}{3}$|⁠), mixed (⁠|$\frac{1}{3}\le ||\pmb {v}_{0}||\cdot t^{c}/L_{c}\le \frac{2}{3}$|⁠) and unilateral (⁠|$||\pmb {v}_{0}||\cdot t^{c}/L_{c}\ge \frac{2}{3}$|⁠). As is made clear in Fig. 3, almost no ruptures fall in the bilateral category, suggesting that the vast majority of these ruptures have a large component of unilateral directivity, which is consistent with previous results (McGuire et al. 2002; Ross et al. 2020). But, while many of these ruptures fall in the unilateral category, suggesting these ruptures favour unilateral directivity over bilateral directivity, many of the ruptures also fall in the mixed category. Mixed directivity in ruptures may point to a more complicated rupture process and are difficult to achieve by chance. McGuire et al. (2002) computed the expected distribution of directivity ratios given a 1-D predefined fault with uniform slip and hypocentre location and found that the distribution is heavily skewed towards unilateral ruptures. We reproduce this distribution and compare it to our observed distribution of directivity ratios. For this comparison, we choose the median as the test statistic, because the median is a good estimate of the central tendency of skewed distributions. We create a distribution of possible median values for our observed distribution by computing them for random ensemble member draws from the source directivity distributions and comparing these values to the distribution of medians from random draws of the same size from the simulated distribution. This workflow is summarized in Fig. S7. We find that the vast majority of the ensemble draws from our observed distribution are unlikely to be drawn from the simulated distribution, suggesting the difference between these distributions is significant. This would suggest that the simple hypothesis that ruptures start at random positions on a predefined fault with uniform slip cannot explain our data. In reality, fault ruptures are more complex than this simple model, and a number of factors contributing to this complexity may explain the difference between these distributions. For example, ruptures may preferentially start near high strength regions (Mai et al. 2005; Melgar & Hayes 2019) which may be biased towards the centre of ruptures (e.g. Manighetti et al. 2005) and thus bias ruptures towards bilateral behaviour. With additional observations and tests against simulations of more complicated ruptures, this type of data may supplement arguments for specific rupture behaviours that prefer bilateral propagation.

Finally, in Fig. 3, we plot twice the characteristic spatial deviation of the rupture against magnitude to observe the scaling relationship present in this catalogue. Additionally, we plot best-fitting solutions for L model (Scholz 1982) and W model (Romanowicz 1992) scaling relationships and empirical scaling relationships from previous studies that use independent slip distribution catalogues (Wells & Coppersmith 1994; Mai 2000; Blaser et al. 2010; Thingbaijam et al. 2017). We find the trend of our L_c values, in general, matches the trends of the scaling relationships, but there is not sufficient coverage along the magnitude axis to determine with which scaling relationship this catalogue is best aligned. However, among the intraplate oceanic events, which have more diverse observations with magnitude, the trend between M_w and L_c appears to deviate from these of the scaling relationships; the significance of this is described subsequently.

We compute ensembles for the quantities described in eqs (5)–(8) for all the events considered in this study and plot the kernel density estimates of the underlying probability distributions of these ensembles in Fig. 4. Almost all rectilinearity values for these events are well-above 0.5, suggesting that most of these ruptures are elongated along a single dimension. Large strike-slip earthquakes are large-enough to be constrained by the width of the seismogenic zone, and thus elongation of the spatial deviation ellipsoid in the strike-parallel direction is expected. The directivity ratio, also represented in Fig. 3 and discussed in the previous section, yields a wide distribution that is biased towards values closer to one. Average moment densities for these events vary by over two orders of magnitude (1–100 MPa) with a median value near 10 MPa. The distributions of these values should be interpreted carefully, as the average moment density estimates are computed using the smallest eigenvalue, λ₃, of the spatial covariance matrix; this value is poorly constrained by our data. We expect λ₃ to be small for most ruptures, but since the data do not constrain λ₃ well, the posterior distributions of λ₃ can be very broad above zero, and thus skew the distributions of source volume upward. This issue may be addressed by assigning a well-justified prior distribution on the size of λ₃, but we leave this for future work. The distribution of vertical deviations is mostly below 20 km with a noteworthy tail extending up to 40 km. The ensembles of these parameters for each event are shown in Figs S4, S5 and S6 and are summarized in Table 2.

We separate these events into three categories: continental strike-slip, oceanic interplate and oceanic intraplate. For oceanic events, we make classifications by cross-referencing the event locations with the plate boundary types outlined in Bird (2003). This classification is not always simple. For example, the 2003 Carlsberg Ridge event ruptured close to the active spreading centre and is expected to be aligned with the spreading-associated transform faults in the vicinity (Antolik et al. 2006), but it has an opposite sense of slip as expected for these faults. We nonetheless classified this event as interplate because of it’s proximity to the plate boundary assuming the true nodal plane is that aligned with the nearby active transform faults. Classifying this event as intraplate does not change the conclusions of this study. We classify all continental strike-slip events together to ensure that there are enough events in the group to evaluate potential systematic behaviour. These classifications are given in Table 1. The rectilinearity and directivity of these events are too variable within groups to make statements regarding differences between groups, but these ensembles demonstrate the potential of this second-moment framework to objectively obtain quantities that potentially illuminate rupture features that have historically been difficult to obtain independently of the necessary assumptions associated with fault slip distributions. Some differences in average moment densities are evident in these distributions; namely, the ensembles of average moment densities for oceanic interplate earthquakes create a distribution that is much narrower than that of other strike-slip earthquakes, because the other types of events include the contributions of several earthquakes with ensembles of exceptionally high average moment density values.

The depth-dependence of slip distributions has historically been difficult to constrain because surface observations typically have poor sensitivity to variability with depth (e.g. Duputel et al. 2014). As is shown in Fig. 4, the uncertainty of the vertical deviation parameter in this study is quite high. However, because we estimate the full posterior of our solutions, we have a good estimate of the uncertainty present in these solutions. We can thus make error-informed conclusions while accounting for the high uncertainty inherent to resolving parameters related to the source distribution with depth. The distributions of vertical deviations of these events suggest a clear difference between oceanic intraplate earthquakes and other strike-slip earthquakes. The oceanic intraplate earthquakes exhibit much wider ruptures than other strike-slip earthquakes; this has been suspected before, particularly for Wharton Basin events (Aderhold & Abercrombie 2016). The observations in this study provide some evidence for this difference globally. Mindful that the vertical deviation measure only accounts for ±1σ of the source distribution’s vertical width, it is reasonable to expect that the source distribution extends over an even larger width than the estimated vertical deviation of the rupture. Given that these intraplate oceanic earthquakes generally rupture colder lithosphere, some difference in vertical deviation between rupture types is to be expected. For earthquakes that ruptured oceanic crust with measured age, we can account for the expected thermal differences by comparing the age of the oceanic crust encompassing the centroid, taken from Seton et al. (2020), to the expected thermal profile of a half space cooling model with an ambient mantle temperature of 1350° and thermal diffusivity of 10⁶ m² s⁻¹ as in Aderhold & Abercrombie (2016). Noting that the vertical extents are estimates of the variance of the stress glut from the centroid, we can estimate the depth extent of these earthquakes by summing the centroid depth with quantities proportional to the vertical deviation of the source.

We plot these distributions of depth extent assuming the depth extents given by the conservative estimates of |$\pmb {\xi }^{c}(z)+\frac{1}{2}Z$| and |$\pmb {\xi }^{c}(z)+\frac{\sqrt{2}}{2}Z$| in Fig. 5. The |$\pmb {\xi }^{c}(z)+\frac{1}{2}Z$| values assume slip only occurred within a single standard deviation of the source distribution below the source centroid, and the |$\pmb {\xi }^{c}(z)+\frac{\sqrt{2}}{2}Z$| values would yield the exact depth extents for source distributions with uniform slip. We note here that there may be some inaccuracies in the Green’s functions due to the crustal model when generating sources in oceanic crust. These inaccuracies would increase the modelled excitation of the surface waves (Abercrombie & Ekström 2003) which we would expect to increase the radial spatial gradients incorporated in eq. (3). This bias would serve to decrease the estimated depth extent of the earthquakes, and thus would not change subsequent conclusions of this study. Additionally, because the Wharton Basin (3) event was an aftershock that occurred only 2 hr after the Wharton Basin (2) event, the Wharton Basin (3) event waveforms show long-period noise that is likely causing the high uncertainty for that event in Fig. 8, this solution should thus be considered in that context. An estimate of how well resolved the vertical extent values are can be made by comparing the ratio of λ₂/λ₃. Since λ₃ is expectedly small, larger values of λ₃ may be interpreted as an approximate noise level. These ratios are plotted in Fig. S8. In general, while some ensemble members for some sources are lower, these ratios greatly exceed 2; however, these values are not purely estimates of the signal-to-noise ratio, as physical characteristics of the source can contribute to elevated λ₃ values.

As is apparent in Fig. 5, for interplate events these values generally fall above the 800° isotherm and for intraplate events these values generally fall below the 800° isotherm. While there is some variability in estimates of the limiting temperature for slip on oceanic transform faults (e.g Wiens & Stein 1983; Abercrombie & Ekström 2001), the estimates of depth extent of the slip for oceanic intraplate earthquakes frequently exceed even the highest expected slip-limiting temperatures. These observations are supported by the L_c values, which are better constrained, for the intraplate oceanic ruptures shown in Fig. 3. These values shows that the intraplate events, L_c, when plotted against magnitude, follow a trend that is distinct from expected scaling relationships for strike-slip earthquakes which are expectedly restricted by their widths. These results motivate the consideration of a mechanism for coseismic slip below the brittle–ductile transition zone for large intraplate earthquakes. One potential mechanism for such a large amount of moment released below the brittle-ductile transition zone, proposed as an explanation for the slip in 2012 Wharton Basin Earthquake, is a large component of deep viscous failure, rather than frictional failure, resulting from a runaway feedback system due to heating of the shear zone (Kelemen & Hirth 2007; McGuire & Beroza 2012). This mechanism, however, would only account for moment released between 600° and 800°. Another potential explanation is that diffuse deformation zones in the interior of oceanic lithosphere may be hydrated at depth, altering the rheological profile of the fault (Bishop et al. 2023) and thus allowing for deeper rupture.

4.2 Orientations of ruptures and crustal fabric

For large strike-slip earthquakes that occur on plate boundaries or in continental crust, the fault associated with the majority of slip is typically well resolved or easy to infer from the focal mechanism, and in most cases, the processing involved in computing a fault slip distribution helps discern the true nodal plane (Hayes 2017). However, for intraplate oceanic earthquakes, the true fault plane is often not easily resolved (e.g. Nettles et al. 1999; Abercrombie et al. 2003; Robinson 2011; Meng et al. 2012; Lay et al. 2018). The determination of the primary fault plane for these earthquakes (if one exists) is an important problem, because the enigmatic rupture of these very large events may broaden our understanding of the strength of the oceanic lithosphere and the nature of these earthquakes. A frequent assumption is that these events rupture fossil fracture zones, crustal fabric imparted by ocean spreading. However, some detailed studies of several events suggest rupture plane strikes that disagree with those of nearby fracture zones or rupture several nearly orthogonal fault planes (Nettles et al. 1999; Meng et al. 2012; Lay et al. 2018), making the fracture zone hypothesis an unlikely universal explanation for weakening of the oceanic lithosphere. Fault slip distributions are often used to resolve nodal plane ambiguity by comparing the data fitting between different fault plane parametrizations (Hayes 2017), and this too has been used to suggest that oceanic intraplate ruptures are often unaligned with fossil fracture zones in Wharton Basin, for example (Aderhold & Abercrombie 2016). However, resolving fault plane ambiguity this way requires initial assumptions about the fault geometry and rupture process and ultimately a subjective, albeit data-informed, decision about the true fault plane. Backprojection of the energy radiated from the source has been used to successfully image the rupture surfaces of earthquakes (e.g. Ishii et al. 2005; Meng et al. 2012; Ruppert et al. 2018), but estimation of uncertainty with this method is difficult, and the quality of backprojection recovery is heavily dependent on the array used for backprojection’s density of coverage and distribution of azimuths and distances (Kiser & Ishii 2017). The second moment formulation can be used to supplement these prior studies by providing an objective and probabilistic assessment of the orientation of the rupture.

The second moment formulation can be used to resolve the nodal plane ambiguity of a point source moment tensor solution. Since we expect the slip variance to be maximized along the rupture plane, the principal eigenvector of the spatial second moment should be approximately aligned with the strike of the true nodal plane. We can verify this by comparing the strike of the principal eigenvectors of the ensembles of second spatial moments with true nodal plane of earthquakes for which this can be easily inferred. For all continental strike-slip earthquakes and interplate oceanic events, we can be reasonably confident of the true nodal plane in the presence of an observed surface rupture or nearby plate boundary information. For each of these events, we plot (Fig. 6) the difference in angle between the strike of the inferred nodal plane and the strike of the principal eigenvectors for the corresponding ensemble of spatial second moments. The median and vast majority of the angular differences of each ensemble are less than 45°, suggesting that this methodology is an effective tool for approximating the strike of the true rupture plane. We avoid interpreting the strike of the principal eigenvectors at a more granular level, because, as is apparent in Fig. 6, some misalignment of the principal eigenvector and the true nodal plane is common. The threshold of 45° used to determine alignment is large enough to exceed potential errors in the principal eigenvector strike.

Intraplate oceanic ruptures have been observed to rupture multiple, near-orthogonal faults throughout the course of the rupture (Meng et al. 2012; Lay et al. 2018; Yamashita et al. 2021). As stated earlier, for these ruptures, the source mechanism remains approximately constant because the rupturing faults are nearly orthogonal, and thus the geometric assumption of this methodology remains valid. The stress glut covariance will track the distribution of moment release throughout the course of the rupture, and so for multi-fault ruptures the principal eigenvector may be unaligned with the rupture planes. For an intraplate oceanic earthquake, if this results in the principal eigenvector being unaligned with fossil fracture zones, we consider this a valid example of weakening resulting in a rupture that cannot be explained solely by fossil fracture zones, because such a rupture necessarily propagates on multiple faults, some of which will be oriented at high angle to local fossil fractures.

We apply this technique to the intraplate oceanic earthquakes that are situated near mapped fracture zones. We estimate the fracture zone strike in the area by considering the strikes of the fracture zones within 2° the event centroid. We compute fracture strikes using the shapefiles provided by Seton et al. (2014) and Wessel et al. (2015). For each event, we take the angular difference between the entire ensemble of principal eigenvectors and every potential fracture strike to create a distribution of potential angular differences. We plot these distributions in Fig. 6. We find that only 3 of the 7 distributions have any members within 45° of the fracture zone strike, suggesting that, more often than not, these ruptures are unaligned with the surrounding ocean fabric. Indeed, this methodology cannot decisively say whether or not the ruptures with principal eigenvectors that fall within 45° of the fracture zone strike actually ruptured the fracture zone, but the large disagreement between principal eigenvector strike and fracture zone strike for the other ruptures suggest that it is very unlikely that fossil fracture zones dictated the primary rupture plane. This motivates new explanations for weakening in the interior of the oceanic lithosphere that could host these large events, such as the explanation that diffuse deformation away from the plate boundary may induce faulting independent of pre-existing ocean fabric.

4.3 Comparisons with previous studies

Many of the events included here are well-studied, and we thus devote considerable space to evaluating our solutions in the context of prior literature. In Fig. 7, we show projections of the ensembles of spatial deviation ellipsoids for a subset of the events used in this study. These deviation ellipsoids are defined such that the axes are the principal components of the second spatial moments scaled by the standard deviation along that component. This is equivalent to finding an ellipsoid that captures a single standard deviation of the source distribution in all directions. We selected these events for display here because there exist readily available fault slip distributions that were constrained using high-quality geodetic observations and thus have a data constraint that is independent from seismic observations (Mai & Thingbaijam 2014; Socquet et al. 2019). We plot the ellipsoid representation and expected vertical extent for a larger subset of events against all the events with available peer-reviewed fault slip distributions in the SRCMOD database in Figs S9 and S10. We compute the expected spatial deviation ellipsoids for these fault slip distributions and compare them to our solutions to assess the strengths and shortcomings of our inversion. In general, there is good agreement between our solutions and that of the fault slip distributions in terms of the rupture azimuths, the ruptures’ vertical deviations, and, factoring in the subsequently discussed uncertainties, the overall spatial deviations of these ruptures. Notably, this agreement holds for the 2012 Wharton Basin earthquake for which, despite high variability between slip distributions, a large component of approximately east–west moment release is resolved in each study. All of these comparisons are made more difficult by the intrinsic variability between slip models due to the use of different data sets, inversion techniques, fault parametrization and assumptions about the rupture process. For example, in the case of the 1999 Izmit Earthquake, the difference between the largest and smallest characteristic spatial deviations of the models included in this study is over 70 per cent. Given these factors, the solutions presented in this study provide a valuable contribution, as this methodology avoids the a priori parametrizations and assumptions present in slip distributions, and is applied uniformly to all of these earthquakes.

One key difference between our results and those of the fault slip distributions is that our ellipsoids are much wider in the fault normal direction than those computed from the slip distributions; this is in fact expected behaviour. This difference is a consequence of the fact that the fault slip distributions parametrize slip on planar/curved surfaces rather than 3-D volumes, resulting in spatial deviation ellipsoids that are artificially much narrower than they would be if the inversions were solved without this constraint. We also note that for the 1998 Ceram Sea Earthquake, the strike of our characteristic ellipsoid is at high angle to the solution from Hayes (2017). Selection of the true nodal plane for this event which, as reported in Hayes (2017), was an intraplate earthquake in a poorly instrumented area with localized slip near a single asperity, is difficult. Our methodology does not require us to define a rupture plane a priori, and instead characterizes a moment distribution directly from the waveforms. We expect this moment distribution to be aligned with the true nodal plane. As we discussed previously, this methodology is particularly useful for resolving nodal plane ambiguity for such challenging cases.

Another important difference is that for the events that are expected to have the greatest spatial length, our inversion resolves characteristic spatial deviations that are sometimes lower than those of the fault slip distributions. There are several potential reasons for these discrepancies, including spatial smoothing due to regularization of the fault slip distributions and imaged slip artefacts due to inaccuracies in the a priori assumptions of the inverse model (Graves & Wald 2001; Gombert et al. 2018; Ragon et al. 2019b; Ortega-Culaciati et al. 2021). Both of these factors would contribute to increases in the dimensions of spatial second moments computed from slip distributions because they can artificially increase the spread of inferred slipping regions. These factors are likely compounded for large ruptures because the a priori fault surface dimensions of the slip model are set to be conservatively large to ensure the incorporation of the entire rupture process, and thus extend inversion artefacts farther from the centroid. This is acceptable for a qualitative description of the rupture process; however, this slip can significantly bias the spatial second moments of these slip distributions. For example, in the case of the 1998 Balleny Islands Earthquake, the authors note that deep slip in this model is poorly resolved and potentially overestimated (Antolik et al. 2000). We thus excise a segment of deep slip to the west of the centroid of their slip model, which greatly reduces the size of the spatial ellipsoid to be much closer to our second moment solution. We can systematically alleviate this issue by removing slip patches with less than 10 per cent of the peak slip, a very conservative estimate of the amplitude of poorly constrained slip in these models, from each of these slip distributions prior to computing their second moments. We plot these corrected second moment ellipsoids and vertical deviations in Figs S11 and S12. This accounts for a large amount of the characteristic spatial deviation discrepancy for some of these events, particularly the Denali and Balleny Islands Earthquakes.

Another potential issue with some of these slip distributions is that, if they are performed using only geodetic data with sparse temporal resolution, such as Interferometric Synthetic Aperture Radar (InSAR), they will include post-seismic slip into their solutions and resolve coseismic slip with high uncertainty (Weston et al. 2012; Twardzik et al. 2019; Ragon et al. 2019a). This is potentially the explanation for a large discrepancy between our solution for the 2001 Kunlun Earthquake and that of Lasserre et al. (2005), who perform their inversion using coseismic InSAR data that incorporate between 1 and 6 months of post-seismic slip. For example, there is a large amount of slip along the Taiyang Lake Fault that is resolved using data far from the fault with a large temporal window between frames (3.7 yr before the earthquake and 6 months after the earthquake). Excising the slip from this fault segment almost completely removes the difference in second moments between their solution and our solution. We are not suggesting that there was no slip on the Taiyang Lake Fault during this earthquake, as other slip distributions suggest a low-amplitude subevent to the west (e.g. Antolik et al. 2004), but this comparison illustrates the high degree of sensitivity these second moments computed from slip distributions have to poorly constrained and potentially overestimated slip far from the centroid.

There are also many slip models not included in this direct comparison. While we cannot discuss all of them here, we note several key differences between these stress glut variances and finite fault slip distributions that are important in contextualizing these results. The stress glut covariance only yields an estimate of a parameter that summarizes the distribution of slip, and does not estimate the full slip on the fault. For example, although several hundred kilometres of the Bering Fracture Zone are expected to have ruptured during the Komandorsky Islands earthquake (Lay et al. 2017), the slip distribution suggests highly concentrated moment release along a much shorter segment of the fault. This type of distribution would yield a source spatial variance that is much smaller than the full length of slip on the rupturing fault. The temporal deviations should be interpreted similarly. Additionally, the average centroid velocity describes the integrated velocity behaviour of the entire rupture, and thus should not be confused with instantaneous rupture speed. For example, a rupture that experienced supershear propagation could have a much lower average centroid velocity if the rupture had a bilateral component or if the rupture was slower before or after the supershear period.

Finally, solutions for the Tibet, Balleny Islands, Ceram Sea and Izmit earthquakes were computed with stress glut second moment inversion frameworks in prior studies (McGuire et al. 2000, 2002; Clévédé et al. 2004). The characteristic quantities from these studies are shown in Table S2. In general, while these solutions agree for many parameters, there are several differences between these sets of solutions. The Izmit Earthquake is the most useful earthquake for comparison, because second moments were computed for this event in three different frameworks, and this event is the best documented of the overlapping events. In short, the spatial deviations are largely consistent between studies, McGuire et al. (2002) reports a much narrower temporal deviation than Clévédé et al. (2004) and this study, and this study reports a much lower average centroid velocity than McGuire et al. (2002) and Clévédé et al. (2004). For this event, the longer temporal deviation estimates and the expectation of near-bilateral directivity reported in this study are consistent with studies that investigated the time-dependent slip of this earthquake (Yagi & Kikuchi 2000; Delouis 2002; Bouchon 2002; Li 2002), so we believe our framework yields the preferred solution for this event. Again, for the other events, there are many commonalities but several differences between solutions. These differences may arise from the differences in priors and parametrizations between studies. For instance, McGuire et al. (2000) and McGuire et al. (2002) allow the lower order moments to vary, constrain the radial component of the spatial covariance, use a fixed frequency band for all events, and use a different earth model and distinct phases. We choose to fix the lower moments and allow the second moments to vary freely, because we consider the centroid moment tensor solutions for these large events to be very well-constrained relative to our knowledge of the finite properties of these events. Our selection of the frequency band and earth model are meant to improve the isolation of the contribution of the second moments of these ruptures.

5 CONCLUSIONS

By computing ensembles of stress glut second moments for all large strike-slip earthquakes of the past few decades, we illuminated several patterns in these ruptures that suggest predominant behaviours and variability with tectonic environment. Our results show that large strike-slip ruptures usually have a large component of unilateral behaviour, with many ruptures exhibiting evidence for complicated rupture propagation sequences. We also observed that intraplate oceanic events have wider ruptures than other large strike-slip earthquakes, potentially illustrating a systematic behaviour of rupturing below the expected brittle-ductile transition zone depth. Finally, we show that by resolving the approximate rupture plane strike of major oceanic intraplate earthquakes, most of the earthquakes show no alignment with the fossil fracture zones in the area. This suggests that the assumption that large oceanic intraplate earthquakes reactivate fossil fracture zones is usually incorrect.

ACKNOWLEDGMENTS

This work was partially funded by the National Science Foundation’s Graduate Research Fellowship Program (GRFP) with grant number DGE-1745301. The authors would like to thank Dr Jeffrey McGuire and an anonymous reviewer for their insightful comments which greatly improved this paper and the editor Dr Sidao Ni for facilitating the review process. The authors would like to thank Dr Théa Ragon for providing very useful suggestions regarding the uncertainties of fault slip distributions. The authors would also like to thank Dr Joann Stock for her insightful comments regarding intraplate oceanic ruptures.

DATA AVAILABILITY

The map shown in Fig. 1 was created using The Generic Mapping Tools (GMT), version 6 (Wessel et al. 2019), which is available at https://www.generic-mapping-tools.org/. The centroid moment tensor solutions used in this study are from the Global Centroid Moment Tensor (gCMT) catalogue (Ekström et al. 2012) which is accessible online at https://www.globalcmt.org/. The theoretical Green’s functions were computed using Salvus (Afanasiev et al. 2019), which is available at https://mondaic.com/. The waveform data in this study are from the Global Seismographic Network operated by both the Albuquerque Seismological Laboratory (IU: IRIS/USGS; https://doi.org/10.7914/SN/IU) (Albuquerque Seismological Laboratory (ASL)/USGS 1988) and the Scripps Institution of Oceanography (II: IRIS/IDA; https://doi.org/10.7914/SN/II) (Scripps Institution Of Oceanography 1986). These waveforms may be accessed through the IRIS Data Management Center (DMC).

References