Directivity Modes of Earthquake Populations with Unsupervised Learning

TL;DR

This paper introduces an unsupervised learning method using Gaussian mixture models to identify rupture directivity modes in earthquake populations, revealing distinct rupture directions and fault planes in large datasets.

Contribution

It presents a novel spectral decomposition and clustering approach to determine earthquake rupture modes without fitting kinematic models, improving understanding of rupture directions.

Findings

Earthquake datasets decompose into distinct rupture modes

Fault planes are unambiguously identified for all cases

Small earthquakes exhibit unilateral ruptures 53-74% of the time

Abstract

We present a novel approach for resolving modes of rupture directivity in large populations of earthquakes. A seismic spectral decomposition technique is used to first produce relative measurements of radiated energy for earthquakes in a spatially-compact cluster. The azimuthal distribution of energy for each earthquake is then assumed to result from one of several distinct modes of rupture propagation. Rather than fitting a kinematic rupture model to determine the most likely mode of rupture propagation, we instead treat the modes as latent variables and learn them with a Gaussian mixture model. The mixture model simultaneously determines the number of events that best identify with each mode. The technique is demonstrated on four datasets in California with several thousand earthquakes. We show that the datasets naturally decompose into distinct rupture propagation modes that correspond to different rupture directions, and the fault plane is unambiguously identified for all cases. We find that these small earthquakes exhibit unilateral ruptures 53-74% of the time on average. The results provide important observational constraints on the physics of earthquakes and faults.

Figures

Tables

Table 1. Table 1: Modal results for the Cahuilla swarm.

Mode	$w\left(K=2\right)$	$w\left(K=3\right)$
NW	0.55 [0.49, 0.51]	0.40 [0.30, 0.43]
SE	0.45 [0.41, 0.60]	0.34 [0.30, 0.45]
Bilateral		0.26 [0.21, 0.33]

Table 2. Table 2: Modal results for the Creeping section of San Andreas.

Mode	$w\left(K=2\right)$	$w\left(K=3\right)$
NW	0.50 [0.46, 0.56]	0.30 [0.26, 0.38]
SE	0.50 [0.44, 0.54]	0.24 [0.20, 0.40]
Bilateral		0.46 [0.26, 0.49]

Table 3. Table 3: Modal results for the Parkfield cluster.

Mode	$w\left(K=2\right)$	$w\left(K=3\right)$
NW	0.63 [0.56, 0.68]	0.21 [0.17, 0.42]
SE	0.37 [0.32, 0.44]	0.32 [0.28, 0.38]
Bilateral		0.47 [0.26, 0.51]

Table 4. Table 4: Modal results for the Hayward fault.

Mode	$w$	95% CI
NW	0.38 [0.33, 0.67]	0.37 [0.16, 0.47]
SE	0.62 [0.33, 0.67]	0.28 [0.17, 0.60]
Bilateral		0.35 [0.06, 0.48]

Equations

$$d_{ij}(f)=s_i(f)+p_{ij}(f)+s{t}_j(f)+{ϵ}_{ij}(f).$$

$$s{a}_{ij}(f)=d_{ij}(f)-p_{ij}(f)-s{t}_j(f).$$

$$E_R=4{\pi }^2{\int }_{f_{min}}^{f_{max}}{f}^2|sa(f){|}^{2}df.$$

$$p(x)=\sum _{k=1}^K{w}_k\mathcal{N}(x|{\mu }_k,{\Sigma }_k)$$

$$\sum _{k=1}^K{w}_k=1$$

Full text

1 Introduction

The kinematics of the earthquake source process has a first-order effect on seismic ground motions. In particular, it is well-known that for moving sources, rupture directivity induces azimuthally dependent changes in the frequency content of seismic waves (Haskell, 1964). Often earthquake ruptures are simplified as falling into one of two end-member modes: bilateral (symmetric) ruptures and fully unilateral ruptures. These rupture modes produce distinct signatures in the seismic radiation. Some studies discuss a continuum of scenarios in between these modes for which there is some degree of asymmetry in the distribution of seismic moment (e.g. Boatwright, 1984).

