Earthquakes unveil the global-scale fractality of the lithosphere

The relationship between the magnitude of earthquakes and their spatial and temporal distribution has been observed to exhibit a scale invariance hypothesised to originate from self-organized critical regimes. However, the fractality of earthquake distributions has been mostly established in circumscribed areas, despite the fact that the self-organized criticality of the lithosphere should only emerge at global or continental level. Here, we analyze seismic observations occurring over the whole Earth between 2004 – 2020 to investigate the fractal correlation dimension of earthquakes distribution. We ﬁ nd that the distribution of earth-quakes is fractal on a global scale, as well as approximately magnitude-independent and stationary over decadal time scales. Our results set a primary constraint on the spatial scaling properties of lithosphere dynamics. We suggest that macroscopic models should ful ﬁ l this constraint to correctly replicate the features of seismicity, and potentially improve seismic hazard assessment.

T he Earth's lithosphere is a complex dynamical system whose evolution covers different temporal and spatial scales, ranging from the million-year-long reorganization processes of tectonic plates down to the almost continuous occurrence of earthquakes 1 .The statistical features of earthquake occurrence provide experimental signatures of the lithosphere's complex, nonlinear dynamics [2][3][4] .The best-known of these features is the Gutenberg-Richter (G-R) law 2,5 , which corresponds to a scale invariance, i.e., a power-law scaling, of the distribution of the released energy.The observation of this scaling law prompted the idea that the lithosphere is in a regime of self-organized criticality (SOC) [6][7][8] .SOC is generally considered to be a spontaneous, collective organization of an externally-driven system towards a stationary state that shows power-law-like, i.e., selfsimilar, distributions of observable quantities 9 .A scale invariance also exists in the time domain: the Omori law 10,11 describes the power-law decay, and thus the persistent autocorrelation, of the rate of aftershocks that ensue from a main event.While basic models of SOC such as the sandpile model 6 do not predict aftershocks 3 , extensions of these models were shown to be able to reconcile the SOC paradigm with the experimentally-observed Omori law 12,13 .Remarkably, by means of natural time analysis 14 to extract information from time series, it was recently proven that SOC does not imply unpredictability of earthquakes 15 .
A third feature that is commonly ascribed to the complexity of the Earth's lithosphere, and thus related to SOC as well, is fractality, which corresponds to a scale-invariance in the spatial domain 2,3 .Fractality, which is typically assessed by estimating a noninteger dimension 16 referred to as correlation dimension, was applied to describe geological data 4,17 and, specifically, the spatial distribution of earthquakes 18 .It is worth noting that other fractal dimensions exist, most notably the Hausdorff dimension and the information dimension.However, the correlation dimension is the most informative one in the context of seismicity because it describes the probability that, given an earthquake, another earthquake occurs within a given distance 16,18 .
Understanding the origin of the lithosphere complexity is ultimately linked to the possibility of earthquake forecasting 3 .Describing the geometry of a system basically means to determine its spatial extent.This simple statement is also true for fuzzy systems like fractal structures: in a system having fractal dimension close to 1, events occur along a not-necessarily-straight line.Similarly, once an event belonging to a system with dimension close to 2 is observed, the probability of observing a second one essentially depends only on the distance from the first occurrence.
On the other hand, a thorough understanding of the fractality of earthquake occurrence, and thus of the lithosphere, is still lacking.Existing mechanical and statistical models, like the Burridge-Knopoff spring-block model 19,20 and its cellular automaton transposition 21,22 , provide a partial description of the scale invariance features introduced above 23 .More recently, sophisticated models including friction heterogeneities by means of both a brittle and a ductile region within a fault allowed to capture more in detail the complex features of instrumental earthquakes 24 .Chaos was proposed as a key to describe complexity in the lithosphere, both in numerical models [25][26][27] and experimental studies [28][29][30] .An attempt of unifying G-R law, Omori law and fractality within a single model makes up a promising approach, but is still incapable to predict the scaling law parameters 31 .
Whereas the experimental determination of G-R and Omori law follow from the application of standard statistical tools, the fractal geometry of earthquake distributions is less straightforward to evaluate 18 .In the last two decades, this issue has been tackled with increased attention and by means of different approaches.Some works [32][33][34][35][36][37][38][39] deal with events occurring over time spans of many years and within bounded geographical regions of size ranging from ~50 km up to ~3000 km.Another set of works [40][41][42][43][44] concern sequences of small-magnitude aftershocks following a main event: these studies analyze regions of size up to ~500 km and cover shorter time spans, with some works considering time intervals as short as a few months 41,43 .Noteworthily, fractal dimension was mostly evaluated out of epicenter locations.While the choice of analyzing epicenters can be justified by high uncertainty affecting depth estimations 45 , considering epicenters or hypocenters can lead to different results 18 .In terms of the assessed fractal dimension, most works provide noninteger values between 1 and 2, thus suggesting a fractal structure that is halfway between a linear and a planar one.Many studies also found sub-linear dimensions (<1) that are associated to clustered distributions, and even dimensions greater than 2 in cases in which hypocenters are considered 36,40 .A fractal dimension close to 2 was found by analyzing H/V spectral ratio of seismic noise 46 .Occasionally, unphysical dimensions greater than 2 are reported in works dealing with epicenters [47][48][49] .
The analysis of the existing scientific literature shows that the evaluation of fractal dimension was mostly carried out on a local scale, typically not exceeding 1000 km, and with poorly populated samples of just a few hundred events.Despite its practical interest, for example, to assist the creation of risk assessment maps 50,51 , this approach does not provide a characterization of the lithosphere in its entirety.On the other hand, the regime of SOC might hold only on a global or continental scale 52 , so that the understanding of the lithosphere dynamics cannot occur without a whole-Earth analysis of fractal dimension.
A few works examined fractality at the global level.In two of them 53,54 , the fractal dimension was estimated out of the b value of the G-R law.However, the resulting, and bereft of any uncertainty, assessment of 1.78 depends on the model used to link the frequency-magnitude distribution to the spatial distribution 55,56 .In 1980, by analyzing the spatial distribution of "only" 429 events, Kagan and Knopoff 57 provided a fractal dimension of 1.5, though admitting that this value was affected by "large uncertainties".In a later work, Kagan 18 suggested values in the range 1.5-2.2.Similarly, by looking at large-scale regions of order ~5000 km, ref. 58 estimated the fractal dimension to be in the range 0.75-1.59.
In the present work, we estimated the fractal dimension of the spatial distribution of earthquakes occurring on the whole Earth.Due to the extent of the dataset, consisting of ~2 ⋅ 10 5 events, and the statistical analysis implemented, our estimate of ν ¼ 1:30 ± 0:03 is the first one to provide robust, model-free evidence of global-scale fractality.Records were extracted from the earthquake catalog of the United States Geological Survey (data selection is discussed in the Methods section).Because the catalog provides also depth estimates, we considered here both hypocenters and epicenters, and compared the respective results.To estimate fractal dimensions, we computed the correlation dimension 59,60 out of sets of points that correspond to the locations of events.Rather than following the widespread approach of least-squares fitting a power-law on the correlation integral 59 , we carried out the evaluation by relying on Takens estimator 16 , which is more efficient and easier to automatize 61,62 .The estimator was combined with a bootstrapping approach to evaluate the related uncertainty.To strengthen our results, we also compared estimates obtained by means of different distance metrics.While our main goal concerns the estimation of fractal dimension on a global scale, we also performed the same analysis on a local scale (~1000 km) as a mean of comparison with existing estimates.Finally, we investigated possible dependencies of fractal dimension on time and magnitude.

