Clustering · Part B: the Geysers catalogue

The Geysers has recorded tens of thousands of earthquakes with focal mechanisms. Do they fall into groups — in space, in mechanism, in time — and which of those groups are real?

Part A taught the methods on planted examples, where every answer could be scored against the truth. Here there is no truth. We take the whole focal-mechanism catalogue of the Geysers, every event with a first-motion solution from 1969 to 2026, and ask the questions Part A prepared: what counts as close, what counts as a group, at what scale — and is it real?

You can run it without running Part A first.

The questions this notebook answers

1. Where are the earthquakes, and do they form groups at all? (B.1)
2. What do the five methods of Part A find on these data, and what does each one's parameter do? (B.2–B.3)
3. Which groups are real: do they survive a change of scale, the location errors, and a comparison with a catalogue that has no small-scale structure? (B.4–B.5)
4. Do the focal mechanisms form groups of their own, and do groups found from position alone have their own mechanism? (B.6–B.7)
5. In time, which earthquakes were triggered by an earlier one, and are the Geysers bursts aftershock sequences or swarms? (B.8–B.10)

B.1 · The catalogue, and where the earthquakes are

The focal-mechanism catalogue lists, for every earthquake that FPFIT (Reasenberg & Oppenheimer 1985) could solve, its origin time, its hypocentre, its magnitude and its mechanism as strike, dip and rake. FPFIT searches for the mechanism whose predicted first motions (up or down at each station) best match the observed ones; its misfit is the weighted fraction of first motions the mechanism gets wrong. FPFIT sometimes reports more than one solution for the same earthquake; we keep the one with the smallest misfit, so each earthquake appears once. We apply no other cut: every magnitude, every year, every quality — including the solutions whose misfit is above 0.2, which a stricter analysis would drop. Positions are converted to kilometres east and north of the field centre, and depth is in kilometres.

46,395 earthquakes with a focal mechanism, 1969–2026
magnitude: median 1.46, range 0.0–5.7
first motions per solution: at least 25, median 36; misfit above 0.2: 18% of the events
from 2005 on: 92% of the events; before 2003: at most 184 events a year

Most of the catalogue is recent: before 2003 no year has more than 184 events with a mechanism.

To look inside the cloud we need a depth section as well as a map. The field is elongated from north-west to south-east, so the section line A–A′ runs along its long axis: the direction in which the epicentres spread most (the eigenvector of their covariance matrix with the largest eigenvalue). Every event is projected onto that line. Map and section are drawn at the same scale — a kilometre is the same length on both — so a cluster keeps its true shape. Rather than one dot per event, which would hide the densest places under overlapping dots, each 100 m cell is coloured by how many earthquakes it holds.

80% of the events are between 0.7 and 3.2 km deep (median 1.9 km)

The earthquakes are not spread evenly through a cloud. They concentrate in compact patches, each about a kilometre across, separated by quieter rock; four in five events are between 0.7 and 3.2 km deep. In the section, patches at different distances across the field are projected onto the same line and overlap; the map shows them apart. This is what clustering should find — and the rest of this notebook asks whether it does, and whether what it finds is real.

B.2 · Two definitions of a cluster on the same earthquakes

Part A's first lesson was that a clustering method decides what a group is. Here are two of its definitions side by side, on the same 3-D hypocentres (east, north, depth, all in kilometres, so no feature needs rescaling):

• k-means (A.2): a cluster is the set of events nearest to one centre, and every event belongs to some cluster;
• HDBSCAN (A.6): a cluster is a region denser than its surroundings, and events in sparse places belong to no cluster. We use cluster_selection_method="leaf", which keeps the finest dense patches, with min_cluster_size=50 and min_samples=10.

To compare the definitions rather than the numbers, k-means is given the number of clusters HDBSCAN finds.

HDBSCAN: 87 clusters, 31% of the events in a cluster
k-means with k = 87: every event in a cluster

k-means cuts the whole field into tiles of similar size, sparse fringes included, and draws its boundaries straight through the dense patches of B.1 wherever the tiles happen to meet. HDBSCAN keeps those patches — each coloured cluster on its map sits on a dark patch of B.1 — and leaves the sparse events between them in no cluster (grey). Neither number of clusters came from the data: it came from min_cluster_size. The next section looks at the parameter of each method in turn.

B.3 · Every method has a parameter the data do not settle

