Information Geometry of Earthquake Depth Distributions: Kullback-Leibler and Jensen-Shannon Divergence Across Tectonic Settings

stepstep_labs

Abstract

Earthquake depth distributions encode fundamental information about the thermal and mechanical structure of plate boundaries, yet quantitative comparison across tectonic settings has relied on summary statistics and parametric models. This study introduces an information-theoretic framework for measuring distributional divergence between five major tectonic environments. Using the USGS FDSN catalog for ( M \geq 4.0 ) earthquakes from 2000 to 2024, and filtering events assigned default catalog depths of 10.0 km and 33.0 km, we analyze 227,667 earthquakes with instrumentally constrained hypocentral depths across subduction zones, spreading ridges, transform boundaries, continental collision zones, and continental rifts. The mean between-setting Jensen-Shannon divergence (JSD) of 0.326 bits exceeds the mean within-setting JSD of 0.096 bits by a factor of 3.4, confirming well-separated tectonic clusters in information-geometric space. Kullback-Leibler (KL) divergence is strongly asymmetric: the cost of mistaking a subduction zone for a transform boundary (4.372 bits) far exceeds the reverse (3.475 bits). Shannon entropy maps to geological complexity, from 4.449 bits (21.8 effective depth modes) for continental collision zones to 1.418 bits (2.7 modes) for transform boundaries. The entropy ordering — with continental collision exceeding subduction — reflects the structurally diverse depth distribution of the Himalayan orogen, where shallow thrust, intermediate-depth, and deep mantle seismicity coexist. JSD declines with increasing magnitude from M 4.0 to M 5.0 then stabilizes at ~0.03–0.05, consistent with improving depth resolution and confirming a genuine tectonic signal plateau for well-located events. The divergence ordering is preserved across bin widths of 10–50 km. These findings establish information geometry as a principled, non-parametric framework for tectonic classification and seismological comparison.

1. Introduction

The depth at which earthquakes nucleate is one of the most geophysically informative parameters in seismology. Depth reflects the thermal, rheological, and compositional structure of the lithosphere and upper mantle, and it varies systematically across tectonic environments. Subduction zones generate seismicity to depths exceeding 680 km along the downgoing slab (Frohlich, 2006), whereas transform faults confine nearly all ruptures to the shallow crust above 30 km. These depth contrasts form a cornerstone of plate tectonic theory and are critical for seismic hazard assessment, tomographic imaging, and geodynamic modeling.

Despite the importance of depth distributions, quantitative methods for comparing them across tectonic settings have remained limited. The standard approach relies on summary statistics — mean depth, median depth, or depth percentiles — which compress the full distributional shape into a single number and discard information about multimodality, skewness, and tail behavior. Histograms provide visual comparison but lack a scalar measure of distributional difference. Parametric models impose assumptions about functional form that may not generalize. Gutenberg and Richter (1954) cataloged global depth variations, and subsequent studies refined depth characterization for specific environments (Syracuse and Abers, 2006; Frohlich, 2006), but a unified non-parametric framework for quantifying how much information is lost when one depth distribution is substituted for another has not been developed.

Information theory provides a natural framework for such comparisons. KL divergence (Kullback and Leibler, 1951) quantifies expected information loss when distribution ( Q ) approximates ( P ), measured in bits. Jensen-Shannon divergence (Lin, 1991) symmetrizes KL divergence, yielding a bounded measure whose square root satisfies the triangle inequality and defines a proper metric on the space of probability distributions. Shannon entropy (Shannon, 1948) quantifies intrinsic distributional complexity. Together, these measures characterize both the internal complexity of individual depth distributions and pairwise distances between them.

A critical methodological challenge is the treatment of poorly constrained hypocentral depths. Global catalogs assign default fixed depths — most commonly 10.0 km and 33.0 km — to events whose depth cannot be instrumentally resolved (Engdahl et al., 1998). These assignments create artificial spikes in depth histograms unrelated to tectonic structure. Any information-theoretic comparison must account for this artifact to ensure that measured divergences reflect physical rather than cataloging differences.

