Dike intrusions during rifting episodes obey scaling relationships similar to earthquakes

As continental rifts evolve towards mid-ocean ridges, strain is accommodated by repeated episodes of faulting and magmatism. Discrete rifting episodes have been observed along two subaerial divergent plate boundaries, the Krafla segment of the Northern Volcanic Rift Zone in Iceland and the Manda-Hararo segment of the Red Sea Rift in Ethiopia. In both cases, the initial and largest dike intrusion was followed by a series of smaller intrusions. By performing a statistical analysis of these rifting episodes, we demonstrate that dike intrusions obey scaling relationships similar to earthquakes. We find that the dimensions of dike intrusions obey a power law analogous to the Gutenberg-Richter relation, and the long-term release of geodetic moment is governed by a relationship consistent with the Omori law. Due to the effects of magma supply, the timing of secondary dike intrusions differs from that of the aftershocks. This work provides evidence of self-similarity in the rifting process.

T here is a general consensus that earthquakes exhibit self-similar behavior over a wide range of magnitudes on both local and global scales. Tectonic earthquakes satisfy two power-law scaling relationships to a reasonable approximation: the Gutenberg-Richter relation 1,2 , which describes the frequency-magnitude distribution, and the modified Omori law 3,4 , which describes the temporal decay of aftershocks. Driven by the accumulation of tectonic stresses, the earthquake rupture process is controlled by local stress conditions and the frictional and material properties of faults.
Self-similarity is also observed in volcanic systems. Both the global rate of volcanic activity 5 and the frequencyvolume distributions of major volcanic eruptions 6 have been shown to follow power laws. Most volcanic eruptions are initiated by dike intrusions, which nucleate when the pressure inside an inflating magma reservoir exceeds a threshold value. Once a dike is nucleated, propagation is controlled by the reservoir pressure as well as viscous, elastic and thermal stresses 7 . Magma replenishment is necessary for additional intrusions to be sourced from the same reservoir 8 .
Along late-stage continental rifts, plate boundary separation is accommodated primarily by dike intrusions, and to a lesser extent by faults, during discrete rifting episodes 9 . Tensile stresses accumulated over decades to centuries exert the dominant controls on the deformation cycle, although crustal magma accumulation influences the timing of individual dike intrusions [9][10][11][12][13] . Discrete rifting episodes have been observed along two subaerial divergent plate boundaries 9 -the Krafla segment of the Northern Volcanic Rift Zone in Iceland (1975Iceland ( -1984 and the Manda-Hararo segment of the Red Sea Rift in Ethiopia (2005)(2006)(2007)(2008)(2009)(2010). In both episodes, the spatial and temporal patterns of dike intrusions share many similarities with tectonic earthquakes in mainshock-aftershock sequences. The initial and largest dike intrusions were followed by clusters of smaller dike intrusions. Additionally, in a process that closely resembles the triggering of aftershocks by static stress changes, later dike intrusions were emplaced in areas where local tectonic stresses had been either increased or not completely relieved by earlier dike intrusions 13,14 .
In the East African Rift System, the spatial distributions of earthquakes and volcanic vents are known to follow power laws 15,16 . Here, we examine the scaling relationships for dike intrusions associated with the Krafla and Manda-Hararo discrete rifting episodes. We perform a detailed statistical analysis of the dimensions and timings of individual dike intrusions and assess whether their behavior follows power laws analogous to the Gutenberg-Richter relation and the modified Omori law. The Krafla rifting episode is composed of 20 major dike intrusions associated with subsidence of the Krafla caldera. Accounting for the fact that reservoirs expand as their internal pressures decrease, the volume of each dike intrusion is estimated to be twice the subsidence volume measured from elevation and gravity data (Fig. 1d) 17,18 . The lengths of the Krafla dike intrusions are constrained by field measurements of fractures and the hypocenters of triggered earthquakes 17,18 .
The dimensions of the 14 dike intrusions of the Manda-Hararo rifting episode are constrained well by seismic and geodetic observations 13,14,19 . We use the MHA 19 (and I. Hamling personal communication 2011), and MHB 13 datasets for estimates of the lengths and volumes of each dike intrusion. The MHB dataset does not include information for the most recent intrusion, but does provide estimates of the dislocated areas of the first thirteen dike intrusions. For the volume of the first intrusion in the MHB dataset, we exclude contributions from magma chambers beneath Dabbahu and Gabho volcanoes, and only consider the volume sourced from the magma reservoir associated with the Ado'Ale Volcanic Complex (Fig. 1a), which was active during the entire rifting episode 9,20 .
In addition to examining the dimensions of the dikes associated with each rifting episode separately, we also combine the volume information to produce two global rifting datasets -KMHA (Krafla and MHA) and KMHB (Krafla and MHB). where M C is the seismic moment corresponding to the magnitude of completeness for a given seismicity catalog and b 5 2/3 for the typical situation where the b-value is one 22,23 . The probability that an earthquake meets or exceeds a given size is therefore described by the survivor function F(M 0 ) 5 (M C /M 0 ) b . Because seismic moment is proportional to rupture area, if earthquakes can be approximated by slip on circular fault patches, similar power-law relationships can be defined for rupture area and length 24 .
The size of a dike intrusion is often described by the geodetic moment, which is proportional to the volume of magma intruded in an opening tensile crack. If dike intrusions along divergent plate boundaries obey self-similar scaling relationships like earthquakes, and if dike intrusions can be approximated by circular cracks, then the lengths, dislocated areas and volumes of dike intrusions should all follow power law distributions.
Here, we examine each of our rifting datasets for power-law behavior. Assuming that dike intrusions obey self-similar scaling www.nature.com/scientificreports SCIENTIFIC REPORTS | 4 : 3886 | DOI: 10.1038/srep03886 relationships above some threshold size, x min , we model each survivor function as N( is the probability that a dike intrusion meets or exceeds a given size and x represents either the volumes, dislocated areas, or lengths of dike intrusions in a given rifting dataset. We estimate the parameters b and x min following Clauset et al 25 and assess the goodness-of-fit using the Kolmogorov-Smirnov statistic. Using the value of x min determined from the power-law parameterization, we also model each rifting dataset with exponential and log-normal distributions and use the likelihood ratio test to determine which right-skewed distribution provides the best fit to the data (see Method and Tab. 1).
For the two global rifting datasets (KMHA and KMHB), we find that intrusion volumes are modeled well by power laws and poorly by exponential and log-normal distributions (Fig. 2). In most cases, power laws also outperform the alternative distributions for individual volume, length and area datasets (Tab. 1 and Fig. 3), although the small sample sizes preclude a statistically robust comparison. Nonetheless, our observation that a power law consistently performs well across the various rifting datasets suggests that dike intrusions have self-similar scaling relationships.
Power-law parameterizations for the global rifting datasets are characterized by b-values of ,1.3 to ,1.5, which are significantly higher than the theoretical value of b 5 2/3 for earthquakes as well as the empirical value of b found for most tectonic earthquake sequences 2,26,27 . Earthquake sequences with elevated b-values (i.e. 1 , b , 2), which have a larger proportion of small events, have been observed along divergent plate boundaries as well as in volcanic and geothermal environments 23,28,29 , and have been attributed to low values of confining stress 30 , material heterogeneity 31 , large thermal gradients 32 , and the presence of fluids 33 .
We note that the value of x min for both global rifting datasets is ,40310 6 m 3 to ,60310 6 m 3 . Due to differences in how dike intrusions were detected at Krafla and Manda-Hararo, however, it is difficult to determine if x min results from catalog incompleteness 34 or whether this quantity is associated with some threshold property of the source magma chambers, such as a minimum volume or pressure required for dike nucleation 7,35 .
Individual dike intrusions of the Krafla rifting episode were identified by rapid subsidence of Krafla caldera caused by the deflation of a 3-5 km deep magma chamber 17 . There, we note that the volume rifting datasets (Fig. 3c) could be fit by two power laws with different values of b corresponding to intrusions smaller and larger than x min .    The smallest intrusions occurred later in the sequence within 10 km of the center of the caldera (Fig. 1d), and were accompanied by fissure eruptions. Our observation of a potential piecewise power-law relationship may be due to the fact that lateral propagation played a lesser role in the later dike intrusions, or because these events were controlled to a greater degree by local stresses rather than tectonic extension 9 . At Manda-Hararo, individual dike intrusions were identified by local seismic and satellite geodetic data rather than monitoring the 8-10 km deep source magma chamber directly 9 , and we find that the volume, length and surface area rifting datasets are well modeled by a single power-law relationship. The values of x min for these datasets are close to the smallest dimensions observed, although we cannot exclude the possibility that additional small dike intrusions are missing from the catalog since the resolution of geodetic data decreases with depth and the maximum magnitude of earthquakes induced by dike intrusions is correlated with the volume of magma injected in the rift 36 .
A modified Omori law for rifting episodes. For earthquakes, the decay of aftershock activity with time follows an empirical relationship known as the modified Omori law 3 , n(t)~K (tzc) p , where n(t) is the rate of aftershocks in a given magnitude range occurring at time t since the mainshock, and K, c, and p are constants. Assuming that aftershocks follow the Gutenberg-Richter relation and the modified Omori law, the cumulative seismic moment released in a time interval can be expressed as the product of the average seismic moment released per event and the total number of events in the time interval, the integral form of the modified Omori law. It follows therefore that the seismic moment release rate also exhibits powerlaw time dependence 37,38 . These simple relationships result from selfsimilarity in the processes associated with static and dynamic earthquake triggering 39,40 , and deviations are often linked to the occurrence of strong aftershocks that generate secondary sequences of dependent events 3 . A modified Omori law has also been identified for earthquakes preceding and following volcanic eruptions, which suggests that the brittle crust may respond similarly to tectonic and volcanic stresses 41 . Here, we evaluate whether the rate of volume emplaced in the crust during the rifting sequences can be modeled with power laws similar to the modified Omori law for earthquakes.
Analogous to earthquakes in a mainshock-aftershock sequence, the long-term volumetric rate during rifting episodes decays with time whereas the time between consecutive dike intrusions increases with time (Fig. 1). Secondary diking sequences are also observed  following the largest dike intrusions. We quantify this behavior by calculating rifting volumetric rates using the equation v r 5 K/t p , where t is the time from the onset of a rifting episode, p is the power-law exponent and K is a constant (Fig. 4). From a simple linear regression in logarithmic space, we find that ,40-60% of the variance in the datasets can be explained by power-law relationships with p-values of 0.73 for Krafla and 0.84-0.89 for Manda-Hararo. These parameterizations closely resemble the form of the modified Omori law, which is predicted by rate-state theory for any impulsive stress change 42 . The relatively short emplacement times of the dike intrusions, which were on the order of a few hours to days, compared to interevent times of a few months to years may explain why we find consistency between the temporal behavior of earthquakes and dike intrusions in rifting episodes. For earthquake sequences that follow the Gutenberg-Richter relation and the modified Omori law, it has been shown that aftershock productivity is a function of mainshock magnitude 43 . If rifting episodes are governed by similar physics, this could explain why the initial dike intrusions of the Krafla and Manda-Hararo sequences, which had volumes larger than ,100310 6 m 3 , were followed by several smaller diking events. Sequences with multiple dike intrusions and large cumulative volumes have also been recognized during past rifting episodes in the Northern Iceland Volcanic Zone 17 . Conversely, single dike intrusions are associated with rifting episodes at Dallol in Ethiopia in 2004 and Lake Natron in Tanzania in 2007. The volumes emplaced during the Dallol and Lake Natron episodes, 60310 6 m 3 and 90310 6 m 3 respectively 44,45 , are comparable to the values of x min for the Krafla and Manda-Hararo rifting datasets.
At volcanoes, the interevent repose time, or the time between magmatic eruptions, is correlated with erupted volume 46 . In open conduit systems characterized by short repose periods, volcanic eruptions generally follow the time-dependent model, where interevent time is dependent on both the size of the last eruption and the magma recharge rate [47][48][49] . Here, we investigate the role of the magma recharge rate in controlling the timing of individual dike intrusions at Krafla and Manda-Hararo by examining the time intervals between successive dike intrusions in each rifting episode. We use the Bayesian hierarchical time-predictable model of Passarelli et al 49 , assuming that the interevent time, r i 5 t i 1 1 2 t i , follows the empirical equation where v i is the volume of the ith dike intrusion, and d and c are constants. If c 5 0 the system is not time-predictable, whereas if c 5 1 the magma recharge rate is constant and the interevent times and volumes of dike intrusions scale linearly 50 . We perform 10,000 simulations of the statistical model to determine the value of c for each rifting dataset.
We find that the average values of c are 0.84 for the Krafla dataset and 0.63 and 0.69 for the Manda-Hararo datasets (Fig. 5). These values of c , 1 demonstrate a non-linear relationship between dike intrusion volume and interevent time, such that the interevent times following large dike intrusions are longer than those following small dike intrusions, but shorter than those predicted by the classical time-predictable model. Because the inflation of magma reservoirs is observed to decay exponentially after dike intrusions 8,20,51 , this   observation suggests that the shallow magma chambers beneath Krafla and Manda-Hararo are hydraulically connected to larger and deeper reservoirs 9 that allow magma recharge to occur at faster rates for larger pressure drops. The residuals about the best-fit regression lines for the Krafla and Manda-Hararo datasets are larger than those obtained for similar analysis completed for Kilauea and Etna volcanoes 49 , which may suggest that magmatic stresses play a lesser role in controlling the timing of individual dike intrusions during rifting episodes relative to purely volcanic systems. Influenced by the magma recharge rate, the time-predictable behavior of dike intrusions at Krafla and Manda-Hararo seems to be in direct contrast to earthquakes, where the time interval between two events is inversely correlated with the magnitude of the first event, such that shorter interevent times are observed following larger earthquakes 52,53 .

Discussion
We have demonstrated that the Krafla and Manda-Hararo diking episodes follow self-similar scaling relationships analogous to the Gutenberg-Richter relation and the modified Omori law for tectonic earthquakes. High p KS -values demonstrate the goodness of fit of the best-fit power-laws at high significance levels, and low p LR -values prove that power-laws reliably outperform the best-fit log-normal and exponential distributions even for our small sample sizes. Our statistical analysis demonstrates that the dimensions of dike intrusions follow a power-law model with small levels of variability around the best-fit parameters. With the possible existence of a pressure or volume threshold for dike nucleation, the major processes influencing the dimensions and timings of individual dike intrusions appear to be the accumulation of tectonic stress caused by long-term plate motion, static stress transfer between intrusions, and magma recharge.
Along continental rifts as well as mid-ocean ridges, dike emplacement is controlled by tectonic stresses and the magma supply rate. In areas with limited magma supply, where the pressure inside a crustal magma chamber is reduced during dike propagation, a sequence of dike intrusions rather than a single diking event may be necessary to fully release the accumulated tectonic stress 12 .
If the first dike intrusion encounters a structural barrier, it may continue to widen rather than lengthen until the driving pressure reaches a critical level, resulting in an aspect ratio that differs from later, smaller dikes intruded at lower magma pressures. Analogously, the thickness of the seismogenic zone limits the downdip rupture width, resulting in different aspect ratios for small and large earthquakes. The largest earthquakes are known to deviate substantially from the typical Gutenberg-Richter relation 54 , which is similar to our observation that the largest dike intrusions from the Krafla and Manda-Hararo rifting episodes do not obey the scaling laws that fit smaller dike intrusions well. We suggest that this discrepancy may be the result of factors such as the limited magma supply rate, threshold pressure conditions for dike propagation that change over time, and the structural constraints of lithospheric thickness and rift segment length.
Given that we are currently limited to studying rifting episodes along subaerial divergent plate boundaries, and since rifting episodes at a given location have recurrence intervals on the order of hundreds of years, we are restricted to this relatively small dataset. Nonetheless, the power-law behavior identified here is surprisingly robust, and with time and improvements in geophysical instrumentation, it may be possible to evaluate whether these observations apply more broadly to divergent plate boundaries in general.

Methods
In order to estimate the parameters b and x min we follow the methodology presented by Clauset et al 25 and select values that jointly minimize the Kolmogorov-Smirnov (KS) statistics. Errors for b and x min are estimated from 1000 bootstrap iterations 25 . Given the best-fitting values of b and x min , we calculate 1000 synthetic power-law distributed datasets and their KS statistics. The p KS -value associated with each rifting dataset is obtained by counting the fraction of the time that the KS statistics for the synthetic datasets are larger than the KS statistic for the empirical dataset. We rule out a power law distribution for datasets with p KS -value # 0.1 25 .
We also model each rifting dataset using exponential and log-normal models lefttruncated at x min , with the functional forms f(x) 5 lexp(lx min )exp(2lx) and tively. Once the best-fit distributions are found and the fits are assessed by the KS test as for the power law model, we use the likelihood ratio test for model selection. The logarithm of the likelihood ratio, LR, has a positive sign if the power-law distribution provides a better fit to the data than the alternative distribution. The associated p LRvalue describes the statistical significance of LR, that is whether it is robustly far from zero. For large p LR -values, the likelihood ratio test cannot discriminate between the two distributions, however, for p LR -value # 0.1, the sign of LR is a reliable indicator of which model provides a better fit to the data 23 .
Analogous to the classical Omori-decay fitting procedure, we bin the time intervals of each dataset (40 bins for MHA and MHB and 45 bins for Krafla) in order to obtain the rate of volume emplacement in each time interval. Bins with zero events are discarded. We use classical linear regression analysis to solve the equation v r 5 K/t p for the best-fit parameters and their associated errors, which are reported in Fig. 4 along with the R 2 value. We perform an F-test on the slope of the regression line where the null hypothesis is that the slope p 5 0. The p-values are less than 0.02 in all three cases, which means we can reject the null hypothesis at the 95% confidence level.
We evaluate a time-predictable model for dike intrusions using the Bayesian hierarchical model developed in Passarelli et al 47 . We assume that the errors for the interevent times and volumes are one day and 25%, respectively. The posterior distributions of the c and d parameters are obtained by performing 10,000 simulations using the Metropolis-Hastings integration technique. Using these posterior distributions, we select the average values of the parameters as the best-fit values of the model and plot the 5th and 95th percentiles to demonstrate model variability.