Part A (A.7) showed each method's parameter on a map where the truth was known, and the diagnostic one would use to choose it. Here are the same five methods on the Geysers. Each figure shows four settings of the parameter that matters most, each as a map with its section under it, and on the right the diagnostic. The settings are chosen to bracket the value the diagnostic recommends.

k:           [2, 3, 5, 10, 20, 50, 100]
silhouette:  [0.54 0.46 0.5  0.39 0.32 0.32 0.32]
two halves:  [0.82 0.45 0.83 0.8  0.65 0.78 0.69]

k-means. The silhouette is highest at k = 2: the split of the field into a north-western and a south-eastern part, and no larger k is better separated. Larger k only tiles the cloud more finely. The two random halves agree about equally well on k = 2 and k = 5, and much less on k = 3; no larger k is clearly more reproducible. The catalogue reaches tens of kilometres beyond the field, and k-means spends clusters on those few distant events: at k = 5 only two clusters are inside this map.

BIC (spherical), k = [2, 5, 10, 20, 50, 100, 150]: lowest at k = 150
BIC (full), k = [2, 5, 10, 20, 50, 100, 150]: lowest at k = 150

Gaussian mixture. The covariance shape decides what a cluster can look like: spherical clusters are round, like soft k-means; full clusters are ellipsoids of any orientation and can follow an elongated patch. BIC is meant to choose k, and on Part A's map it did. Here it keeps falling up to the largest k tried, for both shapes: it never chooses, which says the cloud is not a few Gaussian blobs. A diagnostic works only when the data match its assumptions.

Ward: the largest relative jump in merge distance among the first cuts is at 3 clusters

single linkage, 20 clusters: the largest holds 100% of the events
average linkage, 20 clusters: the largest holds 100% of the events

Hierarchical clustering. Each event is joined only to its ten nearest neighbours, which keeps the tree computable (A.4). The linkage matters more than where the tree is cut. Single linkage joins two groups through their closest pair of events, so one chain of events joins almost the whole field into one cluster; average linkage does nearly the same. Ward linkage joins the pair that least increases the spread inside clusters and gives compact clusters. Cut where the merge distance jumps, it returns three clusters: the same coarse north-west/south-east split as k-means, and a third made of the distant events off this map.

k-distance knee (10th neighbour): 1.09 km

eps 0.1 km, min_samples 10: 160 clusters, the largest holds 11% of the events
eps 0.2 km, min_samples 10: 68 clusters, the largest holds 48% of the events
eps 0.2 km, min_samples 80: 13 clusters, the largest holds 8% of the events
eps 1.09 km, min_samples 10: 36 clusters, the largest holds 93% of the events

DBSCAN. eps and min_samples act together as one density threshold (A.5): 10 events within 0.1 km and 80 within 0.2 km are the same density in 3-D, eight times the volume and eight times the events. The usual rule reads eps from the knee of the k-distance plot; here the knee is near a kilometre and joins almost the whole field into one cluster. On Part A's map the knee landed in DBSCAN's narrow working window; here the density varies by orders of magnitude across the field, and no single threshold suits it — a small eps finds only the densest patches, a larger one merges the dense centre into one cluster.

excess of mass, size 50: 15 clusters, the largest holds 93% of the events
leaf, size 50: 87 clusters, the largest holds 1% of the events

HDBSCAN. It builds DBSCAN's clusters at every density at once and then selects from that tree (A.6). The selection method matters first: "excess of mass" keeps the most persistent branches, which here puts most of the field in one cluster; "leaf" keeps the finest dense patches. With leaf selection, min_cluster_size — the smallest group we would call a cluster — then sets the scale, which is the subject of B.4.

Together the five figures say what A.7 said on the planted map, with one difference: there the diagnostics worked, and here most of them do not. They find only the coarse north-west/south-east split (k-means, Ward) or fail outright (BIC, the k-distance knee). The one parameter with a physical meaning is HDBSCAN's min_cluster_size, and the next section asks what answer survives a change of it.

B.4 · How many clusters? The answer depends on the scale — but some patches do not

HDBSCAN's min_cluster_size is the smallest group we would call a cluster. We run it at four sizes and ask which events are in a cluster at every one of them. Those events form the robust patches: they do not depend on our choice of scale.