Theoretical studies of earthquake physics have suggested various links between rupture propagation and properties of fault zones. For example, large earthquakes that saturate the seismogenic zone are often viewed as being constrained geometrically to have elongated ruptures (McGuire et al., 2002). The bimaterial hypothesis (Weertman, 1980; Andrews & Ben-Zion, 1997) predicts that for faults with a velocity contrast across them, a rupture direction that is aligned with the direction of slip in the more compliant medium will be favored due to a dynamic reduction of normal stress at the crack tip. Other models may anticipate earthquakes to rupture with unilateral directivity as a frictional phenomenon, such as when the nucleation size is much smaller than the size of a seismogenic patch (e.g. Michel et al., 2017; Lin & Lapusta, 2018). It is desirable to better understand whether any of these models are relevant to earthquakes in nature, and whether there are resolvable differences in the properties of faults that can lead to different physics. Providing observational constraints on some of these physical models would require assembling directivity measurements for large numbers of earthquakes in order to perform a statistical analysis.

The extent of rupture directivity is readily determined for M $\geq 6$ earthquakes with various seismological methods (e.g. McGuire et al., 2001; Ye et al., 2016; Van Houtte & Denolle, 2018). These events are recorded typically by many stations and are sufficiently large that a slip model can be obtained using seismic waveforms recorded at teleseismic distances. However, the number of these large events on any single fault is too few to determine the statistical properties of rupture directivity patterns. Small earthquakes on the other hand are plentiful on individual faults due to the well-known power law magnitude distribution. There have been various approaches to estimating rupture directivity in small earthquakes, such as the use of higher-order moment tensors (McGuire, 2004), fitting directivity functions (Boatwright, 2007; Tan & Helmberger, 2010; Wang & Rubin, 2011; Kane et al., 2013; Abercrombie et al., 2017), measuring spectral splitting (Ross & Ben-Zion, 2016; Calderoni et al., 2015), measuring variations in apparent duration among repeating earthquake sequences (Lengliné & Got, 2011), and azimuthal analyses of ground motion (Kurzon et al., 2014). For most of these approaches, empirical Green’s functions are needed to correct the data for propagation and site effects (e.g. Mueller, 1985; Hough & Dreger, 1995; Prieto et al., 2004), and then some kind of geophysical inverse problem is solved to obtain measurements about directivity.

In this paper, we present a data-driven approach to learning modes of directivity in earthquake populations that can scale to large datasets (Fig. 1). We demonstrate that for a population of earthquakes, the azimuthal distribution of radiated energy naturally decomposes into distinct modes of rupture propagation without assuming a fault plane geometry or solving a geophysical inverse problem. The modes generally represent unilateral ruptures in different directions as well as bilateral ruptures. This decomposition simultaneously determines the fraction of events in each mode. We demonstrate the method on four datasets from California with thousands of events each, and clearly identify the fault plane and rupture directivity modes for all four cases. We show that the earthquakes in these clusters have predominantly unilateral ruptures, with bilateral ruptures being relatively infrequent.

2 Methods

Our approach to resolving directivity modes is based on the idea that an individual earthquake within a population can be represented as a random sample from one of a handful of end-member rupture scenarios. For strike-slip faults, the number of observable cases might be three: a bilateral rupture with symmetric rupture propagation, and two unilateral rupture modes that propagate in opposite directions along a fault. Alternatively, one might want to only consider unilateral directivity modes, with bilateral ruptures corresponding to a superposition of unilateral modes. The problem can then be formulated as one of recovering the distinct rupture modes that exist in the data. Our approach can be summarized with the following steps:

A cluster of earthquakes is identified for detailed analysis 2. 2.

The seismograms from these earthquakes are used with a seismic spectral decomposition algorithm to obtain apparent source spectra for as many events and stations as possible (Fig. 1, step 1). 3. 3.

The radiated energy at each station is determined by integrating the square of the velocity spectrum (Fig. 1, step 2). 4. 4.

We apply an algorithm for filling in missing data values by learning a low-rank approximation to the entire dataset. 5. 5.