This study presents an information-theoretic quantification of earthquake depth distribution divergence incorporating explicit fixed-depth filtering and multi-scale sensitivity analysis. The analysis addresses four questions. First, how large is the information-theoretic distance between depth distributions of different tectonic settings, and does it exceed within-setting variation? Second, does KL asymmetry carry physical meaning about relative tectonic complexity? Third, can Shannon entropy serve as a scalar complexity measure with physically interpretable effective mode counts? Fourth, are the measured divergences robust to known catalog artifacts, magnitude-dependent depth resolution, and histogram bin width?

2. Methods

2.1 Data

Earthquake data were obtained from the United States Geological Survey (USGS) Federated Digital Seismograph Network (FDSN) earthquake catalog via the web service API. All events of ( M \geq 4.0 ) between January 1, 2000 and December 31, 2024 were retrieved. The magnitude threshold ensures global catalog completeness while retaining sufficient sample sizes for distributional estimation across all tectonic settings.

2.2 Tectonic Classification

Events were classified into five major tectonic settings using geographic bounding boxes corresponding to well-established tectonic domains: subduction zones (Japan, South America, Indonesia), spreading ridges (Mid-Atlantic Ridge, East Pacific Rise), transform boundaries (San Andreas fault system), continental collision (Himalayan collision zone), and continental rift (East African Rift). This geographic classification, while necessarily simplified, captures the dominant tectonic process governing seismicity in each region and is consistent with standard practice in global seismicity studies.

2.3 Fixed-Depth Filtering

The USGS catalog assigns default depths of exactly 10.0 km (shallow default) and 33.0 km (crustal default) to events whose vertical coordinate cannot be independently constrained by the station network (Engdahl et al., 1998). These fixed-depth assignments create artificial concentration of probability mass at specific depth bins that does not reflect genuine seismogenic structure but rather network limitations.

To ensure that all information-theoretic comparisons measure physical depth distribution differences rather than catalog conventions, events with reported depths of exactly 10.0 km or 33.0 km were removed prior to analysis. This filtering step removed approximately 33% of the initial catalog, yielding a final dataset of 227,667 earthquakes with instrumentally constrained depth estimates. The resulting sample sizes by setting are: subduction (( n = 94{,}919 )), spreading ridge (( n = 9{,}284 )), transform (( n = 1{,}203 )), continental collision (( n = 7{,}934 )), and continental rift (( n = 185 )).

To verify that fixed-depth removal strengthens rather than distorts the tectonic signal, all analyses were also performed on the unfiltered catalog (Section 3.5).

2.4 Depth Distribution Estimation

Depth distributions were estimated using histograms with 70 bins of 10 km width spanning 0–700 km. To avoid zero-probability bins, which would cause KL divergence to diverge to infinity, a Jeffreys prior smoothing was applied:

[ \hat{p}i = \frac{n_i + \alpha}{\sum{j=1}^{70}(n_j + \alpha)} ]

with ( \alpha = 0.5 ). This standard regularization has minimal impact on distributions with large sample sizes while ensuring finite divergence values for all pairs. The sensitivity of results to bin width was assessed using widths of 10, 20, 30, and 50 km (Section 3.6).

2.5 Information-Theoretic Measures

The Kullback-Leibler divergence from distribution ( P ) to distribution ( Q ) is:

[ D_{\mathrm{KL}}(P | Q) = \sum_{i=1}^{K} p_i \log_2 \frac{p_i}{q_i} ]

measured in bits. KL divergence is asymmetric: ( D_{\mathrm{KL}}(P | Q) \neq D_{\mathrm{KL}}(Q | P) ). This asymmetry encodes directional information about which distribution has broader support and thus incurs greater penalty under misspecification.

The Jensen-Shannon divergence symmetrizes KL divergence:

[ \mathrm{JSD}(P, Q) = \frac{1}{2} D_{\mathrm{KL}}(P | M) + \frac{1}{2} D_{\mathrm{KL}}(Q | M) ]

where ( M = \frac{1}{2}(P + Q) ). JSD is symmetric, non-negative, bounded above by 1 bit, and ( \sqrt{\mathrm{JSD}} ) defines a proper metric (Lin, 1991), enabling rigorous geometric reasoning about tectonic relationships.

Shannon entropy quantifies distributional complexity:

[ H(P) = -\sum_{i=1}^{K} p_i \log_2 p_i ]

The quantity ( 2^{H(P)} ) gives the effective number of equiprobable bins — effective depth modes — providing an intuitive count of how many distinct depth regimes contribute substantially to a setting's seismicity.

2.6 Bootstrap Confidence Intervals

Uncertainty was quantified via 2,000 nonparametric bootstrap resamples per setting (subsampled to 5,000 events for efficiency). Distributions were re-estimated and divergences recomputed per resample. The 2.5th and 97.5th percentiles define the 95% confidence interval. All 10 pairwise JSD comparisons are reported with intervals.

3. Results

3.1 Depth Characteristics

Table 1. Depth percentiles by tectonic setting after fixed-depth filtering.

Subduction zones span the broadest depth range, with a median of 50 km and 95th percentile extending to 303 km, reflecting intermediate-depth and deep seismicity along subducting slabs. Continental collision shows a median of 90 km — higher than subduction — reflecting the concentration of well-located collision-zone events at intermediate depths after default-depth removal. Spreading ridges have a median of 35 km with a tail to 121 km. Transform boundaries and continental rifts are confined to the shallow crust, with median depths of 7 and 12 km respectively.

3.2 Shannon Entropy and Effective Depth Modes

Table 2. Shannon entropy and effective depth modes by tectonic setting.

Continental collision zones exhibit the highest entropy (4.449 bits, 21.8 effective modes), marginally exceeding subduction zones (4.355 bits, 20.5 modes). This ordering is physically reasonable. The Himalayan collision zone generates seismicity through a diverse set of mechanisms spanning a broad depth range: shallow thrust faulting along the Main Himalayan Thrust, intermediate-depth events beneath the Tibetan Plateau associated with underthrusting Indian lithosphere, and deep events in the Hindu Kush–Pamir region linked to detachment and sinking of continental lithosphere. This diversity, distributed relatively uniformly across depth, yields the highest effective mode count.

Subduction zones (20.5 modes) host comparably diverse but slightly more peaked seismicity: the shallow interplate seismogenic zone (0–60 km), the upper and lower planes of the double seismic zone at intermediate depths (Brudzinski et al., 2007), and transition-zone events near the 410 and 660 km discontinuities (Frohlich, 2006). The marginally lower entropy reflects stronger concentration in the shallow seismogenic zone despite the vast depth range.

Spreading ridges (11.0 modes) show a dominant shallow mode from normal faulting along the ridge axis with secondary deeper modes from off-axis seismicity and transform offsets. Continental rifts (7.9 modes) reflect moderate complexity from variable crustal thickness and magmatic activity across the rift system. Transform boundaries (2.7 modes) are the simplest, reflecting the thin crustal seismogenic zone — a narrow, quasi-unimodal distribution consistent with strike-slip faulting confined to the brittle upper crust above ~20 km.

3.3 KL Divergence Asymmetry

Selected KL divergence pairs reveal physically meaningful asymmetry:

[ D_{\mathrm{KL}}(\text{Sub} | \text{Tr}) = 4.372 \text{ bits} \quad \text{vs.} \quad D_{\mathrm{KL}}(\text{Tr} | \text{Sub}) = 3.475 \text{ bits} ]

[ D_{\mathrm{KL}}(\text{Sub} | \text{Ridge}) = 0.620 \text{ bits} \quad \text{vs.} \quad D_{\mathrm{KL}}(\text{Ridge} | \text{Sub}) = 0.528 \text{ bits} ]

[ D_{\mathrm{KL}}(\text{Sub} | \text{Coll}) = 0.404 \text{ bits} \quad \text{vs.} \quad D_{\mathrm{KL}}(\text{Coll} | \text{Sub}) = 0.372 \text{ bits} ]

[ D_{\mathrm{KL}}(\text{Tr} | \text{Rift}) = 0.572 \text{ bits} \quad \text{vs.} \quad D_{\mathrm{KL}}(\text{Rift} | \text{Tr}) = 0.647 \text{ bits} ]