min_cluster_size  25: 135 clusters, 29% of the events in a cluster
min_cluster_size  50:  87 clusters, 31% of the events in a cluster
min_cluster_size 100:  53 clusters, 33% of the events in a cluster
min_cluster_size 200:  34 clusters, 33% of the events in a cluster
20% of the events are in a cluster at all four sizes; they form 34 robust patches

From 25 to 200 the number of clusters falls by a factor of four while the share of events in clusters barely changes: a larger min_cluster_size merges small patches rather than adding new events. About a fifth of the events are in a cluster at every size, and each of the large clusters contains such a core. These cores — the robust patches — are what we will interpret; a cluster found at only one size should not be.

B.5 · Are the patches real?

Part A (A.8) showed that every method returns clusters even from random events, so a cluster count proves nothing by itself: it must be compared with a catalogue that lacks the structure in question. We keep the large-scale shape of the seismicity and move the events instead, in two ways that ask two different questions:

• each event moved by its own location error. The catalogue gives each hypocentre a horizontal and a vertical error; we move every event by a random amount of that size (in each horizontal direction by the horizontal error divided by $\sqrt{2}$, vertically by the vertical error). This is a catalogue we could equally have recorded, and it asks: do the patches survive the uncertainty of the locations?
• each event moved by about half a kilometre in each direction. This keeps the outline of the field but erases structure smaller than about a kilometre, and asks: are the patches more than chance lumps in a smooth cloud?

The errors are in the location catalogue, joined to the mechanism catalogue by event id.

location errors found for 98.5% of the events; median 0.15 km horizontal, 0.25 km vertical

             data: clusters at min_cluster_size [25, 50, 100, 200] = [135, 87, 53, 34]
   location error: clusters at min_cluster_size [25, 50, 100, 200] = [104, 63, 40, 24]
 half a kilometre: clusters at min_cluster_size [25, 50, 100, 200] = [83, 33, 19, 8]

Moved by its location error, the catalogue keeps most of its clusters: the patches are much larger than the location errors and survive them. Moved by half a kilometre, it keeps far fewer, and the gap widens with the cluster size — at the smallest size every catalogue gives many clusters, which is why a count at one setting proves nothing, while at the largest size the data have several times as many clusters as the blurred catalogue. The large, persistent patches are more than chance lumps.

B.6 · Do the focal mechanisms form groups?

Now the feature is the mechanism. A mechanism is reported as strike, dip and rake, but these three angles are a poor feature for clustering: the same double couple has two descriptions, one for each nodal plane, and they can be far apart in strike, dip and rake. Session 3, Focal mechanism & moment tensor, wrote a double couple as a moment tensor, six numbers that are the same whichever plane describes it. We use the tensor scaled to unit size, so that the distance between two mechanisms depends only on their orientation: two identical mechanisms are at distance 0 and two opposite ones at distance 2.

To name a mechanism's style we use its three axes: the P (pressure) and T (tension) axes, and the B (null) axis perpendicular to both. A mechanism is normal if its P axis plunges more steeply than 60°, reverse if its T axis does, strike-slip if its B axis does, and oblique if none of the three is that steep.

FPFIT searches a grid of mechanisms, and reports the best grid point: strike and dip in steps of 5°, rake in steps of 10°. Many earthquakes therefore share exactly the same reported mechanism. A clustering method sees those repeated values as dense points.

17,641 distinct reported mechanisms among 46,395 events; 84% of the events share theirs with another event

each cluster of the reported mechanisms spans a median of 27 grid points
HDBSCAN on the mechanisms as reported: 93 clusters, 42% assigned; largest 4,900 events
HDBSCAN on the mechanisms spread within the grid cell: 2 clusters, 69% assigned; largest 32,106 events

On the reported values, HDBSCAN finds dozens of clusters. Each groups a few dozen neighbouring grid points: the repeated values make the reported mechanisms lumpy, and the method cuts between the lumps. Spread each solution within its grid cell, which is as much as the grid tells us, and only a few clusters remain, one of them holding most of the events: the dozens of clusters were structure in the method that produced the catalogue, not in the Earth. Before clustering a feature, find out how it was produced.

Is that region one group, or several that touch? k-means with k = 2 splits it, and Part A's stability check asks whether two random halves of the catalogue agree on the split.

k:                   [2, 3, 4, 5, 6, 8, 10]
two halves, Geysers: [0.87 0.59 0.81 0.81 0.5  0.46 0.46]
two halves, random:  [0.05 0.13 0.16 0.2  0.24 0.2  0.22]