Directivity modes are learned with a Gaussian mixture model. (Fig. 1, step 3) 6. 6.

Each mode can be visually interpreted as an azimuthal directivity function for a rupture mode that exists repeatedly in the data. The modes are obtained without solving a geophysical inverse problem.

We now discuss each of these steps in detail.

2.1 Data Pre-processing and Quality Control

The first step of the procedure involves selecting a cluster of earthquakes for analysis of rupture propagation. In this study, we use spatially-compact clusters for which the source-receiver azimuths are roughly the same for all events. We restrict the minimum source-receiver distance to 20 km to help ensure this approximation is valid. For each event, we obtain all available automated and manual P-wave picks, and for each associated waveform, we use these picks to define a 1.5 s signal window starting 0.1s before the listed P-wave arrival, and a background noise window of equal length immediately preceding it. We use a multitaper algorithm (Prieto et al., 2009; Krischer, 2016) to compute seismic amplitude spectra on all available vertical-component channels with sampling rates $\geq$100 Hz. We convert each spectrum to units of displacement and interpolate where necessary to a uniform frequency spacing of 2/3 Hz with minimum and maximum values of 0 Hz and 50 Hz (the Nyquist frequency for 100 Hz sampling).

For the purposes of quality control, we further consider only spectra with signal-to-noise ratio (SNR) greater than 5 within three frequency bands spanning the 3 Hz to 30 Hz range. Clipped waveforms are commonly observed on short-period and broadband stations for moderate earthquakes, and to mitigate this we further exclude spectra flagged by an automated clipped detection algorithm based on the fourth moment of the time domain waveform amplitudes observed in the signal window (Trugman & Shearer, 2017). The remaining displacement spectra after these waveform pre-processing steps and quality control steps are the inputs for the spectral decomposition method described below.

2.2 Spectral Decomposition

Spectral decomposition is a technique designed to separate the observed displacement spectra $d(f)$ of earthquakes $i$ recorded and a set of seismic stations $j$ into source, path, and site terms. As described in detail by Trugman & Shearer (2017), for densely-recorded datasets, each earthquake will be recorded by many stations, each approximate source-station path will be sampled many times, and each station will record many earthquakes. If these conditions hold, then by working in the log frequency domain, relative source terms $s_{i}$, path terms $p_{k(i,j)}$, and site terms $st_{j}$ can be estimated at each frequency point as part of an overdetermined inverse problem defined by the linear equation:

[TABLE]

Following Shearer et al. (2006) and Trugman & Shearer (2017), we assume azimuthally isotropic path terms $p$ that depend only on the source-receiver travel-time. We limit our analysis to stations with epicentral distances less than 100 km, and to ensure that the basic assumptions of spectral decomposition hold, we require each station used in the analysis to have recorded at least 20 earthquakes. We estimate the source, path, and site terms defined by equation (1) using an iterative, robust least-squares inversion that uses Huber-norm weighting to suppress the influence of outliers on the final solution (Shearer et al., 2006).

Our primary observational constraints come not from the source spectra $s(f)$ themselves, but instead the apparent source spectra $sa(f)$ observed at each station and how they vary as a function of azimuth:

[TABLE]

The apparent source spectrum $sa(f)$ is thus a path- and site-corrected form of the observed displacement spectrum, while the source spectrum $s(f)$ itself is the azimuthal average of apparent source spectrum across all stations (to good approximation). Rupture directivity can cause significant azimuthal variations in apparent spectra, with higher amplitudes along azimuths aligned with rupture direction. One general limitation of the spectral decomposition method is that the estimated terms only resolve relative differences between each source, each path, and each site, as one can add an arbitrary function to all source terms and subtract the same function from all site or path terms without effecting the misfit. While this presents a challenge for estimating earthquake source parameters like corner frequency and stress drop (Shearer et al., 2019), it does not affect the results presented in this study, which are based upon on the relative variations in spectral energy as a function of azimuth, and not the absolute energy values.

2.3 Building the Data Matrix