Global fractal dimension.
We considered here four distance metrics (see Methods): besides the spherical (orthodromic) distance (SP), namely the shortest distance measured along the surface of a sphere, and the Euclidean 3-d distance (E3), we also used two locally-defined distances, henceforth referred to as flat 2-d (F2) and flat 3-d (F3) distance, to highlight the essential independence of the results on the choice of the metric.In the case of epicenters, the metrics SP, F2, E3 were considered; for hypocenters, F3, E3 were used instead.
Indeed, a positive bias of less than 0.1 is present between fractal dimensions estimated with the E3 metric and with other metrics.The presence of a bias in estimates of the fractal dimension of an object residing on a spherical surface is theoretically motivated: the inequality ν E3 > ν SP was proven to always hold, and the extent of the bias was quantitatively assessed for the integer dimensions 1 and 2 63 .As discussed in Supplementary Note 2, similar discrepancies are also observed in three synthetic datasets: a stochastic set having integer dimension 2; a chaotic attractor lying on a sphere; a more realistic dataset generated by an epidemictype aftershock sequence model 64 seeded with the same chaotic attractor.In all cases, the metric E3 provides the correct estimate, while SP and F2 yield values that are systematically lower by approximately 0.1.
It is worth remarking that, considering the E3 metric, hypocenters yield a slightly higher fractal dimension than epicenters, due to the hypocenter distribution extending over an additional dimension, i.e., depth.The same discrepancy is also observed by comparing the F2 metric for epicenters with the F3 metric for hypocenters (F3 is indeed an extension of F2 that includes depth).These differences might have a geophysical origin, possibly linked to the depth distribution of events, that could prompt further investigations.Nevertheless, the differences between hypocenters and epicenters are of the same order of the respective uncertainties.This result suggests that the global analysis queries the distribution of earthquakes over distances that are larger than the typical depth of events, namely tens or hundreds of kilometers.To support this interpretation, we performed an auxiliary analysis-discussed in Supplementary Note 3-by increasing the depths of events through the addition of uniformly-distributed random depth displacements.As shown in Supplementary Fig. 10, the larger the displacement, the more the increased-depth hypocenter fractal dimension differs from the epicenter one.The difference overcomes the uncertainties affecting the two assessments when the displacement range is about 10% of the Earth radius, i.e., approximately 600 km, which also occurs to approximately correspond to the maximum depth observed in real earthquakes.
A final consideration that is worth highlighting is the consistency between the uncertainties estimated from the results of Fig. 1 and those assessed in Supplementary Note 2 for the three synthetic datasets mentioned above (see Supplementary Figs.6-9).The extent of these uncertainties, although stemming from independent datasets, are comparable, ranging between 0.01 and 0.02.The similarity also concerns the shape of the histograms.These facts hint at a numerical origin of the experimentally observed distributions of the fractal dimension ν.
Local fractal dimension.The estimation of local fractal dimension is carried out on 10 ∘ × 10 ∘ sectors (latitude × longitude, see Methods), so that the Earth curvature becomes negligible.The different metrics, which provide-besides a 0.1 bias-consistent results on a global scale, are then expected to provide a fortiori similar results also on a local scale.This expectation is confirmed by the analysis reported below.For the sake of simplicity, we therefore discuss here local fractal dimension estimates only for the E3 metric, which also allows direct comparison between epicenters and hypocenters.It should be noted that the evaluation was carried out only on those sectors containing at least 500 events (see Methods).
The local fractal dimension of epicenters is shown in Fig. 2a as a color map on the surface of the Earth, whereas the distribution of local ν estimates is shown, as a histogram, in Fig. 2b.The local fractal dimension ν takes on values that range from ~0.5 to ~2.According to the criterion based on the distance to the closest integer value (see Methods), out of the 480 subsectors for which ν is available, 434 are deemed to provide a noninteger value of fractal dimension at a 95% confidence level.These subsectors contribute to the colored portion of histogram bars in Fig. 2b, while the 46 subsectors compatible with an integer (i.e., not fractal) dimension account for the gray portion of the same bars.In summary, most subsectors (90%) are incompatible with an integer dimension, thus providing a strong evidence that the fractality of the distribution of epicenters is a local property, in addition to being a global one.As mentioned above, the local fractal dimension of epicenters computed with the SP and F2 metrics, respectively, shown in Figs. 3 and 4, provide similar results.
The local fractal dimension of hypocenters is shown in Fig. 2c, d as a color map on the surface of the Earth and as a histogram of local ν estimates.Again, the colored portion of histogram bars are contributed by subsectors that are deemed to provide a noninteger value of fractal dimension at a 95% confidence level: of the 424 subsectors for which ν is available, 344 (81%) provide a noninteger value of fractal dimension.The results obtained with the E3 metric are similar to those obtained with the F3 metric, as shown in Fig. 5.
The results obtained from the analysis of hypocenters are qualitatively similar to those obtained for epicenters.As discussed in Section "Global fractal dimension", this outcome is indeed expected for global properties, as they are weakly affected by the introduction of depth due to the difference of orders between typical depth values and the Earth radius.On the other hand, the local landscape could in principle be significantly modified when depth is taken into account.To address this point, we carried out a comparison between local fractal dimension estimates for epicenters and for hypocenters.Figure 2f plots, for each location, the local fractal dimension value obtained for the hypocenter versus the corresponding epicenter value.In both cases, the metric E3 was used.A color map of the differences of local fractal dimension Δν ¼ ν hyp;E3 À ν epi;E3 is shown in Fig. 2e.In most cases, the difference is positive.This fact is consistent with epicenters being projections of hypocenters onto a twodimensional manifold, an operation that can reduce dimension 18 .Indeed, the set of points colored in red in Fig. 2f, for which 1:5 ≤ ν epi;E3 ≤ 2 and 2 ≤ ν hyp;E3 ≤ 2:5, can be interpreted to correspond to hypocenters distributed on a nonplanar manifold, while the corresponding epicentral locations, which are projected on a locally planar spherical surface, are constrained to distribute with a dimension ν ⩽ 2.
The sporadic observation of sectors having ν epi;E3 > ν hyp;E3 (blue points in Fig. 2f) can be qualitatively explained as follows.One can consider a finite, worm-like set of points, distributed in a more-than-linear fashion, so as to yield an experimentally assessed dimension between 1 and 2. If the set is elongated along a given direction, projecting it on a plane perpendicular to that direction can yield a planar random distribution, whose fractal dimension can be estimated as approaching 2. This example shows that projecting a set of points distributed on a trench The gray portions of the histogram bars count the subsectors for which the estimate of fractal dimension is deemed to be compatible with an integer value (see Methods).In a, c, the same subsectors are also reported in gray.e Map of the difference Δν ¼ ν hyp;E3 À ν epi;E3 of local dimension between hypocenters and epicenters.f Scatter plot of the values of local fractal dimensions for hypocenters versus those for epicenters: each data point corresponds to a subsector in which both assessments are available.The black, dashed lines highlight integer dimensions.In panels e, f the color scale is the same: red hues denote a higher dimension for hypocenters, while blue hues correspond to a higher dimension for epicenters.In a, c, e the Earth surface is mapped via an equirectangular projection.
inclined with respect to the Earth surface on the surface itself can occasionally induce the sporadic behavior observed, namely ν epi;E3 > ν hyp;E3 .It is worth noting that the occurrences of ν epi;E3 > ν hyp;E3 correspond to both fractal dimensions being within 1 and 2, and that most of the related sectors (Fig. 2e) are located in correspondence of trenches (e.g., South America trench, Aleutian trench, Japan trench).
A thorough understanding of these discrepancies unavoidably calls for a geophysical modelization of the underlying earthquake distributions.Time dependence of global fractal dimension.The above analysis was carried out on the period between 2004 and 2020.A point worth investigating is whether and how the fractal dimension changes in time.The global fractal dimension as a function of the year is shown in Fig. 6.In most cases, the yearly values of fractal dimension are compatible, within the respective uncertainties, with the value for the whole period, which can be seen as an evidence of the constancy of the global fractal dimension all along the time span considered (17 years).Three points exhibit significantly lower fractal dimension, namely the years 2005, 2011 and, to a lesser extent, 2010.
These three cases are affected by a significant, low-dimensional contribution.Each year is indeed characterized by a largemagnitude event: the m = 8.6 Indonesia earthquake of 28 March 2005 65 ; the m = 8.8 Chile earthquake of 27 February 2010 66 ; the m = 9.1 Japan earthquake of 11 March 2011 67 .These events were followed by a large number of smaller-magnitude aftershocks, which account for a relevant fraction (>15%) of that year's total number of events.It is, therefore, reasonable to hypothesize that the global fractal dimension was affected by the point-like distribution-as viewed on a global scale-arising from the aftershocks of these large events.In these cases, despite the exclusion of correlated events (see Methods)-which relies on an algorithm for cluster identification but does not collapse a cluster to its mainshock-aftershocks associated to large events are so many that their distances with other, uncorrelated events contribute significantly to the estimation of fractal dimension.
To verify this interpretation, we carried out an additional, trimmed computation of the global fractal dimension on the years 2005, 2010 and 2011.Prior to the analysis, and for each of the three years, we excluded the respective large-magnitude event from the corresponding dataset.Moreover, all the successive events of that year that occured within a 10 ∘ radius around the large event were also excluded (see Eq. (3) in Methods).To prevent analysis biases, we applied the same exclusion approach to years 2004 and 2012, both containing events of comparably large magnitude, i.e., larger than 8.5 (the m = 9.1 Indonesia earthquake on 26 December 2004 68 and the m = 8.6 Indonesia earthquake on 11 April 2012 69 ).
The global fractal dimensions resulting from the trimmed computations for the years 2004, 2005, 2010, 2011, and 2012 are shown in Fig. 6 alongside the originally computed values.The dimensions for the years 2005, 2010, and 2011 indeed turn out to increase and, in the case of 2005 and 2011, become consistent with the values of the years devoid of large events and of the whole period.In the case of 2004 and 2012 the values are essentially unchanged.This outcome can be explained by considering the fractions of events within each year that were excluded in the trimming procedure: these fractions are equal to 5% in 2004, 16% in 2005, 16% in 2010, 31% in 2011 and 6% in 2012.In other words, a reduction of fractal dimension due to aftershocks following a large event is observed only if the corresponding number of events accounts for 15% of the total number of events in the dataset.
Considering the trimmed analysis, the global fractal dimension is essentially stationary, although some residual variability persists.The time dependence of this variability does not appear to follow any obvious law.Attempting to fit it with a model bereft of a theoretical justification-for example, a linear law-can lead to unreliable results.Also, any time dependency might occur over time scales that are significantly different from the ones analyzed in Fig. 6.Addressing these issues would require higher statistics, i.e., a lower magnitude of completeness, and/or a longer temporal coverage as compared to the two decades analyzed in the present work.
One could argue whether the yearly residual variability observed in Fig. 6 might have a statistical origin, for example, a dependence on the number of events included in the estimation.To test this possibility, we analyzed correlations between yearly global fractal dimension and the related number of events.A scatter plot of the two quantities is reported in Fig. 7.The Spearman correlation coefficient r between number of events and global fractal dimension was computed and tested against the null hypothesis of uncorrelated data.The resulting correlation coefficients r and the corresponding p-values are r = − 0.27, p = 0.31 for epicenters and r = − 0.30, p = 0.23 for hypocenters.In both cases, as p-values are well above the 5% significance level, the null hypothesis of uncorrelated data is not rejected.
The residual dimension variability observed here cannot be straightforwardly explained: tackling this issue, e.g., by correlating it with other lithosphere observables, can provide an additional key to understand the inner dynamics of the Earth.