The two halves agree on the k = 2 split much better than they would on random mechanisms, which have no preferred split at all. But the right panel shows what the split is: projected onto the line between the two centres, the events form one hump with no gap, and k-means cuts it in the middle. The split is reproducible because the cloud of mechanisms is elongated — it runs from normal faulting at one end toward oblique and strike-slip faulting at the other — and k-means cuts any elongated cloud across its long axis, in the same place every time. There are not two separate families. This is A.8's lesson on real data: a reproducible clustering is not proof that clusters exist. The mechanisms form one continuum, from normal to oblique and strike-slip faulting.

B.7 · Do the patches have their own mechanism?

The robust patches of B.4 were found from position alone, so the mechanism is a feature the clustering never saw. Part A (A.8) used such a held-out feature to test whether clusters mean something. If a patch is a real structure — one fault zone, one volume under similar stress — its earthquakes should have more similar mechanisms than a random group of events.

We measure how similar the mechanisms in a group are by the length of their average unit tensor: it is 1 when they are all the same and near 0 when they point every way. Then we compare each patch with groups of the same size drawn not from the whole field but from its own neighbours: events between 0.5 and 2 km from the patch that belong to no robust patch. Neighbours share the patch's regional stress field, and they are recorded by the same stations along nearly the same ray paths, so they share its errors in FPFIT too. Part A's toy compared with random groups from the whole map; on real data the local comparison is the fairer one, because a patch could beat the whole field simply by sitting in a region where every mechanism is alike. A patch passes if it is more uniform than 95% of the neighbour groups. We also resample each patch's own events 200 times and keep its style only if the dominant mechanism's style is the same in at least 80% of the resamples.

31 of 34 patches have enough neighbours to compare
more uniform than their neighbours: 20; and their style survives resampling: 18
their styles: {'strike-slip': 9, 'normal': 6, 'oblique': 3}
uniformity, patches with their own mechanism: median 0.47; their neighbours: median 0.35

More than half of the patches have a mechanism of their own: more uniform than their neighbours, with a style that survives resampling. The effect is real but modest — the typical uniformity inside such a patch is only somewhat above that of its neighbours, so many events in a patch still differ from its dominant mechanism. On the map, normal faulting dominates the patches along the northern edge of the field, and strike-slip those along its axis toward the south-east. The patches are structures in the Earth, not only in the catalogue: a feature the clustering never saw agrees with them.

B.8 · Space, time and magnitude: the parameters, measured

So far every event was a point in space or in mechanism. Now origin time joins, and with it the question Part A (A.9–A.10) asked: which earthquakes were triggered by an earlier one? The nearest-neighbour method of Zaliapin & Ben-Zion (2013a) links each event to the earlier event with the smallest nearest-neighbour proximity

$$\eta = \Delta t \; r^{d_f} \; 10^{-b\,m},$$

where $\Delta t$ is the time since the earlier event (in years), $r$ the distance between their epicentres (in kilometres), $m$ the earlier event's magnitude, $d_f$ the fractal dimension of the epicentres and $b$ the Gutenberg–Richter b-value. In Part A both were planted; here they are measured.

• The fractal dimension. The fraction of pairs of epicentres closer than $r$ grows as $r^{d_f}$, so $d_f$ is the slope of that fraction against $r$ on logarithmic axes.
• The b-value. The maximum-likelihood estimate from Session 1, Regression & uncertainty, above a cutoff magnitude. It needs a catalogue that is complete above the cutoff: one that contains every earthquake above it.

cutoff M1.0: b = 0.68 (focal-mechanism catalogue), 1.02 (location catalogue)
cutoff M1.4: b = 0.82 (focal-mechanism catalogue), 1.06 (location catalogue)
cutoff M2.0: b = 0.96 (focal-mechanism catalogue), 1.12 (location catalogue)
cutoff M2.8: b = 1.24 (focal-mechanism catalogue), 1.18 (location catalogue)

fractal dimension of the epicentres d_f = 1.60; b used = 1.06

The fractal dimension is about 1.6: the epicentres fill more than a line and less than the plane. The b-value needs care. For the location catalogue, b stays between 1.0 and 1.2 whatever the cutoff. For the focal-mechanism catalogue it keeps rising as the cutoff rises, because that catalogue is not complete. A mechanism needs many clear first motions (every solution here has at least 25), and a small earthquake is recorded clearly at fewer stations. So the smaller the earthquake, the more likely it is to have no mechanism — not only near the smallest magnitudes, but well into the range where the catalogue looks complete. (Both catalogues mix magnitude types: most events have duration magnitudes, the largest have local or moment magnitudes.) The b-value is a property of the earthquakes, not of which ones got a mechanism, so we use the location catalogue's value. Part A's lesson about features applies to parameters too: find out how a number was produced before using it.