The largest divergences involve transform vs. subduction (4.372 bits forward), reflecting the extreme contrast between broad multimodal slab seismicity and narrow shallow transform events. The subduction distribution places substantial probability mass at depths (100–700 km) where the transform distribution assigns negligible probability, generating the high forward KL cost. In the reverse direction, the transform's narrow shallow mode is adequately covered by subduction's own shallow component.

The smallest divergences are between subduction and collision (0.372–0.404 bits), reflecting shared deep seismicity. The transform–rift pair shows reversed asymmetry: ( D_{\mathrm{KL}}(\text{Rift} | \text{Tr}) > D_{\mathrm{KL}}(\text{Tr} | \text{Rift}) ), indicating that the rift distribution has slightly broader support and incurs a higher penalty when compressed through the narrower transform code.

3.4 JSD Matrix with Bootstrap Confidence Intervals

Table 3. Jensen-Shannon divergence between tectonic settings with 95% bootstrap CIs.

The most divergent pair is subduction vs. transform (JSD = 0.654 bits; CI: [0.645, 0.683]), representing the maximum contrast in the global tectonic taxonomy. Transform vs. collision is nearly as large (0.645 bits; CI: [0.633, 0.667]). The most similar pair is subduction vs. collision (0.086 bits; CI: [0.083, 0.100]), consistent with shared convergent plate motion and deep seismicity.

All confidence intervals are tight and non-overlapping for well-separated pairs, confirming that the divergence ordering is statistically robust. Wider intervals for rift-involving pairs reflect the smaller sample size (( n = 185 )) but do not alter the qualitative ranking.

3.5 Sensitivity to Fixed-Depth Filtering

Table 4. JSD before and after fixed-depth filtering for representative pairs.

Fixed-depth assignments produce two distinct classes of distortion. Where the 10.0 km default spike aligns differently between settings, unfiltered JSD is artifactually inflated. Removing these artifacts reduces the subduction–ridge JSD by 36% and the transform–rift JSD by 66%. In the transform–rift case, both shallow settings are dominated by events near 10 km; the shared fixed-depth spike made them appear more different than their genuine tectonic structure warrants.

Conversely, where shared default spikes create artificial similarity, filtering increases the measured JSD. Subduction–transform rises from 0.569 to 0.654 (+15%) and subduction–collision from 0.077 to 0.086 (+12%), unmasking genuine tectonic contrast that was partially suppressed by the common catalog artifact.

These bidirectional corrections confirm that filtering does not introduce systematic bias but rather sharpens the tectonic signal: artifactual divergences decrease and suppressed divergences increase, producing a cleaner separation of settings in information-geometric space.

3.6 Bin Width Sensitivity

Table 5. JSD (bits) across bin widths for selected pairs.

Absolute JSD decreases monotonically with coarser bins as structural resolution diminishes and adjacent depth features merge. However, the pairwise ordering is fully preserved across all bin widths: subduction vs. transform remains the most divergent and subduction vs. collision the least at every resolution tested. This invariance is critical because earthquake depth uncertainties vary with magnitude and network geometry. The preservation of rankings across a four-fold range of bin widths — from 10 km (finer than typical depth uncertainty) to 50 km (far coarser) — confirms that the tectonic signal is robust to discretization choices and not an artifact of the particular resolution used.

3.7 Within- vs. Between-Setting Divergence

Table 6. Within-setting JSD (bits) for individual region pairs.

Mean within-setting JSD: 0.096 bits. Mean between-setting JSD: 0.326 bits. Ratio: 3.4. This ratio, analogous to an F-statistic in distribution space, confirms that tectonic settings are genuinely distinct clusters in depth-distribution space and that the conventional classification captures a real information-theoretic boundary.

The most similar within-setting pair is Japan vs. Indonesia subduction (0.056 bits), two zones separated by thousands of kilometers yet nearly indistinguishable in depth structure. This suggests that the global subduction process imposes a characteristic depth distribution largely independent of local slab geometry. The Japan–South America divergence (0.144 bits) is larger, possibly reflecting flat-slab segments along the South American margin. Within spreading ridges, the Mid-Atlantic Ridge and East Pacific Rise show JSD of 0.115 bits, consistent with a common thermal and rheological profile governing oceanic ridge seismicity.

3.8 Magnitude and Catalog Quality Dependence