By this point, we have obtained apparent source spectra for many events at many stations. From here, we proceed to calculate apparent radiated energy values for the spectra. For an apparent source displacement spectrum, $sa(f)$, the apparent radiated energy (up a multiplicative constant) can be calculated as,

[TABLE]

The optimal values of $f_{min}$ and $f_{max}$ generally depend on the instrument response, sampling rate, and the magnitude range of the events for which the directivity analysis will be performed. They further depend on the signal to noise ratio of the spectra. Our spectra generally have low SNR below 3 Hz and above 30 Hz, so we set $f_{min}$ and $f_{max}$ equal to these values.

Next, we normalize the $E_{R}$ values separately for each event by dividing by the median value calculated over all of the stations. This results in a set of values that indicate whether $E_{R}$ at a given station is greater than or less the median. In doing so, we require a minimum of 10 stations for which $E_{R}$ values are available, or else the event is skipped altogether. We then remove all stations from the dataset which have fewer than 20 spectra over all the events. Then, we identify outlier values by calculating the median absolute deviation (in log units) and looking to see whether there are any values for a single event that are more than 5 deviations away from the median. If so, these values are removed from the dataset. If the outlier removal process results in fewer than 10 values for a single event, the event is skipped.

Applying the steps described in this section results in a set of relative $E_{R}$ values for many events and stations. However, of the complete set of stations available, few if any events will have $E_{R}$ values at all of them. This is especially true for the smallest events in the dataset, which are much more frequent than the larger ones. However, to fit a latent variable model to the data, there cannot be any missing values, and therefore they must be filled.

In applied mathematics there is an extensive literature on data imputation algorithms for filling missing data. The simplest form of these algorithms fills missing values with the mean or median for each variable. Other algorithms are more sophisticated, trying to learn a lower dimensional representation of the data from the observations that can be used to reconstruct the missing values. Algorithms of this type have seen some use in geophysics, for example to fill gaps in GPS data prior to performing principal component analysis or independent component analysis (e.g. Kositsky & Avouac, 2010). In this study, we chose an imputation algorithm that learns a low-rank approximation to the existing data, which can then be used to fill in missing values. This technique is based on thresholding a singular value decomposition of the data (Hastie et al., 2014). Some examples of applying this technique are shown in Figure 2, where real values are shown alongside imputed values. Filled values have a tendency to be generally close to zero (the mean) unless there is significant evidence in the data otherwise.

2.4 Gaussian mixture models

Gaussian mixture models (GMM) are a type of latent variable model that assumes a generative process can be represented as a superposition of two or more Gaussian distributions. Each Gaussian represents a distinct mode of occurrence, and there is a probability $w_{i}$ that a realization $x\in\mathbb{R}^{d}$ will be drawn from mode $i$. Specifically, the model can be written as,

[TABLE]

[TABLE]

where $K$ is the number of components in the GMM and $\mathcal{N}(x|\mu_{k},\Sigma_{k})$ is the multivariate Gaussian distribution. Since $\mu_{k}$ and $\Sigma_{k}$ are not directly observed, they are latent variables; however they can be recovered by fitting a GMM to some data. The most common approach to estimating the parameters of a GMM is with the Expectation-Maximization algorithm (Dempster et al., 1977), which is the method we use in this paper. We also tested the method of moments algorithm (Anandkumar et al., 2014) that matches the empirical moments to the theoretical ones through a tensor decomposition algorithm. We found that the EM algorithm was more effective at separating the modes for the datasets considered in this study.

Our measured $\log{E_{R}}$ values for each earthquake correspond to the $x$ in the GMM, and we therefore wish to obtain the $\mu_{k}$ and $\Sigma_{k}$ for the $K$ modes in the data. For this paper, we will work with strike-slip faults and consider both $K=2$ and $K=3$ cases. The $K=2$ scenario would ideally recover two unilateral directivity modes, while the $K=3$ scenario would correspond to a bilateral rupture mode and two opposite unilateral modes. We will further assume that the covariances are spherical, i.e. $\Sigma_{k}=\sigma_{k}^{2}I$. With this formulation, each $\mu_{k}$ is a vector with dimensionality equal to the number of stations used in the analysis, and the measured $\log{E_{R}}$ of an individual earthquake are modeled as a random sample from the GMM.