B.9 · Which earthquakes were triggered by an earlier one?

Each event's parent is the earlier event with the smallest $\eta$. Splitting $\eta$ into a rescaled time $T = \Delta t\,10^{-b m/2}$ and a rescaled distance $R = r^{d_f}\,10^{-b m/2}$, so that $\eta = T R$, shows where a link sits. The threshold $\eta_0$ separates strong links (triggered) from weak ones (background). Zaliapin & Ben-Zion (2013a) chose $\eta_0 = 10^{-5}$ by inspecting the two modes of their catalogue. We use the rule from Part A instead: a null catalogue keeps every event's position and magnitude and the number of events in each month, but moves each event's time to a random moment within its month, which destroys triggering while keeping the slow changes of the rate; $\eta_0$ is set where the null has 5% as many strong links as the data.

log10 eta0 = -7.56: 7.6% of the events have a strong link (null: 0.38%)
at Zaliapin & Ben-Zion's eta0 = 1e-5: data 62%, null 60%

Most of the Geysers catalogue sits in one cloud that the null reproduces: background events, no closer to an earlier event than chance allows. The clustered population is a weak lobe at small rescaled time, below the threshold line — far weaker than in Part A's synthetic catalogue. Only a small minority of the events — fewer than one in ten — is strongly linked to an earlier one. A fixed threshold would not work here: at Zaliapin & Ben-Zion's $\eta_0 = 10^{-5}$ the null itself has almost as many strong links as the data, because $\eta$ depends on the units and the event density of each catalogue. The threshold has to come from the null.

What do the strong links join? A year of the catalogue, drawn as distance along A–A′ against time, shows them. Events joined by strong links form families; each family is one tree, rooted at its first event.

2016: 1,828 events, 267 of them (14.6%) in families

Families are short vertical streaks: several events in one place within hours to days. The tallest follows the largest earthquake of the year, the Mw 5.0 of 14 December 2016 that Session 3 studied. Everything else, the large majority, is background.

B.10 · Aftershock sequences or swarms?

Part A (A.10) told two kinds of family apart by the shape of their tree. The depth of an event is the number of links back to the family's first event; a leaf is an event with no children; the average leaf depth ⟨d⟩ is the depth averaged over the leaves. A burst-like family, such as an aftershock sequence, has a star-shaped tree — the mainshock is the parent of most events, so ⟨d⟩ is near 1. A swarm-like family is a chain — each event triggers the next, so ⟨d⟩ is large. Zaliapin & Ben-Zion (2013b) found the two kinds as two populations among families of at least ten events, with ⟨d⟩ above about 5 for the swarm-like ones and no sharp boundary, and warned that small families are hard to classify.

families with at least 5 events: 122; with at least 10: 21
average leaf depth of the families with at least 10 events: 1.0 to 3.7
the five largest families: 33 events (largest M4.7, <d> = 1.0), 30 events (largest M4.4, <d> = 1.0), 25 events (largest M4.6, <d> = 1.1), 15 events (largest M4.1, <d> = 1.2), 14 events (largest M4.0, <d> = 1.4)
the deepest of them reaches 7 generations
the 2016 M5.0: its family has 10 events in this catalogue

Every Geysers family large enough to classify is burst-like: ⟨d⟩ stays below 5. The largest are the aftershocks of M 4 earthquakes, nearly every event a direct child of the mainshock (middle). Even the family of the 2016 Mw 5.0 holds only a few events here, because most of its aftershocks were too small to have a mechanism. The deepest among them (right) is a chain of seven generations, but its average leaf depth is still below 5. The focal-mechanism catalogue lacks the small events that make up most swarms, so few families here are large enough to classify; what the figure shows is that the families we can classify are burst-like, not that swarms are absent from the field.

What clustering found at the Geysers

