Micrometre-scale deformation observations reveal fundamental controls on geological rifting

Many of the world’s largest volcanic eruptions are associated with geological rifting where major fractures open at the Earth’s surface, yet fundamental controls on the near-surface response to the rifting process are lacking. New high resolution observations gleaned from seismometer data during the 2014 Bárðarbunga basaltic dyke intrusion in Iceland allow us unprecedented access to the associated graben formation process on both sub-second and micrometre scales. We find that what appears as quasi steady-state near-surface rifting on lower resolution GPS observation comprises discrete staccato-like deformation steps as the upper crust unzips through repetitive low magnitude (MW < 0) failures on fracture patches estimated between 300 m2 and 1200 m2 in size. Stress drops for these events are one to two orders of magnitude smaller than expected for tectonic earthquakes, demonstrating that the uppermost crust in the rift zone is exceptionally weak.


Experiment and Data Analysis
In the afternoon of 30 August 2014, we installed a small profile of three 3-component broadband seismometers (Guralp 6TD 30 s) perpendicular to the graben and inferred dyke (Fig. 1a), with the closest station (DY3) directly at the western shoulder of one of the large graben boundary faults and the other two stations approximately 1 km (DY1) and 2 km (DY2) from DY3. The surrounding area was characterised by several metres of poorly consolidated volcanic ash and sand on top of partially fractured basaltic lava flows 13 , a strongly scattering environment for seismic waves. As strong ground shaking could be felt during the experiment, the operation had to be aborted for safety reasons, resulting in ~26 minutes of synchronous data on all stations. On 5 September 2014, a new fissure opened approximately 600 m east of DY3 and effused lava for 2 days. b a 7,196 7,195 7,194 7,193 7,192 7,191 7,190 7,189 406 408 410 412 4 14 Northing (  The unprocessed vertical velocity seismograms (Fig. 1b) show coherent activity on all three stations. However, the focus of this study lies in five high amplitude signals on station DY3 (red arrows in Fig. 1b), which are not registered by the other stations DY1 and DY2, suggesting that the causative events are relatively small and local to DY3. The velocity seismograms and the corresponding scalograms of these events (event number 3 shown in Fig. 2a,b) show impulsive waveforms with a main frequency peak between 3 and 8 Hz and a secondary peak above 25 Hz. As high-frequency waves are attenuated strongly when travelling through the ground, such high frequencies thus indicate a fracturing process close to the station. In a recent study 14 , we presented a new data processing approach that allows for the recovery of micrometre-scale displacement steps from instrument-corrected seismograms. It is based on long-period noise removal using median filters and its performance was confirmed by laboratory experiments. We apply this method to the event (Fig. 2c-e) and observe displacement steps on all three (orthogonal) components of the instrument, i.e. the station was displaced by approximately 125 μ m in a northwest, slightly upward direction. This represents a motion away from the centre of the graben and the underlying dyke. Applying this procedure to the full-length records on this station reveals similar amplitudes and ratios between different displacement components for all five events (Fig. 2f), suggesting a repetitive process with similar source locations and an apparent average inter-event time of about 4.5 minutes. As horizontal components of seismometers are also susceptible to ground rotation [15][16][17] , possible tilts can be estimated from the data using the tilt transfer function [18][19][20] ; this involves a simple integration of raw data and multiplication with a factor depending on well-known instrument properties. The resulting traces (Fig. 2g) show tilt steps of 1.3-4 μ rad oriented in a northwest direction associated with each of the five events. The tilt step directions coincide with the direction of the displacement steps ( Fig. 2h) and support a repeating source process generating the events, with roughly consistent amplitudes and locations.
Source location. Although one station is not sufficient to fully invert for source locations and mechanisms, we use the observed static deformations from DY3 to explore potential sources with a forward modelling approach. We estimate the source location and magnitude by (i) assuming a plausible source mechanism and (ii) performing a search over a 200 × 200 × 100 m 3 grid around the station, where we match the observed ratios between different deformation components with the theoretical values for a homogeneous, elastic medium 21 . The ratios are defined as where u Z , u N and u E are displacements and t N and t E are tilts. Subscripts Z, N and E denote a vertical, north and east direction, respectively. These ratios are used to compute two misfits, R d (displacements only) and R dt (both displacements and tilts), defined as: The minimum misfits indicate the best location for the chosen source mechanism and the corresponding seismic moment M 0 can be found by a simple least squares inversion.
As the station is located in direct proximity to the faults associated with the graben formation, we suspect that the local events are part of the faulting process. Consequently, we assume a 75° dip-slip mechanism parallel to the N25°E striking boundary fault 9 as the source mechanism. Figure 3 shows the misfits for this normal fault mechanism, where the observed tilt and displacement values of the third step (Fig. 2f,g) are used. Here we assume a medium P-wave velocity of V P = 500 m/s and Poisson's ratio of v = 0.3, consistent with values obtained for unconsolidated upper geological layers at various volcanoes [22][23][24][25][26] (further discussion in supplementary information). For clarity, only misfit values below 0.5 are displayed and all remaining misfits are located in the quadrant south-east of the source. As the displacement-only misfit R d (Fig. 3a) does not converge around a single minimum, it can only indicate the approximate direction of the source with respect to DY3. When tilts are introduced (Fig. 3b), a   Table 1). The source moments are small enough to justify the use of the point-source assumption in our forward modelling approach. Static displacements such as those observed at DY3 are near-and intermediate-field effects and can only be observed within a fraction of a wavelength from the seismic source 27 . The sources inferred above would theoretically cause total static displacements smaller than 1 μm at the next closest station, DY1. Sub-micrometre steps are not detectable with our instruments and methods 14 . The fact that the events are not visible at the other stations also implies that the dynamic seismic signals, i.e. all near-, intermediate-and far-field components, fall under the noise level at these locations, likely due to strong wave attenuation in the unconsolidated surface materials 13 . Source dimensions. For the source location found above, we determine source parameters (size and slip) by removing path effects through deconvolving modelled deformation and seismic radiation (Green's function) from the recorded seismogram shown in Fig. 2c. Here we use the linear relationship between the ground displacement spectrum U(ω) and the source moment spectrum M(ω) 27 : P where the Green's functions G depend on the receiver position r relative to the source, the elastic properties of the medium V P and ν, the density ρ, the radiation pattern RP for a specific source mechanism and the quality factor Q. This simplifies the deconvolution to a simple division M = U/G for each frequency. The resulting source moment spectrum M(ω) is subsequently fit with a Brune ω 2 source model 28 in order to determine the corner frequency. We calculate G(ω) for the inferred normal fault source with the expressions given by Aki and Richards 29 and modified by Lokmer and Bean 27 , using the same parameters as above (V P = 500 m/s, ν = 0.3). Q is varied until we obtain the best fit to the ω 2 -model (Q = 20). The source-time history M(t) resulting from this deconvolution is shown in Fig. 4a. Its spectrum and the ω 2 -model fit are shown in Fig. 4b, resulting in a corner frequency of 4.5 Hz. This frequency is used to determine the source size (and subsequently the slip D using M 0 = μAD, with the shear modulus μ and the source area A): approximating the source as a slipping circular patch 30,31 gives a source radius of approximately 10-20 m with an average slip of 1-4 mm. As the actual source mechanism cannot be inferred from our data and a tensile component could potentially form part of the source process, we additionally performed the location grid search and source slip analysis for a tensile crack (Supplementary Figure 1 and Supplementary Table 1). If a purely tensile source mechanism is considered, the slip displacement on the same patch is reduced by a factor of 2, showing the results are robust for a deviation from the pure normal faulting source. Both results are in agreement with Liu-Zeng et al. 32 , who model the slip-to-length ratio and obtain equivalent values for small faults with rough fault surfaces.
We estimate stress drops of Δ σ = 0.008-0.07 MPa using Δ σ = 7M 0 /16r 3 according to Eshelby 33 . These values are 2-3 orders of magnitude smaller than expected for tectonic seismicity (stress drops typically > 1 MPa 34 ) and point to a very weak uppermost crust. They are consistent with the lack of shallow "standard" earthquakes associated with such a large rifting event. Such small stress drops are in striking agreement with the value of Δ σ = 0.01 MPa obtained for shallow long period seismicity on Mt Etna, Italy 35 , attributed to the presence of exceptionally weak near surface volcanic material that could not sustain high shear or tensile stresses and hence also failed at exceptionally low earthquakes magnitudes.