2.5 Uncertainty Analysis

To estimate the uncertainty in the model parameters, we use a bootstrapping approach. We resample the events with replacement 1000 times and re-fit the model each time. Since the order of the elements of $w$ is randomly determined when the model is fit, we obtain all permutations of the columns of $\mu$ and calculate the $\ell_{1}$ norm of the difference between the best fitting $\mu^{\text{best}}$ and each of the permutations. The permutation with the smallest $\ell_{1}$ value is taken as the $w$ for that particular resampling. Then we calculate a 95% confidence interval (CI) for each of the $w_{k}$.

3 Related Work

In this study, our goal is to estimate the fraction of events for each rupture directivity mode. We are interested in working with hundreds to thousands of earthquakes at a time to ensure that the statistics are robust. To date, there are several studies that have examined directivity on such a scale.

Wang & Rubin (2011) calculated spectral ratios using empirical Green’s functions and fit three different models (bilateral and two unilateral) to each event separately to determine the most likely rupture mode. They perform this analysis in the creeping section of the San Andreas fault and test the method on more than 900 earthquakes. They concluded that roughly 40% had bilateral ruptures, and of the unilateral ruptures, most had southeast directivity.

Kane et al. (2013) used a spectral decomposition technique to obtain apparent source spectra for thousands of events at the Parkfield section of the San Andreas fault. They fit a unilateral directivity model to the apparent spectra for each event individually, and concluded that there was slightly more events with southeast ruptures than northwest ruptures. They did not consider bilateral ruptures.

Building on the results of Calderoni et al. (2015) and Pacor et al. (2016) for the L’Aquila, Italy earthquake sequence, Calderoni et al. (2017) used azimuthal variations in S-wave spectra to demonstrate that along-strike directivity is a common feature of normal faulting earthquakes in the central Appennines. Earthquakes from the Umbria Marche, L’Aquila, and 2016-2017 central Italy sequences exhibit temporally persistent and spatially coherent directivity patterns, which suggests that in-situ fault and geologic properties play an important role in the preferred rupture direction.

A different approach was developed by Lengliné & Got (2011), who exploited repeating earthquake sequences in Parkfield to simultaneously invert for relative perturbations to apparent duration for all earthquake pairs and stations. They used P-wave signals for their analysis, and their kinematic directivity model only allowed for unilateral ruptures. Lengliné & Got (2011) concluded that the majority of earthquakes analyzed exhibited southeast rupture signals.

A time-domain method based on ground-motion prediction equations was developed by Kurzon et al. (2014). They proposed an empirical directivity index based on observed azimuthal variations in ground motion amplitude. They applied the method to two clusters of earthquakes in the central San Jacinto fault zone.

4 Experiments

Now, we demonstrate the method on four different datasets from California. These datasets use P-wave seismograms only and were chosen for several reasons. First, each dataset is a spatially-compact cluster with thousands of events. Second, two of them have been used in previous studies to analyze directivity and can serve as a point of reference.

4.1 Cahuilla swarm

The first dataset contains seismograms for 11,631 events that occurred as part of the Cahuilla earthquake swarm in Southern California during 2016-2019 (Hauksson et al. (2019); Fig. 3). The largest of these events is $M4.4$, while most events are considerably smaller. These events were studied in detail by Hauksson et al. (2019), who observed that the events have right-lateral strike-slip focal mechanisms with generally little variation in strike. The waveform and meta data are all publicly available from the Southern California Earthquake Data Center.