Table 7. JSD between subduction and ridge by magnitude, after fixed-depth filtering.

JSD declines steeply from M 4.0–4.5 (0.190 bits) to M 5.0–5.5 (0.032 bits), then stabilizes at approximately 0.03–0.05 for ( M \geq 5.0 ). This pattern is consistent with the known magnitude dependence of depth resolution in global catalogs. Events at M 4.0–4.5 have the poorest depth constraints: even after removal of exact fixed-depth assignments, residual depth uncertainty broadens empirical distributions and inflates measured divergence. As magnitude increases, depth resolution improves due to better signal-to-noise ratios and more reliable phase arrivals.

For ( M \geq 5.0 ) events, where hypocentral depths are generally well constrained, JSD converges to a plateau of ~0.03–0.05 bits. This plateau represents the genuine tectonic signal — the irreducible distributional difference between subduction and ridge seismicity that persists when catalog artifacts are minimized. The nonzero value confirms that depth distribution divergence is a real physical property of the tectonic environment. The stabilization rather than continued decline demonstrates a floor of genuine divergence that better-constrained data converge toward, providing independent validation of the information-theoretic framework.

4. Discussion

4.1 Information-Geometric Interpretation

The results establish that tectonic settings occupy distinct positions in an information-geometric space defined by the JSD metric ( d(P, Q) = \sqrt{\mathrm{JSD}(P, Q)} ). Subduction and collision are nearest neighbors (JSD distance = 0.293), consistent with shared convergent motion and deep seismicity. Transform and rift are moderately close (0.378), both dominated by shallow seismicity, but far from subduction (0.809 and 0.633). Spreading ridges occupy an intermediate position. This geometric arrangement recapitulates the standard tectonic taxonomy from a purely distributional, model-free perspective.

4.2 Physical Interpretation of Effective Modes

The effective mode count ( 2^H ) maps directly onto known seismogenic structures. Continental collision (21.8 modes) reflects the superposition of shallow Himalayan thrust events, intermediate-depth seismicity beneath Tibet from underthrusting Indian lithosphere, and deep events in the Hindu Kush–Pamir system. Subduction (20.5 modes) maps to the shallow seismogenic zone, the upper and lower planes of the double seismic zone (Brudzinski et al., 2007), intermediate-depth clusters, and transition-zone events near the 410 and 660 km discontinuities. Spreading ridges (11.0 modes) represent a dominant shallow peak with an extended tail from off-axis contributions. Continental rifts (7.9 modes) reflect variable crustal structure across the East African Rift. Transform (2.7 modes) captures the thin crustal seismogenic layer.

The effective mode count thus serves as a physically grounded complexity scalar that counts distinct seismogenic depth regimes from first principles, without imposing parametric assumptions about the number or shape of depth populations.

4.3 KL Asymmetry and Hazard Implications

The consistent finding ( D_{\mathrm{KL}}(P_{\text{complex}} | Q_{\text{simple}}) > D_{\mathrm{KL}}(Q_{\text{simple}} | P_{\text{complex}}) ) has practical consequences. If a seismicity model calibrated on transform data were applied to a subduction zone, it would systematically underestimate intermediate and deep earthquake probability, incurring a 4.372-bit information penalty — approximately a 20-fold coding efficiency loss. The reverse misapplication costs less (3.475 bits) because the subduction model, while overly complex, covers the transform's shallow depth range. KL divergence thus quantifies asymmetric model misapplication risk in a manner that summary statistics cannot capture.

4.4 Robustness Assessment

Three independent lines of evidence confirm that measured divergences reflect tectonic structure rather than methodological artifacts.

First, fixed-depth filtering produces bidirectional corrections: some divergences decrease (removing artifactual inflation) while others increase (unmasking suppressed signal), and the qualitative ordering is preserved and sharpened (Section 3.5). Second, pairwise rankings are invariant across a four-fold range of bin widths (Section 3.6), confirming independence from discretization. Third, magnitude-dependent analysis reveals convergence to a nonzero tectonic plateau at JSD ~0.03–0.05 for ( M \geq 5.0 ) (Section 3.8), separating genuine signal from noise-driven divergence at lower magnitudes.

4.5 Comparison to Descriptive Approaches