Magnitude dependence of global fractal dimension.
A further point worth investigating is a possible dependence of fractal dimension on magnitude.To carry out this analysis, events were divided into three subsets by considering the magnitude intervals 4; 4:5 ½ Þ; 4:5; 5:5 ½ Þ, and 5:5; 1 ½ Þ (see Methods).The global fractal dimension for each magnitude range is shown in Fig. 8.While a clear dependence on magnitude does not emerge from these results, an interesting observation is that the highest-magnitude group, though affected by higher uncertainties, exhibits a marked discrepancy between epicenters and hypocenters dimension.This fact suggests that depth might play a role only in the case of higher magnitudes (m ≥ 5.5).

Discussion
In this work we analyzed, in terms of its fractal -or correlation -dimension, the distribution of earthquakes having magnitude ≥ 4 and occurring on the whole Earth over a 17-year period (2004-2020).We found that the distribution of earthquakes exhibits a noninteger fractal dimension on a global scale.To our knowledge, this is the first assessment that provides a robust, model-free evidence of global fractality by taking into account events occurring over the whole lithosphere and including a large number of events (>2 ⋅ 10 5 ).In addition, we evaluated the fractal dimension on a local scale, obtaining noninteger values that are in agreement with regional analyses reported in the literature and are distributed about the value assessed on a global scale.
We showed that global fractal dimension does not significantly change whenever the analysis is restricted to a subset of events, either by time segmentation or by considering different magnitude ranges.This last observation is consistent with the scale invariance described by the G-R law: as no characteristic magnitude determining earthquake occurrence can be identified, the geometry of the distribution of events is essentially the same at any magnitude.This result suggests a possible, future development of the work in which the role of catalog incompleteness is investigated.A possibility could be to introduce artificial incompleteness in a synthetic epidemic-type aftershock sequence model 70 and to assess the resulting effect on the fractal dimension.Along the same line, it would be interesting to compare the present results with those obtained from local catalogs with lower magnitude of completeness.
By analyzing the distributions of both epicenters and hypocenters, without a priori restricting to one of the two kinds of data as typically done in the literature, we found that on a global scale the fractal dimensions of epicenters and hypocenters are statistically compatible.This observation implies that the structures probed by the global analysis mostly concern the distribution over lengths larger than tens or hundreds of kilometers.On the contrary, the assessment of local fractal dimension is affected by the inclusion of event depths.As shown in Fig. 2, local fractal dimensions for hypocenters tend to be higher than epicenter values, although this is not a general rule for all geographical regions.Specifically, larger hypocenter dimensions, which in some cases are greater than 2, are mostly detected in the western Pacific area.These findings highlight the fact that analyzing epicenter or hypocenter locations does not necessarily provide equivalent results.On the contrary, it would be interestingthough beyond the scope of the present work-to investigate the origin of such difference in terms of both the dynamics of the lithosphere and its structural properties.
The local correlation dimension being distributed between 0.5 and 2.5 (see Fig. 2), along with its being well-defined on a global scale as 1.30 ± 0.03 (see Fig. 2), is a marker of multifractality, a property shared by most naturally-occurring and synthetic complex systems 71,72 .In fact, if it were a monofractal, any local subset would have the same structure as the whole object and thus the same local correlation dimension.A multifractal analysis goes beyond the scope of the present work.However, as mentioned in the Introduction, it has to be stressed that it is the correlation dimension ν, rather than other dimensional estimates like the information dimension and the Hausdorff dimension, that is relevant for describing the spatial concurrence of earthquakes.
The found value for the fractal dimension of about 1.3 is characteristic of many other phenomena associated to the system Earth.For example, the coastlines of Britain have a fractal dimension of 1.25 73 , a value that is close to the one of 1.3 predicted by stochastic models of islands' relief 74 .Similarly, topography and bathymetry profiles exhibit a fractal dimension of about 1.2 54 .6][77] .
Interestingly, comparable-though indirect-estimates of fractal dimension were independently obtained by means of natural time analysis of seismicity 14 : an estimate of 1.32 ± 0.06 for the Mediterranean region 78 (corresponding, in the present work, to 33 subsectors for which an estimate of ν is available, and whose average fractal dimension is 1.6 ± 0.3); an estimate of 1.398 ± 0.019 globally, considering events of magnitude m > 4.9 79 ; a slightly smaller value, namely 0.95 ± 0.02, for the Pacific Coast  Scatter plot of the global fractal dimension of epicenters (a) and hypocenters (b), obtained by means of the E3 metric, versus the number of events.Each dot and the related error bar correspond to the value of one year among those shown in Fig. 6, and is reported with the same color: gray dots correspond to the years for which the fractal dimension was computed upon exclusion of large events (see main text).
of North and Central America 80 (in this case, the 9 subsectors for which ν is available correspond to an average fractal dimension of 1.1 ± 0.1).
The similarity of all these values can be attributed to the expression of the underlying internal dynamics of the planet in proximity of its surface, which leads to properties that are fractal in space and chaotic in time.Moreover, the fractality of the earthquake distribution can be explained by the distribution of faults being fractal as a result of plate tectonics, where tectonic forces fragment the continental crust into a fractal distribution of interacting crustal blocks over a wide range of spatial scales; the crustal blocks are bounded via faults, so that a fractal distribution of block size can be associated with a fractal distribution of faults 54 .Finally, these mechanisms of plate tectonics shaping the Earth's surface are likely activated by the chaotic, intermittent energy transport mechanisms that result from mantle convection.
The present, unambiguous estimation of the global fractal dimension of earthquake distribution provides a robust constraint for any dynamical model of the lithosphere, where seismic phenomena have to occur according to well-established scaling laws.