The results of applying the method to the Cahuilla swarm data are shown in Figure 4 and Table 1. In total, there are 829 earthquakes and 86 stations used for the final decomposition. First, we show the results for $K=2$. In this plot, each mode is assigned a different color. The elements of $\mu_{k}$ represent the average value of $\log{E_{R}}$ at a given station for mode $k$. Since there are two modes, each station has two different colored circles present. Value larger than zero indicate that $\log{E_{R}}$ at that station are amplified relative to the median (for that mode). $\mu_{1}$ and $\mu_{2}$ have peaks along the strike of the fault but in opposite directions. There is more than a factor of 15 difference in the $E_{R}$ between these opposing directions. Thus, the method has been able to identify the fault plane, and these modes represent unilateral directivity. For these two modes in the Cahuilla swarm (Table 1), $w_{1}=0.55$ and $w_{2}=0.45$. Since $w_{1}\geq w_{2}$ in 95% of the bootstrap samples, ruptures are statistically more likely to have NW directivity.

For the $K=3$ decomposition (Fig. 5), two of the modes look very similar to the $K=2$ results. The third mode exhibits more complex behavior, with four peaks and four troughs. These peaks are essentially at conjugate orientations and suggest that there may be cross-faulting mixed in. Events have bilateral ruptures between 21-33% of the time (Table 1). Thus it is the least likely of the three modes to occur.

4.2 San Andreas (Creeping section)

The second dataset analyzed in this paper is for the creeping section of the San Andreas fault. This area has a very active seismicity cluster that has produced 4118 events from 2002-2019 (Fig. 6). All of the data are publicly available from the Northern California Earthquake Data Center. This cluster was studied by Wang & Rubin (2011) and serves as an additional baseline for our work.

The $K=2$ results for the creeping section of the San Andreas fault are shown in Figure 7 and Table 2. A total of 969 earthquakes and 90 stations were used to fit the model. As with the previous dataset, the fault plane is clearly identified for this segment, with prominent peaks in the NW azimuths for one mode and the SE azimuths for the other mode. The modes are evenly distributed on average and the confidence intervals further show that there is no evidence of a statistically preferred direction.

For the $K=3$ decomposition, there are two unilateral modes (NW and SE) and one bilateral mode (Fig. 8) the bilateral mode is the most frequent ($w=0.46$), and after factoring in the uncertainty, ruptures are unilateral 51-74% of the time. Wang & Rubin (2011) found that about 40% of the earthquakes on this segment had bilateral ruptures, which is similar to our observations. They also found that of the unilateral ruptures, SE ruptures were more common; however we find no evidence of this that is statistically significant.

4.3 San Andreas (Parkfield)

The third dataset is for the Parkfield section of the San Andreas fault in central California (Fig. 9). We selected all 14,562 events in the NCEDC catalog that have occurred since 2002. This region was used for several directivity studies in the past (Kane et al., 2013; Lengliné & Got, 2011) and serves as a point of comparison for our results. The earthquakes at Parkfield have very homogeneous focal mechanisms (Thurber et al., 2006) which help to simplify the demonstrations of our approach.

Our results for the Parkfield data $(K=2)$ are shown in Figure 10 and Table 3. The model was fit to 618 events at 76 stations. The strike and fault plane are clearly identified, with prominent peaks in in the modes along the NW and SE azimuths. For the NW mode, $w=0.63$ with a one-tailed confidence interval of 0.56 at the 95% level, which indicates that NW ruptures are statistically more frequent than SE ruptures. This finding contrasts the results of Lengliné & Got (2011), who concluded that most of the events had SE rupture directivity signals.

In Figure 11, the results for the $K=3$ decomposition are shown. For this case, the NW and SE modes do not decompose as cleanly as in the $K=2$ case. This mixing makes it more difficult to confidently interpret the results, compared with the $K=2$ results. We find that the mode that appears the most bilateral occurs 26-51% of the time.

4.4 Hayward fault

The final dataset analyzed in this paper is for the Hayward fault. The seismicity here occurs over a long stretch of the fault, and so we chose a compact cluster to work with to ensure that the azimuths were similar to all events. There were 4118 events from 2002-2019 over the whole area and we used these to determine the source spectra (Fig. 12). All of the data are publicly available from the Northern California Earthquake Data Center. The Hayward fault has a velocity contrast in the range 3-8% (Allam et al., 2014).