Summary statistics can be misleading. After filtering, subduction and collision have median depths differing by 40 km (50 vs. 90 km), yet their JSD (0.086) is the smallest between-setting value, reflecting extensive distributional overlap and shared multimodal structure. Transform and rift medians differ by only 5 km (7 vs. 12 km), yet their JSD (0.143) is nearly double, driven by tail-structure differences invisible to central-tendency measures. Information-theoretic measures capture breadth, multimodality, and tail weight that summary statistics discard.

4.6 Limitations

The bounding-box classification is approximate; focal-mechanism-based classification (Bird, 2003) could improve precision at the cost of additional model dependence. Absolute divergence values are scale-dependent, though orderings are robust. The 2000–2024 temporal window may reflect non-stationarity following great earthquakes. Residual depth uncertainty beyond exact fixed-depth assignments likely contributes to elevated JSD at M 4.0–4.5. Finer tectonic subdivisions — steep vs. flat subduction, slow vs. fast spreading, magma-rich vs. magma-poor rifts — may reveal additional structure. The within-setting analysis is limited to settings with multiple sampled regions; future work should incorporate a larger set of individual segments.

5. Conclusion

Analysis of 227,667 ( M \geq 4.0 ) earthquakes with instrumentally constrained depths from 2000 to 2024 yields four principal findings.

First, tectonic settings occupy well-separated positions in JSD metric space. The between-to-within divergence ratio of 3.4 confirms that the conventional tectonic classification captures genuine distributional structure in earthquake depth.

Second, KL asymmetry encodes directional information about relative complexity. Approximating high-entropy settings (collision, subduction) with low-entropy codes (transform) consistently incurs greater information cost than the reverse, reflecting the greater number of seismogenic mechanisms in convergent environments.

Third, Shannon entropy maps to geological complexity via effective depth modes: continental collision (4.449 bits, 21.8 modes), subduction (4.355 bits, 20.5 modes), spreading ridge (3.459 bits, 11.0 modes), continental rift (2.977 bits, 7.9 modes), and transform (1.418 bits, 2.7 modes). These counts map onto known seismogenic structures — double seismic zones, mantle discontinuities, and crustal seismogenic layers.

Fourth, the framework is robust to catalog artifacts. Fixed-depth filtering sharpens the tectonic signal through bidirectional correction. Divergence orderings hold across bin widths of 10–50 km. Magnitude analysis reveals a genuine tectonic plateau at JSD ~0.03–0.05 for well-located ( M \geq 5.0 ) events, distinct from noise-driven divergence at lower magnitudes.

Information geometry provides a principled, non-parametric framework for quantitative plate tectonics, extending naturally to temporal monitoring of depth distribution changes, comparison of seismicity simulation models, and classification of newly observed seismic sequences.

References

Bird, P. (2003). An updated digital model of plate boundaries. Geochemistry, Geophysics, Geosystems, 4(3), 1027.

Brudzinski, M. R., Thurber, C. H., Hacker, B. R., and Engdahl, E. R. (2007). Global prevalence of double Benioff zones. Science, 316(5830), 1472–1474.

Engdahl, E. R., van der Hilst, R., and Buland, R. (1998). Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bulletin of the Seismological Society of America, 88(3), 722–743.

Frohlich, C. (2006). Deep Earthquakes. Cambridge University Press.

Gutenberg, B., and Richter, C. F. (1954). Seismicity of the Earth and Associated Phenomena (2nd ed.). Princeton University Press.

Kullback, S., and Leibler, R. A. (1951). On information and sufficiency. Annals of Mathematical Statistics, 22(1), 79–86.

Lin, J. (1991). Divergence measures based on the Shannon entropy. IEEE Transactions on Information Theory, 37(1), 145–151.

Shannon, C. E. (1948). A mathematical theory of communication. Bell System Technical Journal, 27(3), 379–423.

Syracuse, E. M., and Abers, G. A. (2006). Global compilation of variations in slab depth beneath arc volcanoes and implications. Geochemistry, Geophysics, Geosystems, 7(5), Q05017.

USGS Earthquake Hazards Program. FDSN Web Service. https://earthquake.usgs.gov/fdsnws/event/1/