Discussion
Our data reveal new information about the rifting process, suggesting that it is at least partially discrete, occurring in micrometre scale steps. This raises questions about how these displacements compare to the observed long-term deformation in the area. Combining the time-history of the closest GPS stations with the total graben opening measured from satellite data (see supplementary information), the deformation rate for 30 August 2014 is estimated to be roughly 5 cm/day. Assuming a repeating process with average displacement steps of 133 μ m and average inter-event times of 267 s, observed at DY3, we extrapolate our data and obtain an approximate deformation rate of 4.3 cm/day. Furthermore, accumulating normal fault slip estimates at the source of 1-4 mm yields a horizontal deformation rate of 7-27 cm/day. Although the modelled slip values are approximate, both of our displacement measures are in good agreement with the GPS estimates. The similarity suggests that the satellite and GPS-derived long-term surface deformation associated with Earth surface rifting is a consequence of displacement accumulated through very low magnitude discrete brittle failure at the millimetre scale. The detection of such steps is limited to distances within a few hundred metres from the source, highlighting the rarity of such observations. The similarity also suggests that any aseismic component is small at the spatial and temporal scales captured in this study; it also indicates that fracturing of the weak uppermost crust is limited to microseismic events, consistent with the lack of observed shallow seismicity 4,7 . We conclude that at its smallest temporal and spatial scales, rifting in the uppermost Earth's crust is not a steady state process but rather exhibits transient staccato-like behaviour that yields definable spreading rates only when viewed over longer time scales. Stress drop analysis on the discrete micro-events reveals that the uppermost crust is exceptionally weak in the rift zone.