• In space, the earthquakes concentrate in compact patches about a kilometre across. The number of clusters depends on the scale we choose, but about a fifth of the events sit in patches found at every scale; those patches survive the location errors and, at large scales, are several times as numerous as in a catalogue blurred by half a kilometre.
• In mechanism, the focal mechanisms form one continuum with a preferred style. The dozens of clusters in the reported values come from FPFIT's grid; a reproducible two-way split is a cut through the continuum, not two families.
• Space and mechanism together: more than half of the patches, found from position alone, have a mechanism of their own, more uniform than their neighbours' — normal faulting along the northern edge of the field, strike-slip along its axis.
• In time, fewer than one event in ten is triggered by an earlier one; the families are short-lived, and those large enough to classify are burst-like.

How this compares with published work. Studies of the Geysers usually analyse one such cluster at a time. Kwiatek et al. (2015) followed a single cluster around two injection wells in the north-western field for seven years and found its shape and extent changing with the injection rate; this suggests what the patches here may be — the seismicity around individual injection wells — but the catalogue alone cannot show it, and testing it needs the well data. Martínez-Garzón et al. (2014), studying one such cluster, found its mechanisms mixing normal and strike-slip faulting, with strike-slip mechanisms becoming more common at peak injection: the same two styles that dominate the patches here. Palo et al. (2024) applied unsupervised learning to the space–time density of the whole field and found three concentric domains with different earthquake rates and magnitude distributions, a coarser partition than the kilometre-scale patches here. We have not found published values for the share of triggered events at the Geysers to compare with B.9.

Takeaways

Machine learning

• A clustering result is only as meaningful as the choice of feature, scale and null behind it, and on real data the diagnostics that chose them on a toy (BIC, the k-distance knee) can fail.
• A cluster count means nothing on its own: compare it with catalogues that keep the large-scale shape and remove the structure in question, at several scales.
• Before clustering a feature, find out how it was produced, because a grid in the method that produced it makes clusters of its own.
• A clustering that two random halves of the data agree on can still be a cut through one continuum.
• A feature held back from the clustering is the fairest test that clusters mean something, when it is compared with neighbours that share everything except membership.

Seismology

• Geysers seismicity concentrates in compact patches about a kilometre across that are real, not artefacts of location error or chance.
• The Geysers focal mechanisms form one continuum from normal to oblique and strike-slip faulting, and many patches have a dominant style of their own.
• Fewer than one Geysers earthquake in ten is triggered by an earlier one in the nearest-neighbour sense, and the families large enough to classify are aftershock-like, not swarm-like.
• A b-value from a focal-mechanism catalogue is biased low, because the smallest earthquakes are the ones without enough first motions for a mechanism.

References

• Reasenberg, P., & Oppenheimer, D. (1985). FPFIT, FPPLOT and FPPAGE: Fortran computer programs for calculating and displaying earthquake fault-plane solutions. USGS Open-File Report 85-739. https://doi.org/10.3133/ofr85739
• Kwiatek, G., Martínez-Garzón, P., Dresen, G., Bohnhoff, M., Sone, H., & Hartline, C. (2015). Effects of long-term fluid injection on induced seismicity parameters and maximum magnitude in northwestern part of The Geysers geothermal field. Journal of Geophysical Research: Solid Earth, 120, 7085–7101. https://doi.org/10.1002/2015JB012362
• Martínez-Garzón, P., Kwiatek, G., Sone, H., Bohnhoff, M., Dresen, G., & Hartline, C. (2014). Spatiotemporal changes, faulting regimes, and source parameters of induced seismicity: A case study from The Geysers geothermal field. Journal of Geophysical Research: Solid Earth, 119, 8378–8396. https://doi.org/10.1002/2014JB011385
• Palo, M., Ogliari, E., & Sakwa, M. (2024). Spatial pattern of the seismicity induced by geothermal operations at the Geysers (California) inferred by unsupervised machine learning. IEEE Transactions on Geoscience and Remote Sensing, 62, 1–13. https://doi.org/10.1109/TGRS.2024.3361169
• Zaliapin, I., & Ben-Zion, Y. (2013a). Earthquake clusters in southern California I: Identification and stability. Journal of Geophysical Research: Solid Earth, 118, 2847–2864. https://doi.org/10.1002/jgrb.50179
• Zaliapin, I., & Ben-Zion, Y. (2013b). Earthquake clusters in southern California II: Classification and relation to physical properties of the crust. Journal of Geophysical Research: Solid Earth, 118, 2865–2877. https://doi.org/10.1002/jgrb.50178

The clustering methods and their references are in Part A.