LETTERS
Preseismic velocity changes observed from active source monitoring at the Parkfield SAFOD drill site
Fenglin Niu 1,Paul G.Silver 2,Thomas M.Daley 3,Xin Cheng 1&Ernest L.Majer 3
Measuring stress changes within seismically active fault zones has been a long-sought goal of seismology.One approach is to exploit the stress dependence of seismic wave velocity,and we have inves-tigated this in an active source cross-well experiment at the San Andreas Fault Observatory at Depth (SAFOD)drill site.Here we show that stress changes are indeed measurable using this tech-nique.Over a two-month period,we observed an excellent anti-correlation between changes in the time required for a shear wave to travel through the rock along a fixed pathway (a few microse-conds)and variations in barometric pressure.We also observed two large excursions in the travel-time data that are coincident with two earthquakes that are among those predicted to produce the largest coseismic stress changes at SAFOD.The two excursions started approximately 10and 2hours before the events,respect-ively,suggesting that they may be related to pre-rupture stress induced changes in crack properties,as observed in early labor-atory studies 1,2.
It is well known from laboratory experiments that seismic velocit-ies vary with the level of applied stres
s 3–5.Such dependence is attrib-uted to the opening and closing of microcracks due to changes in the stress normal to the crack surface 6–8.In principle,this dependence constitutes a stress meter,provided that the induced velocity changes can be measured precisely and continuously.Indeed,there were sev-eral attempts in the 1970s to accomplish this goal using either explos-ive or non-explosive surface sources 9–11.The source repeatability and
the precision in travel-time measurement appeared to be the main challenges in making conclusive observations.
With the availability of highly repeatable sources,modern data acquisition systems and advanced computational capability,Yamamura et al.12showed compelling evidence that seismic velocity along a baseline in a vault near the coast of Miura Bay,Japan,responds regularly to tidal stress changes.Silver et al.13found an unambiguous dependence of seismic velocity on barometric pressure from a series of cross-well experiments at two test sites in California.The stress sensitivity depends primarily on crack density and has a strong nonlinear dependence on confining pressure.Consequently,crack density is expected to decrease rapidly with depth,as should stress sensitivity.It is thus unclear whether the stress-induced velo-city variations observed at shallow depths 12,13are still detectable at seismogenic depth.
To explore stress sensitivity at seismogenic depth,we have con-ducted an experiment at Parkfield,California,where adjacent deep wells,namely the SAFOD pilot and main holes (Fig.1),are available.Accurately located seismicity,together with the availability of high-quality geophysical data in the Parkfield region,make it one of the best areas to detect temporal changes related to the earthquake cycle.A specially designed 18-element piezoelectric source and a three-component accelerometer were deployed inside the pilot and main holes,respectively,at ,1km depth (see Methods).The experiment was conducted for ,2months:the first period was 29October to 28
1
Department of Earth Science,MS-126,Rice University,6100Main Street,Houston,Texas 77005,USA.2Department of Terrestrial Magnetism,Carnegie Institution of Washington,5241Broad Branch Road,NW,Washington DC 20015,USA.3Earth Sciences Division,Lawrence Berkeley National Laboratory,1Cyclotron Road,Berkeley,California 94720,
USA.
35º 54′
35º 56′
35º 58′
36º 00′
36º 02′
Middle Mountain SAFOD
Pilot hole
?
1
2
3
a
Receiver
Source
D e p t h  (k m )
b
Figure 1|Map of the experiment site.a ,Circles show earthquakes that occurred during the experiment period.The M 53and M 51events are shown as red and green circles,respectively.Star indicates the Parkfield SAFOD drill site,where the experiment was conducted.Triangle,location of the Middle Mountain creepmeter;squares,locations of the Donalee (DLT)and Frolich (FLT)Gladwin borehole tensor strainmeters.b ,A vertical section (schematic)of the SAFOD main and pilot holes.Red vertical lines
indicate the source and receiver locations.Background image is electrical resistivity 18with blue (red)corresponding to relatively high (low)resistivity.White circles show the seismicity,red dashes indicate the SAFOD
instrumentation,and black short lines represent sub-horizontal holes drilled off the main hole.The profile at the top of the panel is the surface topography with the arrow indicating the surface exposure of the fault.
Vol 454|10July 2008|doi:10.1038/nature07111
204
November2005,and the second was11December2005to10January 2006.We fired a pulse with a width of1ms four times per second and recorded200-ms-long data with a sampling rate of48,000Hz.The waveforms were automatically stacked in groups of100shots,result-ing in one record(Fig.2)acquired every27s(two additional seconds were needed in storing the data).
To enhance the signal-to-noise ratio(SNR)of the data,we further stacked the raw seismograms in sets of100.This stacking procedure reduced the data to one stack every45min.The45-min stacked records were then processed with a bandpass filter of1–5kHz before the travel-time analysis.We used a cross-correlation-based method to estimate the delay time,which permits subsample precision(see Methods).No smoothing or filtering was applied to the measured delay time series.The error in delay time measurement was estimated to be,1.131027s,based on SNR analysis(see Methods),and this estimate was confirmed by comparing measurements from consec-utive recordings.As the nominal travel time of the shear-wave(S-wave)coda along the baseline is about10ms,the detectable threshold of velocity perturbation is,1.131025,
We measured the delay times of the S wave and the S wave plus its coda up to20ms with respect to a fixed reference trace for each period(Fig.3).The measurements show daily cycles that are well correlated with the temperature record(Fig.3).Silver et al.13found that this temperature sensitivity origi
nates from the electronics of the recording system rather than from changes in the subsurface velocity field.We excluded the measurements of the first few days to allow the source and sensor to be stabilized at their locations.We also removed the linear trend from the data as was done by Silver et al.13.In general, the delay times of the coda are about twice as large as those of the S wave,suggesting that they are caused by a change in the velocity of the bulk media,as the coda travels for a longer time in the media and thus is expected to accumulate a larger travel-time anomaly.The delay time closely follows the barometric pressure changes for the first period(Fig.3a).
After removing the temperature effect from the measured delay time variations(Fig.3a),we obtained a delay time change of,3.0m s in the first period.The corresponding velocity perturbation is about 331024,about an order of magnitude higher than the detectable threshold.During the same time period,the change in barometric pressure is,1.3kPa.We used a linear regression to estimate the stress sensitivity of the velocity and obtained a value of 2.431027Pa21.We also calculated the predicted solid Earth tides at the site in the same period and found that the tidal stress varies within240Pa,nearly an order of magnitude smaller than changes in barometric pressure.Thus the travel-time changes induced by tidal stress are of the order of1027s,close to the measurement error and thus are predicted to be undetectable.
The negative correlation between travel time and barometric stress can be further seen in the delay time data through to the end of the ninth day of the second period.After this time the relationship starts to break down,and we observe instead two prominent excursions in the delay time data that are not seen in the barometric pressure record.It is also confirmed that the two excursions were not caused by precipitation or instrumentation.The amplitudes of the two excursions are,5.5m s and,1.5m s,respectively,over the nominal ,10ms coda travel time.Using our measured stress sensitivity of 2.431027Pa21,the corresponding stress changes are2.3kPa and 625Pa for the first and second peak,respectively.
In order to evaluate the possibility of a tectonic cause for the excursions,we examined the seismicity around the SAFOD site occurring in the experiment period(Fig.4a).The first peak appears to correspond to the largest earthquake occurring in this period (date,24December2005;time,10:10:57.21(h:min:s UT);location, 35.9970u N,120.5565u W;depth,3.88km;magnitude3.00,hereafter the M53event),while the second peak corresponds to the second closest(1.5km)event to the experiment site(29December2005, 01:32:50.87UT,35.9788u N,120.5397u W,depth1.82km,magnitude 0.98,hereafter the M51event).The closest event is about1.3km away from our site,but its size is only M50.34and thus should not have a large effect at the site.
We calculated the predicted static stress change at SAFOD assoc-iated with these two earthquakes.The near-field static displacement at a location r with respect to the earthquake is proportional to M o r22,where M o is the seismic moment14.The spatial derivative
of
04080120160200
Time (ms)
Figure2|An example of the raw seismograms obtained from a horizontal component in the two periods.Both are filtered with a band pass filter of 1–5kHz.Inset shows the first30ms of the waveforms.P and S indicate the compressional and shear wave arrivals.
NATURE|Vol454|10July2008LETTERS
205
displacement,strain,thus should be /M o r 23.The static stress
change at r is Ds ~a m L 2D r 3~a m (D =L )(r =L )3
~a Ds o
^r 3,where D s o is the
average static stress change along the fault,^r is the characteristic distance measured in fault lengths (L ),D is slip on the fault,m is the shear modulus,and a is a scaling constant equal to 1/(6p )(ref.14).If w
e assume a static stress change in the range of 3to 10MPa (refs 15,16),then the static coseismic stress change at the SAFOD site is estimated to be ,250–833Pa for the M 53event,which is a few times lower than the total stress change (2.3kPa)calculated from the amount of delay time during the first excursion.The predicted static stress changes at the SAFOD site calculated from the entire local seismicity catalogue are shown in Fig.4b.Here we used all the events that occurred within 10km of the site and made a time series of the coseismic stress changes.The M 53earthquake obviously has the
largest effect at the experiment site.The second largest peak around day 20corresponds to a relatively deep event (22November 2005,03:38:02.13UT ,36.01002120.5692,depth 5.07km,M 52.6),which is not observed in the delay time data.The third peak corresponds to the M 51event.It is not clear to us why the larger M 52.6event is not observed while the smaller M 51event shows clearly in the delay time data.But we noticed that data collected in the second period had a better SNR than those of the first period.The associated stress change of the M 52.6event thus might be below the resolution of the first-month data.
Coseismic change was also observed in other geodetic data.We found a step-function change from the borehole fibre-optic strain-meter data at SAFOD (Fig.4a inset)as well as from the surface creepmeter data at Middle Mountain (Fig.4b).The static strain change observed at SAFOD is ,(20–25)3
1029,corresponding to
a
b
60
56
5248
4421239.59.49.302402Time elapsed since 2 Nov. 2005 (days)
0.1
0.21.00.0–2.0
9.3
9.2
P r e s s u r e (104 P a )
2.0
0.0S  d e l a y t i m e  (µs )
S  d e l a y t i m e  (µs )S +c o d a d e l a y t i m e  (µs )
S +c o d a d e l a y t i m e  (µs )21
23T e m p e r a t u r e (ºC )
P r e c i p i t a t i o n (i n c h e s )P r e s s u r e (104 P a )T e m p e r a t u r e (ºC )P r e c i p i t a t i o n (i n c h e s )
0.060.040
5
10
15
20
25
0.00
0.02–1.01.00.0T e m p . c o r r .d e l a y  t i m e  (µs )Figure 3|Estimated delay times for the two periods.a ,First period;
b ,second period.In each panel,the top two traces are delay times estimated from time windows that contain respectively the S-wave arrival and the S-wave arrival plus the coda;the third trace is barometri
c pressure;the fourth trace is temperature;an
d th
e bottom trace is precipitation.In a ,the second to last trace is a temperature corrected version o
f the second trace.Elapsed time is calculated from 2November 2005,00:00:00UT .
LETTERS
NATURE |Vol 454|10July 2008
206
a coseismic stress change of ,600–750Pa,which is of the same order of magnitude as our estimate.On the other hand,there were no obvious changes in the SAFOD GPS,or the FLT and DLT strainmeter records (Fig 4b).The lack of an observable coseismic signal at these sites is,however,predicted by the theoretical amplitude.
The coseismic offset recorded by the SAFOD strainmeter is not obviously present in the delay time data measured either from the
manually stacked 45-min-per-sample data or from the delay times calculated from the 27-s-per-sample raw data.The derivative of the delay time series (dotted line in Fig.4c),however,does reveal that the largest offset of the entire two-month observing period occurred ,30s after the M 53earthquake.This suggests that there was a small coseismic change in the delay time data.The lack of a stronger coseis-mic signal in the delay time data may imply that the velocity changes we observed here are mainly the result of a poroelastic 17rather than an elastic response to abrupt stress changes.
The two travel-time excursions appear to possess significant pre-seismic components.The first excursion was observed to start at 23:34UT on 23December 2005,while the M 53earthquake occurred at ,10.6h later,at 10:10UT on 24December 2005(Fig.4c).The excursion reached a maximum right after the earthquake,peaking at 21:21UT on 24December 2005.The excursion thus has a clear preseismic component besides the coseismic/postseismic changes.The preseismic and coseismic/postseismic components account for ,46%and ,54%of the total change,respectively.This is also true for the second excursion.Its onset is around 22:59UT on 28December 2005,about 2.5h before the occurrence of the M 51earthquake (01:32UT on 29December,Fig.4c).
With the available geodetic instrumentation,it was impossible to further evaluate the preseismic component.The most direct test would have been with the SAFOD borehole strainmeter data.Unfortunately,the low frequency component is severely contami-nated by surface temperature variations and is unusable for periods longer than a few minutes,and is thus not useful in confirming the two low-frequency excursions (M.Zumberge,personal communication).All other instrumentation is either too far away or not sufficiently sensitive to observe even the coseismic offset.Historically,there has been an absence of preseismic signals in geodetic observations,such as a borehole strainmeter.We suggest that this may be the result of two differences between such instruments and o
ur ‘stress meter’.First,our basic measurement is not strain,but rather a stress-induced change in the effective elastic constants of a poroelastic medium,mediated by variations in crack properties and fluid flow.These changes may regis-ter only weakly on a strainmeter,a GPS,or a creepmeter.Second,a conventional strainmeter measures local change in the volume imme-diately surrounding the instrument while our measurements reflect stress/strain changes occurring over a volume sampled by the coda waves that could be orders of magnitude larger.
We put forward the hypothesis that there is a change in effective elastic moduli before rupture,such as a sudden increase in microcrack density,which is a phenomenon related to dilatancy and observed in many laboratory studies 1,2.As such,further continuous seismic mon-itoring might provide an effective tool for understanding the stress changes that accompany and perhaps precede seismic activity.METHODS SUMMARY
We used a specially built piezoelectric source and a Geode recorder (Geometrics)to generate and record seismic waves travelling along a ,10m baseline near the San Andreas fault at ,1km depth.A cosine fitting method was used to estimate the S-wave travel time to subsample precision.
Full Methods and any associated references are available in the online version of the paper at www.nature/nature.
Received 22October 2007;accepted 12May 2008.
1.Brace,W.F.,Paulding,B.W.&Scholz,C.H.Dilatancy in the fracture of crystalline rocks.J.Geophys.Res.71,3939–3953(1966).
2.Scholz,C.H.Microfracturing and the inelastic deformation of rock I:Compression.J.Geophys.Res.73,1417–1432(1968).
3.Birch,F.The velocity of compressional waves in rocks to 10kilobars,part 1.J.Geophys.Res.65,1083–1102(1960).
4.Birch,F.The velocity of compressional waves in rocks to 10kilobars,part 2.J.Geophys.Res.66,2199–2224(1961).
5.Nur,A.&Simmons,G.The effect of saturation on velocity in low porosity rocks.Earth Planet.Sci.Lett.7,183–193(1969).
6.
Walsh,J.B.The effect of cracks on the compressibility of rock.J.Geophys.Res.70,381–389
(1965).
b
5
10
–10
–5
5
10
D e p t h  (k m )
NW
SE
a
X M M c r e e p (m m )S A F O D G P S (m m )
Along fault distance (km)
20×10–9
10 s time
s t r a i n
SAFOD strain
c
S t r e s s c h a n g e (k P a )
S +c o d a d e l a y t i m e  (µs )E s t i m a t e d d e l a y t i m e  (µs )
Elapsed time since 2 Nov. 2005 (days)
Predicted stress change (KPa)
0.00.20.4Figure 4|A comparison of delay time variations with local seismicity and other deformation measurements.a ,Depth distribution of earthquakes that occurred in the experimental period.Red square,the SAFOD
experiment site;red and green circles,the M 53and M 51earthquake,respectively.Inset,the SAFOD strainmeter record,which shows a step-function coseismic strain change.The low frequency content of the
strainmeter data is severely contaminated by surface temperature variations,and is consequently not suitable for analysis.b ,Top to bottom:creep measurement at Middle Mountain (XMM);GPS measurement of fault-parallel motion at the SAFOD site;calculated static coseismic stress changes at the SAFOD experiment site for all of the earthquakes;and delay times estimated from the S wave plus its coda for comparison.Dashed lines indicate the time when the M 53and M 51earthquakes occurred.Note that the amplitude of the stress change of the M 53event (,0.5kPa)is saturated in this plot.c ,Predicted coseismic stress changes at SAFOD for earthquakes occurring between 22December 2005(day 50)and 1January 2006(day 60)indicated by shading in b are shown with the delay time estimation.Stress changes from the local seismicity between days 55and 60are amplified by a factor of 10.The two filled circles show the stress change of the M 53(red)and M 51(green)event,respectively.The vertical lines indicate the occurrence times of the M 53and the M 51event,and the red and green (upward)arrows show the onset times of the two excursions.Blue dotted line is the derivative of the delay time series.Notice that the largest change occurred about ,30s after the M 53earthquake.
NATURE |Vol 454|10July 2008LETTERS
207
7.Nur,A.Effects of stress on velocity anisotropy in rocks with cracks.J.Geophys.
Res.76,2022–2034(1971).
8.O’Connell,R.J.&Budiansky,B.Seismic velocities in dry and saturated cracked
solids.J.Geophys.Res.79,5412–5426(1974).
9.De Fazio,T.L.,Aki,K.&Alba,J.Solid earth tide and observed change in the in situ
seismic velocity.J.Geophys.Res.78,1319–1322(1973).
10.Reasenberg,P.&Aki,K.A precise,continuous measurement of seismic velocity
for monitoring in situ stress.J.Geophys.Res.79,399–406(1974).
11.Leary,P.C.,Malin,P.E.,Phinny,R.A.,Brocher,T.&Voncolln,R.Systematic
monitoring of millisecond traveltime variations near Palmdale,California.J.
Geophys.Res.84,659–666(1979).
12. al.Long-term observation of in situ seismic velocity and
attenuation.J.Geophys.Res.108,doi:10.1029/2002JB002005(2003).
13.Silver,P.G.,Daley,T.M.,Niu,F.&Majer,E.L.Active source monitoring of cross-active下载
well seismic traveltime for stress-induced changes.Bull.Seismol.Soc.Am.97, 281–293(2007).
14.Aki,K.&Richards,P.G.Quantitative Seismology(Freeman,New York,1980).
15.Abercrombie,R.E.Earthquake source scaling relationships from21to5M L using
seismograms recorded at2.5-km depth.J.Geophys.Res.100,24015–24036(1995).16.Rubin,A.M.&Gillard,D.Aftershock asymmetry/rupture directivity among
central San Andreas fault microearthquakes.J.Geophys.Res.105,19095–19109 (2000).
17.Segall,P.,Jonsson,S.&Agustsson,K.When is the strain in the meter the same as
the strain in the rock?Geophys.Res.Lett.30,doi:10.1029/2003GL017995(2003).
18.Unsworth,M.,Bedrosian,P.,Eisel,M.,Egbert,G.&Siripunvaraporn,W.Along
strike variations in the electrical structure of the San Andreas Fault at Parkfield, California.Geophys.Res.Lett.27,3021–3024(2000).
Acknowledgements We thank the NSF funded SAFOD programme and those involved in providing the experiment site,R.Trautz for supplying the barometric pressure logger,M.Zumberge for providing the SAFOD strainmeter data,and D.Lippert and R.Haught for helping with field work.This work was supported by the NSF,Rice University,the Carnegie Institution of Washington,and Lawrence Berkeley National Laboratory of the US Department of Energy under contract DE-AC02-05CH11231.
Author Information Reprints and permissions information is available at www.nature/reprints.Correspondence and requests for materials should be addressed to F.N.(niu@rice.edu).
LETTERS NATURE|Vol454|10July2008 208
METHODS
Data acquisition system.Our acquisition was conducted with a combination of commercial and specially built equipment.The latter are the piezoelectric source and the high voltage amplifier used to
power it.The source includes18cylin-drical rings of piezoelectric ceramic(lead zirconate titanate)epoxied together and wired for positive and negative voltage on the inner and outer surfaces.The source was fluid coupled to the well casing.A three-component accelerometer was clamped to the well casing to provide coupling and reduce relative motions between the source and receiver.We used a commercial recording system,a ‘Geode’manufactured by Geometrics,which has a24bit analogue-to-digital converter.An air conditioner and heater were used to maintain the recording system electronics within a temperature range of about61u C.
Triggering was used in our data recording system.The digitizer continually samples the data,and receives a trigger that will generally be between two digitized samples.Including a section of pre-trigger data,the time series is interpolated and re-sampled,so that the new time series begins at the time of the trigger.This start time is not exact,and,at a sampling rate of48,000s21,this time is computed to the nearest twentieth of a sample(Geometrics engineering,personal communication).Thus there is a delay time measurement error that will be at most a fortieth of a sample (half-way between samples),and the average error will be an eightieth of a sample, assuming that the errors are uniformly distributed.This corresponds to an average error of260ns per trigger.The error in the stacked data decreases by a factor of N1/2, assuming the errors are uncorrelated.For N5100,we obtain a timing error of26ns. Optim
um experimental design.As shown by Silver et al.13,there is an optimum distance between the source and receiver that minimizes the detectable threshold of subtle velocity changes:
N~Q=pð1ÞHere N is the number of wavelengths between the source and receiver and Q is the quality factor.At the SAFOD site,Q is around200,which gives N564.If we assume the S-wave velocity to be2.8m ms21,then the wavelength of the signal with a dominant frequency of2kHz is about1.4m,so the optimum distance is ,90m.As it was necessary to perform the experiment in the available boreholes, our cross-hole distance was limited to10m,which while not optimal still pro-vided us with a good SNR.
Subsample delay time estimate(DTE).In this study,we employed a cosine fitting method to estimate subsample delay time in the time domain19,20.Given the largest sample of the correlation function,cc(0),and its two neighbours cc(21)and cc(1),the estimated subsample shift is given by following expression:
t~a=arctan cc({1){cc(1)
ð2Þwhere:
a~arccos
cc({1)z cc(1)
2cc(0)
ð3Þ
Error estimation.Silver et al.13derived a low bound of the error in delay time measurements:
s DTE§
1
2p f0SNR
ð4Þ
Here f0is the dominant frequency of the source pulse,and SNR is the signal-to-noise ratio.Equation(4)indicates that SNR is the only parameter that controls the precision in our delay time estimation when the digitizing error is much less than the background noise in this regime.The precision is not controlled by the sampling rate of the digitizer so it is possible to obtain subsample-inte
rval mea-surements of the delay time.The dominant frequency of our data is2kHz and the SNR is around700for the45-min stacked data,resulting in a best achievable precision of,1.131027s,or110ns in the DTE.
We also measured delay time between each two consecutive samples,which follows a gaussian distribution with a standard deviation of,80ns and,50ns for the first and second recording period,respectively.In general they are com-parable to or even better than the theoretical low bound in equation(4).Since there is contribution from the actual stress-induced velocity perturbations in the measurement,our actual precision can be better than the measured standard deviations.Thus the lower bound appears to be larger than the true DTE error. One possible explanation is that the SNR is significantly underestimated,as the noise is estimated from a time window before the first arrival,which actually contains a considerable amount of non-random electronic noise known as cross-talk,and non-random‘wrap-around’noise from the previous shot.
The precision discussed here does not include other systematic non-random noise,such as changes in the source pulses,errors in trigger timing and digitizer’s clock.Such systematic errors could lead to a long-term trend in DTE.To estimate these effects,we also recorded the source pulse waveform in addition to the data. We employed the same method to measure the variation in the source pulse width.
Changes in the source pulse width are between620ns.This indicates that our source pulse generator and recording system were very stable in the two periods and timing error in the digitizer clock was also very small.
19.Cespedes,I.,Huang,Y.,Ophir,J.&Spratt,S.Methods for estimation of sub-
sample time delays of digitized echo signals.Ultrason.Imaging17,142–171(1995).
20.De Jong,P.G.M.,Arts,T.,Hoeks,A.P.G.&Reneman,R.S.Determination of tissue
motion velocity by correlation interpolation of pulsed ultrasonic echo signals.
Ultrason.Imaging12,84–98(1990).
doi:10.1038/nature07111

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。