Methods
The earthquakes considered here were extracted from the earthquake catalog of the United States Geological Survey (USGS) 81 .The analyses presented in this work were carried out on different sets of earthquakes.Earthquakes are highly localized in space and time with respect to the typical scales of the lithosphere size and evolution 1 , and we can therefore consider an earthquake occurrence to be described by a point process.Consequently, throughout the present work, the terms earthquake and event are used interchangeably.Each event E i , where i ¼ 1 I , is labeled by latitude θ i , longitude λ i , depth ζ i , magnitude m i , and occurrence time t i .This last parameter, which is provided in the catalog as a UTC timestamp having millisecond resolution, is expressed as Unix epoch time: t i corresponds to the number of seconds that have elapsed since January 1st, 1970, 00:00:00 UTC up to the occurrence of the i-th event.
Data selection.The USGS catalog was queried for all earthquakes of magnitude greater than or equal to 2.5 and occurring in the years from 1980 to 2020 inclusively.This raw dataset was then suitably trimmed both in time and by determining the so-called magnitude of completeness M c as follows.
First, to assess the time interval in which data are most complete, the yearly number of events having magnitude m ≥ 4.5 was considered, as shown in Fig. 9. Two distinct time periods can be identified, namely years up to 2003 and years from 2004 on, as highlighted in Fig. 9 by different colors: indeed, the second period exhibits a systematically higher event count than the first one.This fact is possibly due to the improvement of instrumental sensitivity and geographical coverage of the detection networks, a trend that is visible also during the first period.Consequently, the period 2004-2020 was selected as the one that provides the most complete data.
In the second step, an accurate value of the magnitude of completeness M c was evaluated on the subset of events corresponding to this last period.To this purpose, the frequency-magnitude distribution of events, namely the histograms of the number NðMÞ of events having m equal to and greater than M were considered: NðMÞ ¼fE i j m i ≥ MÞg.Thereupon, the so-called Goodness-of-Fit Test 82 was applied by fitting a G-R law of the kind NðMÞ ¼ 10 aÀbM on the frequencymagnitude data (the residual threshold parameter for the Goodness-of-Fit Test was set to 5% as in ref. 83 ).The corresponding magnitude of completeness M c turns out to be equal to M c = 4, while the exponent b of the G-R law was found to be b = 1.02 ± 0.01 and a = 9.70 ± 0.06.The parameter M c = 4 is given without error: M c is used here as a threshold value for further numerical assessments, so that an error assignment is unnecessary.The Goodness-of-Fit approach for the evaluation of M c is less stringent than other available criteria, as we aimed at maximizing the size of the analyzed dataset.However, it is worth remarking that, as it can be inferred from the results of fractal dimension estimation as a function of magnitude (see Fig. 8), changes of M c of order ~0.5 would only marginally affect the fractal dimension estimate.
The main dataset analyzed in the present work, referred to as D 0 , is thus given by the set of all events having magnitude m ≥ 4 and occurring between the beginning of 2004 and the end of 2020: where t o and t f correspond to 2004-01-01, 00:00 UTC and 2020-12-31, 23:59 UTC, respectively.The cardinality D 0 , i.e., the total number of events contained in D 0 , is D 0 = 238384.The epicentral locations of all events in the D 0 dataset are shown in Supplementary Fig. 1.To show that the dataset correctly samples events without biases in depths and magnitude, the conditional probability distribution p(ζ|m) of event depths ζ given their magnitudes m is shown in Supplementary Fig. 2. Throughout the work, additional analyses were carried out on subsets of D 0 .More specifically, a time-dependent analysis of fractal dimension relies of yearly sets of events obtained by partitioning the D 0 dataset into 17 subsets D Y , with Y = 2004, …, 2020, according to the year of occurrence of the events: Here, t o;Y and t f;Y correspond to Y-01-01, 00:00 UTC and Y-12-31, 23:59 UTC, respectively.The cardinalities |D Y | of these yearly datasets are of order 10 4 events.Some years, for example 2011, contain a single large event.For each of these years, an additional dataset D 0 Y was analyzed corresponding to the D Y dataset trimmed by excluding the events occurring after the large event i Y and in a 10 ∘ radius around it: Finally, a magnitude-dependence analysis was carried out.To this purpose, events were divided into three subsets by The choice of the intervals aims at approximately balancing the size of the three subsets, thus compensating the effect of the power-law magnitude distribution of earthquakes.
Distance metrics.The fractal dimension ν was estimated through the evaluation of the correlation dimension by applying the Takens estimator 16 .As for other correlation dimension estimators, the Takens estimator relies on the computation of pairwise distances between events in a given embedding space.Results are expected to be independent of the metric used to compute distances.
The following four metrics, henceforth referred to as spherical surface (SP), flat 2-d (F2), flat 3-d (F3), and Euclidean 3-d (E3) were considered: where R 0 is the Earth radius.Longitude differences Δλ are computed according to the following rule: The variables x i , y i , z i are computed by mapping spherical coordinates to Cartesian coordinates, namely Finally, the angular depth ξ i is computed as 180 . The SP and F2 metrics are both acting on points that lie on a two-dimensional manifold-the surface of an R 0 -radius sphere or a twodimensional flat space with periodic boundary conditions, respectively-whereas the E3 and F3 metric act in a threedimensional space.Consequently, the metrics SP and F2 provide dimension estimates bounded between 0 and 2, while the E3 and F3 metrics provide values between 0 and 3.
The choice of the E3 metric follows from its immediate interpretation as the physical distance between points in threedimensional space (for example, it is the distance that seismic waves have to travel within the Earth's interior).On the other hand, because epicenters are located on a spherical surface, another natural choice is the SP metric, which is the geodetic distance between two points on the Earth surface.The F2 metric is chosen as a mean of comparison because it includes information both on latitude and longitude (contrarily, for example, to Chebyshev metric) and because it approximates the SP metric at low distances (contrarily, for example, to the taxicab metric).Finally, the F3 metric corresponds to an extension of F2 to three dimension by means of the inclusion of an additional angular quantity that incorporates information on depth.
Dimension estimation.Given a set of events, the Takens estimator is implemented as follows.
1. Pairs of events i, j are randomly extracted with replacement: a given pair is accepted only if the two events are uncorrelated 60,84 (see Section "Details on avoiding correlated events").For each accepted pair, and given a metric, the spatial distance d ij is computed.The extraction is repeated until N P accepted distances are collected.2. The set of N P accepted distances is ordered.The set can be expressed as {r k } with k 2 ½1; N 0 P , where N 0 P ¼ N P unless at least two equal distances d ij exist: r k thus corresponds to the k-th smallest distance d ij .3. The following function of the distance r is computed: where k runs from 2 to N 0 P .4. An interval of r, henceforth referred to as a plateau, on which the curve u(r k ) is constant is searched for, as explained in Section "Details on plateau assessment".A successful search, corresponding to the identification of a plateau, yields an estimate of fractal dimension νn and its uncertainty δν n . 5. According to a standard bootstrapping procedure, steps 1-4 are repeated N B times, thus yielding a set fν n g of estimated dimensions.On this set, the average ν, the standard deviation σ ν , the median ν, the 15th percentile P 15 ν , the 85th percentile P 85 ν , are then computed.To ensure the reliability of the method, the number of extracted pairs has to be much less than the number of available ones: N P ≪ D 2 /2, where D is the size of the set considered 60,85 .For global estimations, where D 0 2 /2 ≃ 3 ⋅ 10 10 , N P = 5 ⋅ 10 4 .For local estimations, which operate on sets of minimum size 500 corresponding to D 2 /2 ≃ 10 5 , we set N P = 2 ⋅ 10 4 .The number of bootstrap iterations was set to N B = 10 3 and N B = 10 2 for global and local estimates, respectively.Details on avoiding correlated events.The estimation of correlation dimension requires the evaluation of distances between pairs of points that are to be uncorrelated with each other 84 .This condition of uncorrelation is required by any method that relies on the evaluation of the correlation integral out of a set of points, regardless of the physical origin of the set itself 84 .Typically, this requirement reflects into the rejection-during the random extraction procedure-of those pairs of points that are deemed to be correlated according to some practical criterion, e.g., their closeness in time with respect to the typical autocorrelation time of the underlying process 60,[84][85][86] .
In the present work, the selection criterion applied to earthquake events is based on the identification of clusters of earthquakes according to the window-based procedure by Gardner and Knopoff 87 .The procedure considers a time interval w(m) and a distance d(m) that depend on magnitude.Given an event E i , all subsequent events within a time window w(m i ) and occurring within a distance d(m i ) from E i are deemed to be aftershocks of E i and thus assigned to the same cluster.Thereupon, the space-time window-based collection of aftershocks is applied again to the event in the cluster having the largest magnitude (if it differs from E i ).The parameters w(m), d(m), originally provided as tabulated values 87 , were here computed by means of the interpolating formulas proposed in ref. 88 .
The application of this clustering procedure allows to assign each event to a cluster (possibly containing a single earthquake).The selection criterion then consists in discarding a pair of events E i , E j if they belong to the same cluster.It is worth highlighting that the criterion does not correspond to the selection of mainshocks only, because aftershocks are retained in the dataset.The estimation of dimension neglects pairs of aftershocks belonging to the same cluster, but includes pairs belonging to different clusters.In this sense, the present correlation criterion allows to study the fractal geometry contributed by both mainshocks and aftershocks, while avoiding spurious results due to the presence of correlated distances 84 .
An analysis on the role of the cluster identification method is reported in Supplementary Note 1 (Supplementary Figs.3-4).Other two cluster identification methods were considered, namely the one proposed by Reasenberg 89 , and the one proposed by Baiesi and Paczuski 90 and further improved by ref. 91 .As shown in Supplementary Note 1, the three methods yield equivalent results in terms of global fractal dimension (Supplementary Fig. 5).Because the Gardner-Knopoff method provides the most conservative approach for the purpose of the present analysisi.e., it identifies larger clusters-it was selected as the cluster identification method throughout the work.
Details on plateau assessment.As described in the previous section, the Takens estimator implemented in the present work is based on the assessment of the curve u(r) given by Eq. ( 4).The curve u(r) should ideally converge to the value ν of the correlation dimension 16 .In order to identify this convergence, i.e., to obtain an estimate of the dimension ν, a possible approach is to look for a plateau within the curve, namely a region of r on which u(r) is approximately constant.
As explained above, a curve u(r) is sampled by a set of N 0 P values, with N 0 P at least equal to 2 ⋅ 10 4 .For each k value between ℓ/2 and N P /2 − ℓ/2 + 1, we consider the running window centered in k and encompassing ℓ = 10 3 point from k The choice of the maximum value of the center k corresponds to limiting the search only at distances below the median distance.This choice is justified by the fact that, in the standard framework of correlation dimension estimation, the power-law scaling of the correlation integral only holds for small distances 59,60 .
Given a running window, a linear regression is carried out on the corresponding ℓ points: the slope β and intercept α of the regression are computed as where all averages run on the (r k , u k ) pairs contained in the given window.In order to estimate confidence intervals, the following quantities are also computed: The 95% confidence intervals for α and β are estimated by considering the value t 0.025,ℓ−2 such that the cumulative t-distribution with ℓ − 2 degrees of freedom equals 0.975 (a two-tailed confidence interval is considered).In detail: The 95% confidence intervals for slope and intercept are finally determined as Among all the available windows, only those such that the slope is compatible-at 95% confidence-with zero, i.e., those for which 0 2 β À δ β; β þ δ β h i , are selected.If none is available, the search is deemed to be failed.Otherwise, the (up to) ten windows whose slope is compatible with zero and have the smallest β are considered.Within this set, the window having the smallest sample variance s 2 u , namely the one minimizing the fluctuations of u, is finally chosen as the plateau.The intercept α associated to this plateau is taken as the estimate ν.
Earth surface partition.The present work deals with both global and local estimates of fractal dimension.In the global case, events on the whole Earth were included in the estimation of dimension, thus yielding a single, global evaluation of ν.On the other hand, local estimates are inferred from events occurring within small contiguous subsets into which the Earth surface was partitioned and that are henceforth referred to as sectors.
We consider here partitions in 648 sectors.Let θ o , λ o be the coordinates (latitude and longitude) that identify the origin of a partition Π θ o λ o .Each sector belonging to Π θ o λ o has a size 10 ∘ × 10 ∘ and is defined as follows: where the indexes h, l take on values in the intervals [ − 9, 8], [ − 18, 17], respectively.Only the sectors that contain a sufficient number of events were considered for the evaluation of fractal dimension: the sufficiency threshold was set to 500.
To improve spatial resolution, we considered here four partitions, namely Π 00 , Π 05 , Π 50 , Π 55 .The set of all intersections of the sectors of these partitions corresponds to another partition Π 0 of the Earth surface consisting of 2592 subsectors each having a size 5 ∘ × 5 ∘ .A diagram illustrating this partitioning is shown in Fig. 10.
Let Σ be the 5 ∘ × 5 ∘ subsector given by the intersection of four sectors S 1 ; ; S 4 , each belonging to one of the four partitions, and let νðS k Þ; σ 2 ν ðS k Þ be the mean and variance estimated as in step 5. above in the case of the sector S k .The dimension ν corresponding to the subsector is then computed as the weighted average of the four mean values ν k ¼ νðS k Þ, where the weights are given by the inverse of the variances It is worth noting that, occasionally, the sum includes less than four terms because some of the sectors contributing to it might not provide an estimate of ν k due to a low event count or a failed estimate.The uncertainty δν on ν is evaluated as Once the local ν and its uncertainty δν are estimated, the underlying geometry can be considered fractal only if these values are incompatible with an integer value.To address this point we considered the distance of ν from its closest integer value b νe, normalized by the standard deviation δν: if the resulting quantity ν À b νe=δν is greater than 2, i.e., if the distance of ν to the closest integer is greater than two standard deviations, the value ν is deemed to be noninteger.

Fig. 1
Fig. 1 Global fractal dimension estimate.Normalized histograms of the fν n g sets (see Methods) of the estimated global fractal dimension of epicenters and hypocenters for different metrics.Points and error bars correspond to the median ν and the percentile ranges ½P 15 ν ; P 85 ν of each histogram.The black, dashed line corresponds to the dimension ν = 1.

Fig. 2
Fig. 2 Local fractal dimension estimates and comparison between epicenters and hypocenters.Maps of estimated local fractal dimension ν on the surface of the Earth and related histograms of ν for epicenters (a, b) and hypocenters (c, d), obtained by means of the E3 metric.In b, d, the vertical solid lines and shaded areas correspond to the median ν and the percentile range ½P 15 ν ; P 85 ν of the global estimate, respectively represented by the blue and red points and error bars of Fig. 1.The black, dashed lines highlight the integer values ν¼ 1; ν ¼ 2.The gray portions of the histogram bars count the subsectors for which the estimate of fractal dimension is deemed to be compatible with an integer value (see Methods).In a, c, the same subsectors are also reported in gray.e Map of the difference Δν ¼ ν hyp;E3 À ν epi;E3 of local dimension between hypocenters and epicenters.f Scatter plot of the values of local fractal dimensions for hypocenters versus those for epicenters: each data point corresponds to a subsector in which both assessments are available.The black, dashed lines highlight integer dimensions.In panels e, f the color scale is the same: red hues denote a higher dimension for hypocenters, while blue hues correspond to a higher dimension for epicenters.In a, c, e the Earth surface is mapped via an equirectangular projection.

Fig. 3
Fig. 3 Local fractal dimension estimates for epicenters obtained with the SP metric.Maps of estimated local fractal dimension ν on the surface of the Earth (a) and related histogram of ν (b).The color scale is the same in both panels.The Earth surface is mapped via an equirectangular projection.In b the vertical solid lines and shaded areas correspond to the median ν and the percentile range ½P 15 ν ; P 85 ν of the global estimate, namely the gray point and error bar of Fig. 1.The black, dashed lines highlight the integer values ν = 1, ν = 2.The gray portions of the histogram bars contain the estimates of the fractal dimension that are deemed to be compatible with an integer value (see Methods).

Fig. 4
Fig. 4 Local fractal dimension estimates for epicenters obtained with the F2 metric.Maps of estimated local fractal dimension ν on the surface of the Earth (a) and related histogram of ν (b).The color scale is the same in both panels.The Earth surface is mapped via an equirectangular projection.In b the vertical solid lines and shaded areas correspond to the median ν and the percentile range ½P 15 ν ; P 85 ν of the global estimate, namely the orange point and error bar of Fig. 1.The black, dashed lines highlight the integer values ν = 1, ν = 2.The gray portions of the histogram bars contain the estimates of the fractal dimension that are deemed to be compatible with an integer value (see Methods).

Fig. 5
Fig. 5 Local fractal dimension estimates for hypocenters obtained with the F3 metric.Maps of estimated local fractal dimension ν on the surface of the Earth (a) and related histogram of ν (b).The color scale is the same in both panels.The Earth surface is mapped via an equirectangular projection.In b, the vertical solid lines and shaded areas correspond to the median ν and the percentile range ½P 15 ν ; P 85 ν of the global estimate, namely the green point and error bar of Fig. 1.The black, dashed lines highlight the integer values ν = 1, ν = 2.The gray portions of the histogram bars contain the estimates of the fractal dimension that are deemed to be compatible with an integer value (see Methods).

Fig. 6
Fig. 6 Time-dependent global fractal dimension.Global fractal dimension as a function of years computed by means of the E3 metric for epicenters (a) and hypocenters (b).Points and error bars correspond to the median ν and the percentile ranges ½P 15 ν ; P 85 ν .The horizontal dashed lines and shaded areas correspond to the median values ν and the two percentile ranges ½P 15 ν ; P 85 ν (darker shade), ½P 05 ν ; P 95 ν (lighter shade) of the global estimates for the whole dataset (see also Figs. 1, 2).The black, dashed lines correspond to the dimension ν = 1.The difference between the original and trimmed values is thoroughly discussed in the main text.

Fig. 8
Fig.8Magnitude-dependent global fractal dimension.Global fractal dimension of the three magnitude ranges for both epicenters (blue) and hypocenters (red), computed by means of the E3 metric.Points and error bars correspond to the median ν and the percentile ranges ½P 15 ν ; P 85 ν .The horizontal lines and shaded areas correspond, as in Fig.6, to the medians ν and the percentile ranges ½P 15 ν ; P 85 ν of global fractal dimension evaluated on the whole dataset.

Fig. 7
Fig.7Dependence of global fractal dimension on number of events.Scatter plot of the global fractal dimension of epicenters (a) and hypocenters (b), obtained by means of the E3 metric, versus the number of events.Each dot and the related error bar correspond to the value of one year among those shown in Fig.6, and is reported with the same color: gray dots correspond to the years for which the fractal dimension was computed upon exclusion of large events (see main text).

Fig. 9
Fig.9Yearly number of events having m≥ 4.5.The different colors, as well as the separating black, dashed line, highlight the two distinct regimes that can be inferred from the data: gray, 1980-2003; green, 2004-2020.

Fig. 10
Fig. 10 Earth partition diagram.Each of the four 10 ∘ × 10 ∘ squares bounded by solid lines corresponds to the seed sector of the respective partition S 00 ; S 05 ; S 50 ; S 55 (origins are marked by crosses).The point P belongs to the 5 ∘ × 5 ∘ subsector that results from the intersection of the four sectors. ∑