After all of the pre-processing, the Hayward dataset has 170 earthquakes at 82 stations. We show the results for $K=2$ in Figure 13 and Table 4. As with the previous study areas, here the data separate clearly into two modes with directivity behavior in NW and SE directions, which aligns with the strike of the fault. The $w$ values for the modes vary significantly for the best-fitting values, but interestingly, have exactly the same confidence intervals. Thus, they are statistically indistinguishable, albeit with large uncertainties.

When applying a $K=3$ model to the data (Fig. 14, the modes are distinct and form smoothly varying patterns as a function of azimuth. The three modes also have sizable uncertainties as estimated by the bootstrapping, but we can say that unilateral modes occur between 52-94% of the time; thus they are the most frequent mode of rupture.

5 Discussion

5.1 The importance of studying directivity modes

Directivity modes characterize the first-order kinematics of rupture propagation during earthquakes. It is important to understand these modes and their statistical tendencies because they can provide valuable observational constraints on the physics of earthquakes. One of the important findings of this study is that for the four examined fault segments, unilateral ruptures are more frequent than bilateral ruptures. Specifically, the best-fitting rates vary from 53-74% over the four regions tested, and the uncertainty estimates show that this principal conclusion is robust. We find evidence of a statistically preferred rupture direction in two regions: the Cahuilla swarm and Parkfield. For the other regions, the uncertainties on the $w_{i}$ values are larger than the differences between them. This does not mean a preferred direction does not exist; rather the question could potentially be addressed better in the future by incorporating more events, which may tighten the confidence intervals. In all cases, however, there is no evidence of a single unilateral rupture direction being overwhelmingly likely.

Another noteworthy finding is that for these fault segments, there is a clear preference for ruptures to propagate horizontally more often than vertically. As discussed below, along-dip ruptures will generally be included in the bilateral mode because there is little sensitivity to them. This means that our conclusions about the predominance of unilateral ruptures not only applies with respect to bilateral cases, but also along-dip cases. Such a finding was made by McGuire et al. (2002), who observed that most of these were unilateral ruptures. For earthquakes large enough to rupture the full seismogenic zone, along-strike ruptures are expected purely from a geometrical perspective, although this says nothing about whether they should be unilateral or bilateral. However for small earthquakes, our observations may be unexpected. They further suggest that there are attributes of the faults that are breaking the symmetry and leading to horizontal ruptures being so prevalent.

An additional reason that these findings are important is because observational studies of earthquake source properties commonly assume that rupture areas are circular (e.g. Calderoni et al., 2012; Abercrombie, 2015; Ross & Ben-Zion, 2016). However, the geometry of the rupture area has an important influence on quantities like stress drop and radiated energy (Kaneko & Shearer, 2015). Our results suggest that for the examined datasets, unilateral ruptures are more frequent than bilateral ruptures, and this may indicate that some of the large scatter commonly observed in source properties may result from assuming the wrong rupture geometry.

One physical explanation for the systematic tendencies for unilateral directivity in earthquakes is from rate- and state-dependent friction. When the nucleation size $h$ is much smaller than the potential seismogenic region, this results in ruptures that tend to propagate unilaterally (e.g. Michel et al., 2017; Lin & Lapusta, 2018). However the expected rates of unilateral ruptures for this type of physical model have not been analyzed rigorously. Future work in this area could be helpful to compare with observations.

Directivity modes are also important from a hazard perspective because the rupture propagation pattern has a significant effect on the strong ground motion. While earthquake stress drop is widely accepted to influence ground motion amplitudes in certain contexts (e.g. Trugman & Shearer, 2018; Oth et al., 2017; Baltay et al., 2017), directivity has received less attention to date but is readily acknowledged to be one of the few additional source parameters that does have a major impact (Bozorgnia et al., 2014; Douglas & Edwards, 2016). Therefore being able to provide some expectations on the possible range of modal probabilities can better inform hazard estimates.

5.2 Data-driven approaches to studying directivity

There are a variety of advantages to studying directivity with a data-driven approach, as opposed to using techniques that fit kinematic rupture models. The first is that a single model can be fit to the data for all earthquakes simultaneously. This allows for weaker signals to be extracted from the data because the information that is common to all events can be averaged, which acts to suppress noise. Traditional approaches generally fit a model separately to each event, and then if performed for enough events, then the results can be averaged for some statistical estimate. However this type of approach is limited by the capabilities of fitting a model to individual events.

Another reason is that there are generally fewer assumptions necessary with a data-driven approach, compared with traditional model fitting. For example, we do not need to assume that our physical model is correct, or that we understand the statistical properties of the noise. However, as in most analyses of earthquake directivity, we assume that azimuthal variations in seismic spectra are caused primarily by source effects, rather than 3D variations in path effects. While this is likely a valid approximation to first order, in reality the observed spectra will contain hints of both effects.

5.3 Assumptions and interpretability of results

In this study, there are several key assumptions that underlie the analysis. The first, and arguably the most important, is that we assumed the number of distinct modes of rupture propagation a priori. By making an assumption that a certain number of modes exist in the data, we are able to search for and identify their centroids $\mu_{k}$ and variances. If these assumptions are violated, for example in the case that there are more than three distinct modes, then the resulting $\mu_{k}$ vectors will be complicated and uninterpretable. In all four datasets that we tested here, the $\mu_{k}$ are generally simple, smoothly varying functions of azimuth, which suggests that the assumptions are justified. However when applying the data to other datasets, such as clusters with more than one dominant fault strike, more care will be needed to determine the optimal number of modes in the data. There are various strategies for determining this objectively, including the use of information criteria or silhouette analysis. We applied a silhouette analysis to the data and in all regions, the metric favored $K=2$. However, as we have shown, $K=3$ models provide additional insights into the directivity physics and we believe that examining both scenarios together provides a more informed result.

If the recovered modes reflect approximately end-member rupture scenarios, one way of assessing the degree of unilateral directivity for a given event (e.g. Boatwright, 2007) is by calculating the probability that an event belongs to the best-matching mode. This could be useful for identifying the events with the strongest directivity signals in an objective manner. However, it should be mentioned that this strategy will only determine whether an event best matches the centroid of the mode; if an event has even stronger directivity signals than the centroid, the probability will be diminished in proportion to the deviation from the centroid. This simply means that events with the highest probability values will not be the events with the strongest directivity.

Another important aspect of the results is that the modal shapes which we have identified as bilateral are a bit more complicated in practice. In particular, if a rupture is purely along-dip, then without stations right on top of the rupture, it is very difficult to observe any directivity patterns. These ruptures will appear as generally symmetric directivity patterns and will likely be aggregated into the bilateral mode. Thus, the bilateral mode (particularly the $w_{k}$ value associated with it) should be interpreted as an upper bound for the likelihood of bilateral ruptures. It should be understood that this results from a lack of sensitivity to along-dip ruptures. Some studies have included other information, such as depth phases (e.g. He & Ni, 2017), to help constrain these cases, but in general it is a challenge common to all directivity analyses.

The results in the paper do not depend on the particular model used (e.g. GMM). We tested several other latent variable models including BIRCH and spectral clustering, and when setting the number of clusters to $K=2,3$, recover modes that look very similar.

6 Conclusions

We developed a new approach to resolving modes of directivity in large earthquake populations. We formulate the problem as one of recovering $K$ latent variable modes from the azimuthal energy distributions of many earthquakes, where each mode is a distinct state of rupture propagation. A gaussian mixture model is used to determine these modes, which allows for simultaneous estimation of the fraction of events that best align with each mode of rupture propagation. In the process, we have not needed to fit a kinematic directivity model to the data; the decomposition is possible because the data exhibit this type of structure naturally. We applied the method to four large earthquake clusters, and performed an assessment of the uncertainties for each mode of rupture, for each dataset. Our results indicate that unilateral ruptures are more likely to occur than bilateral ruptures for the examined datasets, even after incorporating the uncertainties.

Acknowledgments

A. Anandkumar is supported in part by Bren endowed chair, Darpa PAI, Raytheon, and Microsoft, Google and Adobe faculty fellowships