Shear-Wave Splitting Observed in Local Earthquake Data on the Reykjanes Peninsula, SW Iceland

Examensarbete vid Institutionen för geovetenskaper ISSN 1650-6553 Nr 273

Shear-Wave Splitting Observed in Local Earthquake Data on the Reykjanes Peninsula, SW Iceland

Shear-Wave Splitting Observed in Local Earthquake Data on the Reykjanes Peninsula, SW Iceland

Darina Buhcheva

Darina Buhcheva

Uppsala universitet, Institutionen för geovetenskaper Examensarbete E, 30 hp i Geofysik

ISSN 1650-6553 Nr 273

Tryckt hos Institutionen för geovetenskaper, Geotryckeriet, Uppsala universitet, Uppsala, 2014.

Shear-wave splitting is a phenomenon observed in almost all in situ rocks. Due to propagation through stress-aligned and fluid-saturated microcracks and fractures the initial shear wave splits into two almost orthogonal waves which propagate with different velocities along similar ray paths. The process is characterized by the polarization direction of the faster split shear wave, which is parallel to the orientation of the cracks, and the time delay between the onsets of the two waves. The analysis of shear-wave splitting has been conducted over records of 233 microearthquakes in the vicinity of five seismic stations in SW Iceland. Visual methods have been applied to the data to retrieve the final results for polarization directions and time delays.

The main polarization azimuth for the leading split wave is N30°- 60°E which is in full agreement with the mapped alignments of normal faults and volcanic fissures in the surface. The time delays measured at different sites vary in the range of 10-100 ms for the events of best quality. In general, splitting times do not show a clear pattern at all recording sites with increasing depth. The only firm conclusion that can be drawn from the time delays is that at station BLF in the Brennisteinsfjöll fissure swarm, the time delays are smaller than in the Hengill area and therefore the strength of anisotropy beneath that station appears to be lower.

Examensarbete vid Institutionen för geovetenskaper ISSN 1650-6553 Nr 273

Shear-Wave Splitting Observed in Local Earthquake Data on the Reykjanes Peninsula, SW Iceland

Darina Buhcheva

Copyright © %BSJOB
#VIDIFWB and the Department of Earth Sciences Uppsala University Published at Department of Earth Sciences, Geotryckeriet Uppsala University, Uppsala, 201

Abstract

Shear-wave splitting is a phenomenon observed in almost all in situ rocks.

Due to propagation through stress-aligned and fluid-saturated microcracks and fractures the initial shear wave splits into two almost orthogonal waves which propagate with different velocities along similar ray paths. The process is characterized by the polarization direction of the faster split shear wave, which is parallel to the orientation of the cracks, and the time delay between the onsets of the two waves. The analysis of shear-wave splitting has been conducted over records of 233 microearthquakes in the vicinity of five seismic stations in SW Iceland. Visual methods have been applied to the data to retrieve the final results for polarization directions and time delays. The main polarization azimuth for the leading split wave is N30°- 60°E which is in full agreement with the mapped alignments of normal faults and volcanic fissures in the surface. The time delays measured at different sites vary in the range of 10-100 ms for the events of best quality. In general, splitting times do not show a clear pattern at all recording sites with increasing depth. The only firm conclusion that can be drawn from the time delays is that at station BLF in the Brennisteinsfj¨oll fissure swarm, the time delays are smaller than in the Hengill area and therefore the strength of anisotropy beneath that station appears to be lower.

i

List of Figures

1 Schematic description of shear-wave splitting . . . 1

2 Seismic activity on Reykjanes Peninsula for the period August 2010 - August 2012. . . 4

3 Map of the topography, fissure swarms and volcanoes on the Reyk- janes Peninsula, together with the local seismicity around the record- ing sites. . . 6

4 1-D P-wave velocity model for the area of Hengill volcano. . . 10

5 Coordinate system used in the relocation calculations. NE coordi- nates are transformed into RT. α is the azimuth and θ is the back azimuth. . . 11

6 Typical example of a local event (with θ = 47.7°, i = 11.5°), recorded within the shear-wave window of station HVD. . . 12

7 Simplified hypocentre-recording geometry. . . 13

8 Old and new epicentres of the earthquakes. . . 15

9 Old and new hypocentres of all events used. . . 16

10 Typical examples of events from stations KRI and BJA. . . 19

11 Determination of the polarization direction for two events at BJA and BLF. . . 22

12 Comparison between the initial and second approach for retrieving the polarizations. . . 25

13 Rose diagrams of the fast direction from earthquakes recorded at KRI. 25 14 Map of the Reykjanes Peninsula, showing polarization rose diagrams for the five stations and from four additional stations in the area. . . 26

15 Variation of time delay with focal depth. . . 28

16 Variation of time delay with focal depth for the stations in the Hengill area. . . 29

17 Normalized time-delay distribution. . . 31

18 Geometry of the incidence and resulting rays. . . 38

19 Synthetic data for one-lobed waveform without added noise. . . 43

20 Synthetic data for one-lobed waveform with added noise (SNR ≈ 3.71). . . 44

21 Synthetic data for one-lobed waveform with added noise (SNR ≈ 1.65). . . 45

22 Synthetic data for two-lobed waveform without added noise. . . 46

ii

23 Synthetic data for two-lobed waveform with added noise (SNR ≈ 3.75). . . 47 24 Synthetic data for two-lobed waveform with added noise (SNR ≈

1.60). . . 48 25 Synthetic data for three-lobed waveform without added noise. . . 49 26 Synthetic data for three-lobed waveform with added noise (SNR ≈

3.78). . . 50 27 Synthetic data for three-lobed waveform with added noise (SNR ≈

1.66). . . 51 28 Synthetic data for four-lobed waveform without added noise. . . 52 29 Synthetic data for four-lobed waveform with added noise (SNR ≈

3.77). . . 53 30 Synthetic data for four-lobed waveform with added noise (SNR ≈

1.73). . . 54

iii

Contents

Abstract i

List of figures ii

1 Introduction 1

2 Area of interest 5

3 Data processing methodology 7

4 Event location 10

5 Techniques for measuring of shear-wave splitting 17

5.1 Visual measurements . . . 17 5.2 Difficulties in measuring and scattering effects . . . 23

6 Results 24

6.1 Strike of the fractures . . . 24 6.2 Magnitude of anisotropy . . . 27

7 Discussion and conclusions 32

Acknowledgements 35

References 36

Appendices 38

A Free surface correction 38

B Synthetic data 42

v

1 Introduction

Polarizable waves, such as shear waves, split when passing through an anisotropic medium in the Earth’s crust. The splitting of a near vertically propagating shear wave can be char- acterized using four parameters in the ideal case as follows: orientation of the polarizations of the fast and the slow waves, time delay between the two arrivals and their relative at- tenuation. However, it is possible to estimate only the time delay and the polarization of the first arrival from the recordings of earthquake signals through visual or automatic measurements (Chen et al., 1987). The polarization of the fast shear wave is parallel to the direction of maximum horizontal compression or the local strike of the cracks in the area (Crampin and Lovell, 1991), while the time delay is considered to be related to the magnitude of anisotropy per unit volume along the ray path and the distance passed through the anisotropic medium (figure 1).

Figure 1: Schematic description of shear-wave splitting. Arbitrarily polarized incident shear wave S splits into two waves, S₁ and S₂, which have polarizations parallel and perpendicular to the stress-aligned and fluid-saturated microcracks and fractures in the anisotropic medium (Rial et al., 2005).

Shear-wave splitting due to anisotropy has specific particle-motion patterns that can be revealed through analysing polarization diagrams or rotating the horizontal seismo- grams parallel and orthogonal to preferred polarizations (Volti and Crampin, 2003). The

phenomenon has been observed in sedimentary, igneous and metamorphic rocks below a critical depth around 500-1000 m (Crampin, 1994). Closer to the surface, the splitting could be perturbed by predominant near-surface effects (Volti and Crampin, 2003).

The anisotropy may also be caused by constantly present fractures if they are aligned by the long-term stress field (normal faulting in SW Iceland due to spreading between the tectonic plates). However, the in situ current stress field can modify the anisotropic effect by changing the aspect ratio of fractures and microcracks. The distributions of stress-aligned cracks are known as extensive-dilatancy anisotropy (EDA) (Crampin et al., 1984). The fluid-filled EDA-cracks are aligned along the pressure gradients through fluid migration (Rutter and Brodie, 1991) in a way that the rock turns effectively anisotropic to propagating waves (Chen et al., 1987). The geometry of the cracks is highly compliant to alterations in stress direction or magnitude and increase in the stress results in the alignment of the EDA cracks. Zatsepin and Crampin (1997) and Crampin and Zatsepin (1997) suggest an anisotropic poroelasticity (APE) model which describes the dynamics of fluid-saturated EDA cracks under influence of changing stress field, taking also into account variations of the pore-fluid pressure and its effect on the cracks. APE-model simulations show that an increase of the aspect ratios of the microcracks due to an increase of the stress and an immediate decrease of the initial values after stress release results in changes of the shear-wave splitting, similar to some variations before and after an earthquake. Thus, temporal changes in the intensity of anisotropy are claimed to be a tool for monitoring the earthquake-stress cycle.

Structural geological features such as aligned dikes and aligned heterogeneity (e.g.

volcanic fissures and hyaloclastite ridges) may also contribute to the anisotropy if the wavelengths are larger than the scale of these features. Azimuthal anisotropy, a directional dependence of shear-wave velocity with azimuth, could be observed in the presence of vertically aligned geological elements or formations. Radial anisotropy, a dependence of the velocity on the deviation of the polarization direction from the vertical with the length of the ray path, could be caused by horizontal layers and inclusions.

Usually, sets of small local earthquakes having incidence angles within the shear-wave

window are used to study shear-wave splitting in the upper crust. The shear-wave window is a vertical cone under a seismic receiver with aperture defined by sin⁻¹(VP/VS) = 35.26°

for a Poisson’s ratio equal to 0.25 (Booth and Crampin, 1985). The span of the shear- wave window guarantees that possible splitting would not be obscured by effects related to the interaction of the waves with the free surface and evanescent compressional-wave conversions would not occur due to the small incident angles (Menke et al., 1994). Since the ray paths bend while propagating due to the relatively steep vertical velocity gradient, this window could be expanded to approximately 45° (Volti and Crampin, 2003). It should be noted as well that the orthogonality of the fast and slow shear waves is valid only in case of normal incidence. Otherwise, this characteristic is distorted by near-surface effects (Booth and Crampin, 1985).

Over the past three years Uppsala University (UU) has in collaboration with Mas- sachusetts Institute of Technology (MIT), Reykjavik University and other Icelandic insti- tutions operated a dense network of about 50 seismographs in the area between Hengill volcano and Reykjanes in Southwestern Iceland (figure 2). In this project we apply mea- surements of shear-wave splitting of near-vertically travelling shear waves at a number of these stations. It provides an opportunity to look at this significant property of the upper crust in much more detail than has hitherto been achievable and therefore possibly re- solve differences between the fissure swarms of the rift zone and adjacent areas. Since the network of seismographs is much denser than the networks available for previous studies, earthquake locations are better constrained, including their depth. The idea of this study is to see if it is possible to confine the depth distribution of anisotropy in the upper crust - a very interesting prospect, since it could significantly constrain the crustal accretion process and can be taken as a proxy for fracture permeability in the crust. The depth extent of fracture permeability is a very important parameter in the context of plans to drill for supercritical fluids at depth (3 - 5 km) for geothermal exploitation.

Figure 2: Seismic activity on Reykjanes Peninsula for the period August 2010 - August 2012. The earthquakes are marked with red circles. The triangles present the recording sites (the green ones are temporary placed for that time, the yellow ones are the permanent stations of the Icelandic Meteorological Office (IMO) network). The earthquakes shown here have been detected and located by the IMO network only. Refinements of the seismicity based on the denser project network are still being worked on.

2 Area of interest

The Reykjanes Peninsula in Southwestern Iceland sits astride the Mid-Atlantic ridge rift system. The upper crust is heavily fractured by the rifting process (normal faults).

In addition, volcanic fissures in the area are parallel to the normal faulting and fissure eruptions in the area have during the Quaternary produced confined parallel hyaloclastite ridges. This geological structure is expected to produce strong azimuthal anisotropy in seismic properties. Volcanism in the area during the warm periods of the Quaternary have on the other hand produced large subhorizontal lava flows, which due to their internal structure related to dry or whet cooling are likely to cause radial anisotropy. Earlier studies in the area that are based on near-vertically travelling shear waves from local earthquakes suggest that a significant component of azimuthal anisotropy exists in the upper crust (Volti and Crampin, 2003, Menke et al., 1994). Waves of this geometry are not sensitive to the potential radial component of anisotropy.

This project is mainly focused on the Hengill volcanic complex, which is located where three tectonic systems intersect: the Reykjanes Peninsula, the West Volcanic Zone (WVZ) and the South Iceland Seismic Zone (SISZ) - two spreading zones between the North American and European plates to the north and south and a transform zone. The area of the Hengill volcano is crossed by a fissure swarm which extends from the south coast of Iceland near Selvogur to Lake Thingvallavatn in the north (Jousset et al., 2011).

The stations Bjarnastaðir (BJA), Blafjöll (BLF), Hveradalir (HVD), Litli-Hals (LHA), which have been the main subject of this project, are placed at the four ends of an area in the southern part of the Hengill fissure swarm. Events from station Krýsuvík (KRI), located in the Reykjanes Volcanic Zone, were also taken into account in order to compare results with other research and confirm the reliability of one of the applied approaches.

The station coordinates are BJA: 63.94590N, -21.30258E; BLF: 63.97368N, -21.64893E;

HVD: 64.02451N, -21.39947E; LHA: 64.02405N, -21.04887E; KRI: 63.87810N, -22.07646E.

The five stations together with the fissure swarms, volcanoes and local seismicity in radius of 5 km around the stations in the period August 2010 - August 2011 are shown on figure 3.

BLF is located in the Brennisteinsfjöll fissure swarm, parallel to Hengill’s to the west.

HVD lies in the middle of the Hengill fissure swarm. The location of LHA is near the eastern end of the area, outside the fissure swarm. The seismicity around LHA is caused by an aftershock series after two strong earthquakes with M∼6.0 in 2008 along the two parallel N-S strike-slip faults laying south and south-west of the station. BJA is in the south and is located in the SISZ and near the eastern edge of the Hengill fissure swarm where anomalous seismicity has been observed since recordings started (around 1990).

Figure 3: Map of the topography, fissure swarms and volcanoes (elongated closed con- tours) on the Reykjanes Peninsula, together with the local seismicity of interest. The displayed earthquakes (in red) are in radius of 5 km around each station.

In 1960s, the potential for harnessing natural geothermal energy has been realized in this triple-junction complex. Hot springs and fumaroles are the surface manifestations of the vigorous geothermal activity in the area (Tryggvason et al., 2002). The extracted hot fluids are used for heat production and electrical power generation. Four geothermal power plants are built on the Reykjanes Peninsula, two of which are in the area of the Hengill volcanic complex. The total electrical power generation capacity is about 550 MW and the thermal capacity is comparable.

3 Data processing methodology

A few steps have been followed in order to extract the final earthquake data set used in the project. Each event should fit within all the criteria described below.

• 30°- window: The typical 45°-shear-wave window that should limit the incidence an- gles of the incoming waves is restricted to 30°-window inside of which the hypocen- tres should be located. This requirement leaves just events with incidence angles closer to the vertical so the earthquake records are less contaminated by P-wave conversions. At this stage, catalogue locations from IMO are used.

• M = 0.5: This threshold of the magnitude is chosen with the idea to remove data with unclear P- and S-wave arrivals assuming that the wave onsets will be more visible above the background of the regular noise.

• Signal-to-noise ratio and scattering estimation: The signal-to-noise ratio (SNR) es- timation is applied in order to verify the reliability of the signals. The criterion is used to check the ratio between the amplitudes of the noise (P-wave coda and back- ground noise) before the following S-wave signal based on the greatest amplitude of the waveform. The threshold for the SNR is set to 2.0 for the north and east compo- nents. The maximum amplitude of the S-wave window on the vertical component is compared with the maximum amplitudes on the horizontals. If this ratio is greater than 0.5¹, then the event is discarded as this is an indicator of strong scattering.

A reasonable explanation for the high amplitudes on the vertical component during the S-wave arrival is strong scattering and conversion along the ray paths.

• Location in the shear-wave window: After relocation (as explained in the following section), we check if each earthquake still is located within the shear-wave window.

These criteria describe the general approach to the data and how the final data set was obtained. However, not all stations have been subjected to exactly the same procedure as

1For BJA (because of the importance of the station) this value is raised to 0.67 due to the small number of events which managed to pass the limit.

you can see from the steps included in table 1 which shows the number of events during application of the criteria to the data set. For three of the five stations there are some differences in the finalizing of the events. Since LHA was the first station processed, after step 1 and 2, we were looking at the particle motion and chose those events which have elliptical motion, typical when shear-wave splitting is present. The events related to BJA showed high amplitudes after the S-wave arrival on the vertical component and most of them were supposed to be discarded because they did not pass the 0.5-threshold. BJA has a very important role in a shear-wave splitting study and if one is testing a method, since local earthquakes in its vicinity have been used in previous studies in this field (Volti and Crampin, 2003), giving good results. Therefore, the threshold was increased to 0.67.

After working with HVD, LHA and BJA, BLF had a strategic location for closing the area. Searching for earthquakes around that station when applying steps 1 and 2 on the earthquakes in the catalogue resulted in a small and poor data set. Thus, instead of events with bigger magnitude than 0.5 and located in the shear-wave window, events which have occurred within a radius of 5 km epicentral distance have been considered.

Since the hypocentres are subjected to relocation and qualitative processing, the final set would contain only the events that fit the criterion to be within the shear-wave window after relocation.

During the recording period, there are around 20 or less earthquakes for stations LHA, BJA and BLF and more than 80 for KRI and HVD, located within the 30°- window, that fit the criteria for analysis.

StepCriterionStation LHABLFBJAKRIHVD 1M>0.5and30°-window°(fromthecatalogue)6465379179 Inradiusof5kmaroundthestation114 2Availabledigitalrecordsforanearthquakeatthecloseststation469064224166 3ClearP-andS-arrivals,goodqualityofthedata,etc.464721217147 Ellipticalparticlemotion22 4HighSNR,nostrongscattering,inthechosenshear-wavewindow2115219784 5Extraconsiderationofthewaveformsorotherfactors1282 Final:2112219782 Table1:Numberofeventsforeachstationleftafterapplyingthedifferentselectioncriteria.

4 Event location

For the initial processing of the earthquakes (as explained in more detail in the previ- ous section) catalogue locations have been used. Those are calculated with measured arrival times at the IMO seismic stations. The permanent stations forming the Icelandic Seismological Network are rather sparse in the area of Hengill volcano. Therefore, the catalogue locations are considered not as accurate as needed for the purpose of a research in shear-wave splitting. In particular, the depth determinations are subject to signifi- cant uncertainty. Commonly, the estimated location uncertainty is in the range 0.5 to 1 km in the horizontal directions and 1 to 3 km in depth in between the stations of the network. These uncertainty estimates do not include potential location bias due to the use of a potentially suboptimal one-dimensional velocity model or due to ignored lateral heterogeneity (or anisotropy). In order to refine the location, in particular the depth of each event, a one-dimensional velocity model based on travel time inversion in the region (figure 4) is used. The new hypocentres are calculated based on S-P arrival time differ- ence and estimates of the back azimuth and incidence angle from the observed waveforms, where these three parameters are measured from the recordings at the closest station.

Figure 4: 1-D P-wave velocity model for the area of Hengill volcano based on Tryggvason et al. (2002).

By definition, the azimuth (α) is the angle, measured at the epicentre, from the north in the direction of the station. The back azimuth, respectively, is θ = α + 180°(assuming linear propagation of the waves (no lateral refraction) and negligible effects of the spherical Earth) and is measured at the station towards the epicentre, as shown on figure 5. θ is estimated as the angle where the transverse component has minimal energy, after the horizontal components are rotated into radial (R) and transverse (T) directions. This rotation is done using a transformation matrix:





 R T





=







cos α sin α

− sin α cos α











 N E





, (1)

so the new components in terms of the azimuth are:

R = cos αN + sin αE T = cos αE − sin αN.

(2)

Substituting α with θ, the radial and transverse components expressed through the back azimuth are:

R = − cos θN − sin θE T = − cos θE + sin θN.

(3)

Figure 5: Coordinate system used in the relocation calculations. NE coordinates are transformed into RT. α is the azimuth and θ is the back azimuth.

Figure 6: Typical example of a local event (with θ = 47.7°, i = 11.5°), recorded within the shear-wave window of station HVD. The seismogram (a) shows a good SNR on all three components. Red lines are used to mark the P- and S-wave picks of the arrivals.

Zoomed section of the initial P-wave motion (b) is shown before and after the rotation into ZRT coordinate system. The window between the coloured lines contains the samples used to estimate θ and i. The compilation of each component against the others after the rotation (c) is displayed and the initial motion is coloured in red.

We use the onset of the P-wave to estimate the angle of incidence, i. Specifically, we use the time window at the beginning of the P-wave while we still see clear linear body-wave motion, as shown on figure 6c. On the plot of the vertical against the radial component, the slope of the initial motion (in red) or the ratio between the amplitudes of the vertical and radial components is used for determining the incidence angle. A correction for the free surface effects (thoroughly explained in appendix A) has been applied for improving the results. The incidence angle for the S-waves is assumed to be the same as for the P-waves. Having estimated back azimuth, incidence angle and the differential time between S- and P-waves, the event is relocated by tracing a ray in the chosen 1D velocity model from the receiver with the initial direction specified by back azimuth and incidence angle. The ray is terminated when the measured S minus P travel time is reached. The end point of that ray is the new location estimate.

Back azimuth and incidence angle may be difficult to estimate accurately. However, the differential time of S- and P-waves is generally estimated with high precision (typical uncertainty of 0.03 seconds). Given that we are working with events near the particular station being analysed, and that although the incidence angle may be relatively poorly determined from the P-wave particle motion the incidence angles are generally small, this relocation method is likely to reduce the uncertainty of the depth determination significantly.

Figure 7: Simplified hypocentre-recording geometry.

Consider the simplified geometry in figure 7: an earthquake occurs at depth z, with distance s between the station and the hypocentre, an incidence angle i for the incoming waves, a time difference t between the arrivals of the S- and P-waves at the station, and wave velocities V_p and V_s. In this case, t can be expressed with the equation:

t = z cos i

1

V_p (γ − 1) , (4)

where γ = Vs

V_p is the Poisson’s ratio. Therefore, the uncertainty of the time δt expressed through the uncertainties of the parameters z, Vp, i and γ is

δt = δz cos i

1

V_p (γ − 1) − z cos i

δV_p

V_p² (γ − 1) + z cos i

sin i cos iδi 1

V_p (γ − 1) + z cos i

1

V_pδγ. (5)

Substituting with t from (4) in (5) leads to the equation

δt t = δz

z − δV_p Vp

+ δi tan i + δγ

γ − 1. (6)

Isolating δz/z gives

δz z = δt

t +δVp

V_p − δi tan i − δγ

γ − 1, (7)

and the variance of the depth can be calculated as:

σ_z z =

s

σ_t t

2

+ σVp

V_p

2

+ (σ_itan i)²+

 σ_γ γ − 1

2

, (8)

assuming that all sources of error are independent. Taking some average, but reasonable, values for the parameters and their uncertainties: V_p ≈ 4.5 km/s with δV_p ≈ 0.3 km/s, i ≈ 10° with δi ≈ 5°, t ≈ 0.5 s with δt ≈ 0.03 s and γ ≈ 1.78 with δγ ≈ 0.02, we estimate the average depth uncertainty to be about 10%, i.e. δz ≈ 300m for an event at 3 km depth. This is an improvement compared to the 1 - 3 km uncertainties reported in the IMO catalogue, which are potential underestimates. The event depths are the location parameter of primary interest here, both for estimating the strength of anisotropy and for selecting the events inside the shear-wave window.

Figures 8 and 9 present comparisons between the old and new locations for epicen- tres and hypocentres, respectively. The new locations are more concentrated around the station they are closest to, more clear for the stations with higher number of events like HVD and KRI. It should be noted that the shown new locations are uncertain, in par- ticular because the incidence angle and back azimuth are uncertain. They are designed to improve the depth estimate. By looking at the hypocentres of events around station HVD before and after the data have been refined, one could see that the average depth of the earthquakes is shallower when constraints on the location from the nearest station are available. For the other stations there are not such severe shifts of the hypocentres.

Figure 8: Epicentres of the earthquakes. The catalogue locations are shown in red and the new estimates in yellow.

Figure 9: Hypocentres of all events used. Each of the five graphs shows the catalogue (empty circles) and the relocated (filled circles) hypocentres of the earthquakes for the indicated station. The red line marks a 30°-window.

5 Techniques for measuring of shear-wave splitting

In this study, only visual measurements of the time-delays and polarizations have been used. At the beginning of the research, the automatic measuring technique proposed by Gao et al. (1995) was applied, but without giving any reliable results due to the pronounced complexity of the data. In it the cross-correlation function (CCF) is used to estimate the correlation between two waveforms obtained after rotating the horizontal components over a range of possible angles and applying different time delays. The combination of rotation angle and time delay that gives the highest correlation is an estimate of the best polarization direction and time difference between the two orthogonal components.

From the point of view of measuring the orientation of the fast polarization and time delay to the slower orthogonal polarization of a split shear wave, most of the noise is signal generated. It consists of the tail of P-wave coda before the arrival of the S-waves in addition to scattered S-waves following the direct S-waves and possibly some S-to- P converted energy just before the S-waves. This noise is not strongly affected by the size of the earthquake. Scattering events are evident in the data. Even if the S-waves rise clearly above the P-wave coda, the S-wave waveform seldom consists of just two orthogonal waves. The application of Gao’s automated method is therefore difficult since the apparent waveforms of the two orthogonal waves do not correlate as they are obscured by scattered phases. We therefore refrain to visual inspection of the particle motion of the shear waves. In this section, we present an overview of the visual methods used for estimating polarization directions and time delays, as well as the difficulties faced in working with shear-wave splitting. We also discuss the scattering effects in the crust causing some of those difficulties.

5.1 Visual measurements

In the case of measuring the parameters of shear-wave splitting, visual methods have proven to give better results. The complicated character of the data recorded at the present seismic stations in SW Iceland has given a reason to change the approach in the

visual measurements.

The strategy we started with is quite straightforward and is widely used in such measurements (for example in Volti and Crampin (2003)). We estimate the polarization of the faster split shear wave based on the initial motion on the horizontal polarization diagram from the onset of the shear wave. Then the horizontal components are rotated into components that are parallel and orthogonal to the polarization direction of the faster split shear-wave. The time delay is measured from the seismograms as the difference between the first and second split shear-wave arrivals. According to the quality of the data and how difficult it is to determine the time delay, each event is assigned a quality factor: 1 - good, 2 - fair, 3 - doubtful.

Examples of the stages of the working process are shown in figure 10 for two events from stations KRI and BJA. The figure shows the original three-component seismograms together with the picks of the S-wave, marked on the east components. The particle- motion diagrams of the two horizontal components are displayed and the typical elliptical pattern, a characteristic for present shear-wave splitting, is observed. The motion has a clear initial direction which evolves into a more elliptical shape when the slow wave arrives at the receiver. After rotation the split shear waves can be clearly identified and the time delay between their onsets them can be measured.

As it will be shown later on, this approach gives consistent results for the general polarization directions for KRI, HVD and BJA. Since the data sets for the first two stations are more than four times bigger than for the other three recording sites, and, furthermore, the polarization diagrams are consistent with the regional tectonic fabric and volcanic setting, we do not consider it necessary to try a different strategy for estimating the polarization direction there. Another reason is the quality of the data for KRI and HVD. The records have better SNR compared to the ones from BJA, BLF and LHA, which makes it easier to determine the polarization direction from the particle-motion diagrams.

If the initial polarization of the shear wave at the source is parallel to either the fast or the slow directions, then no splitting will occur. In such a case, the initial particle motion

(a) KRI (b) BJA

Figure 10: Typical examples of events from stations KRI and BJA. The three-component seismograms are shown in the first row and the picks of the S-wave is marked on the east components. The particle motions of the two horizontal components are also displayed starting two samples before the picks. The first five samples after the S-wave arrival are coloured in red. The direction of the first motion is indicated with a green arrow.

After rotation the split shear waves can be clearly identified as shown on the plots at the bottom. The obtained results are N65°E and N37°E for the fast polarization directions and 0.05 s and 0.06 s for time delays for KRI and BJA, respectively.

at the station could point in either the fast or the slow direction. This dependence on the polarization of the waves at the source clould cause misinterpretation of what is considered to be the fast direction when looking at the particle motion diagram.

False reading of the fast direction could also happen due to the complexity of the signals. The first shear-wave arrival might as well be a scattered wave or converted wave.

Long P-wave coda or high-amplitude noise could mix with a small-amplitude shear-wave and disturb the resulting orientation. We have applied criteria to the data set in order to filter out these complications to a certain degree, but still they cannot be excluded as factors affecting the records. Judging just by the initial motion in the polarization diagrams is difficult and can be misleading in certain situations, and for some stations more often than for others.

The experience gained while working with the events and the patterns seen in the particle motion diagrams led to a strategy where the polarization direction is determined over a longer time window than only the very start of the main shear wave. The start of the fast shear wave can be ambiguous. The chosen time window can begin with a part of the P-wave coda. A linear fast wave interfearing coda or noise can appear as elliptical in the particle motion. As already explained, the initial particle motion could be influenced or disturbed by many causes and thus, this strategy gives an alternative interpretation of the ambiguous polarization diagram. This second strategy turns out to give more consistent results for stations BLF and LHA.

We computed a set of synthetic seismograms and plotted the resulting particle motion to study the structure of the patterns that can occur for different time delays, different noise levels and different waveform complexity. Based on this, we determine the possible polarization directions and time delays for the events for BJA, BLF and LHA. The esti- mation of the time delays in this case is more qualitative than quantitative. The window of the waveform that we focus on is usually of length around 0.2-0.3 s after the arrival of the main S-wave. Since the fast and slow split shear waves are nearly orthogonal, two orthogonal directions are determined which approximately describe the main orientations of the motion by visual inspection, starting from the first few samples assumed to be part of the fast shear wave. Here, the time delay is determined by comparing the polarization diagram of a given event and a set of synthetic polarization diagrams. The synthetic data represent waves with different levels of complexity of their waveform, with a varying time delay applied, and with different levels of noise. The complexity of the waveform is varied by including a successively increasing number of lobes of a 10 Hz sinusoid, ranging from

1 to 4. The time delays range from 0.0 s to 0.12 s and SNR ranges between infinity and 1.75. The synthetics and their particle motions are shown in appendix B where a more thorough description of how they are produced is presented. Each polarization diagram is compared with the synthetic data and the plot which describes the real motion best is selected. The shear waves typically have a period around 0.1 s, the same as the waves used to generate the synthetic data, so these real and synthetic diagrams with similar character of motion are expected to have close to equal time delays.

Examples of the application of the second strategy are displayed in figure 11. The polarization directions are determined for two events at BJA and BLF. The particle- motion diagram for the BJA’s event (same as the sample shown for the initial strategy) can be interpreted as having a clear initial direction for the fast split shear wave and 0.04 s later the motion is opened by the arriving slow wave. In an approximately perpendicular direction and after the fast wave has passed, the slow wave writes a similar pattern as the leading wave. The motion character for BLF’s event is more complicated and can be subject of different interpretations. A possible scenario is that the first five samples are written by the leading shear wave, the slow wave changes the direction after that in an almost orthogonal direction, and then a second fast wave arrives, leaving its sign with a linear motion parallel to the first shear-wave arrival. This repetition in the pattern could be due to a multilobed source signal.

(a) BJA

(b) BLF

Figure 11: Determination of the polarization direction for two events at BJA and BLF.

On the diagrams certain directions of the particle motion can be noticed as marked with green lines and assigned S₁ and S₂ notations. These directions are assumed to be the orientations of the fast and slow split shear waves. For both events the angle of S₁ with the vertical happens to be 34°. The small plots on the sides are taken from the synthetic polarization diagrams in appendix B. In the examples for both stations the waveforms have two and three lobes and SNR ≈ 3.75. Based on similarities between the main plots and the respective small plots, and because of similar periods of the recorded shear waves and the synthetic ones, the time delay is estimated to be 0.04 s for BJA’s event and 0.05 s for BLF’s.

5.2 Difficulties in measuring and scattering effects

Late arrivals (after the direct P- and S-wave) of reflected and scattered waves occur due to the fact that seismic-wave energy can take different paths to travel the distance from source to receiver (Lay and Wallace, 1995). Small-scale heterogeneities in the crust lead to considerable scattering of the waveforms of local earthquakes particularly at short peri- ods. After the direct arrival, scattering causes long lasting coda waves which extends the duration of the oscillation (Sato and Fehler, 2012). The coda may include P-wave rever- berations, generated in the unconsolidated near-surface layers, which could hide or make unclear the arrivals of the shear waves. The incident shear waves may be seriously scat- tered by the free-surface where they may endure amplitude and phase changes. Different arrivals and mode conversions could also be generated in this interaction and complicate the recorded signal around the main S-wave arrival. For example, S-to-P converted waves may arrive earlier than the direct S-wave. The coda of the fast shear wave could obscure the arrival of the slow shear wave (Crampin and Gao, 2006). Therefore, the interpretation of the arrival and polarization of the incident shear wave can be hard or even impossible (Crampin, 1990).

In addition, the two split shear-wave components may differ in their waveforms since waves with contrasting polarizations are maybe exposed to different impedances (Crampin and Gao, 2006). The slow split shear waves are likely to include less energy at higher frequencies because they are subject to more attenuation. This causes difficulties in applying automatic techniques for shear-wave splitting analysis.

6 Results

The analysis of shear-wave splitting is based on over 200 events that have occurred in the area of Hengill and Krýsuvík. The following results have been inferred from the orientation of the leading shear-wave polarizations and time delays for each small earthquake.

6.1 Strike of the fractures

Figure 12 displays a comparison between the two approaches for estimation of the fast shear-wave polarization orientation. The rose diagrams represent the frequency of ap- pearance of a certain direction. The left-hand side of the diagrams result from the initial approach and the right-hand side from the second. The dominant directions of fast shear- wave polarization orientations are obvious and clear for all three recording sites when the second approach is applied. The strike groups are in the range of: N30°- 50°E for BJA, N40°- 50°E for BLF and LHA. However, the initial approach, concluding the polarization of the fast shear wave based on the initial motion of the leading wave, gives less consistent results. The orientation at BJA is still mainly in the N40°- 50°E range, but for BLF and LHA the directions are rather scattered and have peaks at N30°- 60°E and N10°- 20°E, respectively.

The diagrams in figure 12b for station BLF look quite different. This is in part misleading because the peaks of the polar histogram are overemphasized on the expanding circle. The orientation perpendicular to the dominant NE direction, which represents the peak in both histograms, may be the slow direction (if the initial polarization at the source was in the slow polarization). Therefore, they are not so different.

The events recorded at stations HVD and KRI result in rose diagrams with narrow distributions. In figure 13 we compare the estimated directions of the fast shear-wave polarization from this study and another, conducted by Jing Liu, MIT, where both studies work with the same initial catalogue of events. The results are compatible and have narrow concentration in the range of N40°- 50°E. The distributions of orientations on figure 13a is slightly rotated clockwise relatively to figure 13b and somewhat wider.

In figure 14 we show the rose diagrams for the fast polarization from our sites together

(a) BJA

(b) BLF

(c) LHA

Figure 12: Comparison between the initial and second approach for retrieving the po- larizations. The results presented are for stations BJA, BLF and LHA. The rose diagrams obtained after applying the initial approach are shown in the left column and the ones from the second approach in the right column.

(a) UU, October 2013 (b) MIT, April 2013

Figure 13: Rose diagrams of the fast direction from earthquakes recorded at KRI. (a) shows the results from this study and (b) displays the results from by Jing Liu, MIT (personal communication, April 2013).

with the results of a previous study in the area by Volti and Crampin (2003). All sites in- dicate a predominantly NE-SW alignment of fast polarization directions. This alignment is parallel to the normal faults and volcanic fissures in the SW Iceland rift zone, i.e. the over-all maximum horizontal stress. This alignment would also be parallel with the pre- dominant microcracks induced by the same over-all stress field. Note that the alignment is more or less the same for sites inside one of the fissure swarms (KRI, BLF, HVD) of the area, on their margins (HEI, BJA and VOS) and for the sites outside (LHA, SOL, SAN) the fissure swarms. An earlier study conducted by Menke et al. (1994) supports these orientations by estimating N40°- 50°E at SOL and a clear main direction of about N30°E at both BJA and HEI.

Figure 14: Map of the Reykjanes Peninsula, showing the polarization rose diagrams for the five stations in this study (red) and from four additional stations in the area from the results of Volti and Crampin (2003) (green). At stations KRI and BJA, the results of three independent studies agree. At stations HEI and SOL, the results of two independent studies also agree quite well.

6.2 Magnitude of anisotropy

All shear-wave splitting measurements are based on events approximately directly beneath the station. We cannot isolate at which depth the anisotropic effect occurs other than between the earthquake and the station. We have relocated the events to refine the depth with a reduced uncertainty of about 10%. By examining the variation of the time delay with depth, we may be able to constrain the depth distribution of anisotropy. This is made difficult in two ways: each station has a limited depth-distribution of seismicity beneath it and the measured time delays are very uncertain. The time delays may very well be systematically underestimated with erroneous time delays inferred from scattered waves arriving before the slow wave.

We have selected those time-delay data which have been given the highest quality 1 = good for those polarization directions that fall within the main peak of the rose diagram in each case. These have been plotted in figure 15 for the individual stations.

Split shear-wave times vary in the interval 10-100 ms. Normally, if the anisotropy and the average velocity were homogeneous in the area of each station, the time delay would be proportional to the depth. A weak positive correlation is seen at stations BJA and HVD suggesting similar behavior, but the depth range of the events is small in both cases. The time delays range from 50-100 ms for depths around 6.5 km up to 9 km for BJA; 10-50 ms for 3-6 km depths below HVD. For a rather small depth range between 3 and 4.5 km in the vicinity of LHA the time delays range from 10 to 70 ms. The data for KRI and BLF suggest, if anything, a decrease of the time delay with depth which does not have any physical meaning. A small range of time delays are measured at BLF (10-50 ms) between depths of 3.5 and 9 km. The time delays at KRI vary in the range 10-100 ms at depths between 2 and 6 km.

The scatter in figure 15 is considerable. This could in part be caused by an azimuthal variation of the strength of anisotropy around each station. More likely, most of the scatter is due to uncertainty of measurement.

Because of limited depth range for each station we have plotted the results from all stations in and around the Hengill area together in figure 16. This is, of course, of limited

meaning unless the depth behavior of anisotropy is similar beneath all those sites. One of the sites (BLF) lies within the Brennisteinsfjöll fissure swarm. The other three (HVD, BJA and LHA) are near the Hengill triple-junction complex. We have no particular reason to expect the anisotropy to vary uniformly beneath those stations. Together their time delays form a crude linear trend with depth indicating weak anisotropy in the top 3 km and on the order of 7% anisotropy below that depth. This is, however, not a physically likely scenario, since the higher pressure and temperature in greater depths cause the cracks to close and heal faster than in shallower depths.

Figure 15: Variation of time delay with focal depth. The plots include just the events with quality 1 according to the initial approach. The number of earthquakes plotted for each station is marked in brackets.

Figure16:VariationoftimedelaywithfocaldepthforthestationsintheHengillarea.Differentcolorsareusedtodistinguishthedata fromthestationsasfollows:fromBLFinred,BJAingreen,HVDinpurple,andLHAinorange.Again,onlytheeventswithcoefficient 1havebeenincluded.

Figure 17 displays a comparison between the statistical patterns of normalized time delays resulting when applying both methods of visual measurements. By normalized time delay we mean the measured time delay divided by the depth of the event. This is a measure of the average strength of anisotropy above the earthquake. The figure shows results for BLF, BJA and LHA and both methods. Fewer data are displayed than the numbers indicated in table 1 (see the number indicated) because in some cases it was not possible to estimate the fast polarization direction due to waveform complexity and the comparison is restricted only to sets of the same events. The initial approach gives a smaller range of normalized time delays. For BJA - 5-11 ms/km, extended to 1-12 ms/km for the second approach. With two exceptions, the initial approach results in a rather consistent, narrow range of small normalized time delays between 1 and 6 ms/km for station BLF, where as the second approach yields the range 2 to 9 ms/km. Most of the normalized times measured at LHA are in the interval 2-8 ms/km for the initial approach and 5-15 ms/km for the second approach. The histograms in figure 17 are similar, although significant differences are evident. Crampin (1994) reports a smaller range of normalized time delays (1.5-4.5 ms/km) for stations in the early SIL network, primarily in the Southern-Iceland Lowlands transform area. Menke et al. (1994), on the other hand, report time delays in the range of 100 - 300 ms, primarily from stations in the same area. They do not present normalized time delays, but their time delays would result in larger normalized time delays than we obtain. This discrepancy may be caused by changes in the stress field with time, if the anisotropy is primarily caused by transient stress-induced microcracks. That does, however, not explain the discrepancy between the estimated time delays of Crampin (1994) and Menke et al. (1994), because they worked with data from more or less the same time period. Perhaps, the most likely explanation is that time delays are quite imprecisely estimated.

Figure 17: Normalized time-delay distribution. Statistical patterns of normalized time delays [ms/km] are shown for BJA, BLF and LHA, to compare the initial (IA) and the second approach (SA) of measurements. The number indicated in brackets shows the number of events for each station.

7 Discussion and conclusions

Shear-wave splitting caused by propagation of shear waves through fractures and stress- aligned microcracks has been clearly identified for microearthquakes recorded over a two- year period in the area of SW Iceland and more specifically around the Hengill and Krýsuvík volcanoes. The seismicity in the region gives a wide choice for selecting events located right below the recording sites. Having a restricted location is a main requirement for shear-wave splitting analysis due to a strong distortion of the signals when the incidence angles are far from the vertical. All earthquakes are subject to relocation in order to constrain the depth parameter. The study has been limited to earthquakes within the 30°-window, with a minimum threshold applied for the SNR, with clear phase arrivals, and with a maximum threshold applied to the vertical-component amplitude at the onset of the shear waves compared to the horizontal components in order to reject data with strong scattering effects. We observe a prevailing polarization direction of the leading shear wave in the range of N30°- 60°E for all stations. This is consistent with the tectonic fabric of the area as manifest in the orientation of normal faults, volcanic fissures and dikes.

The average polarization direction at all stations agree with three previous studies in the same area. This orientation shows the average, over-all direction of the maximum horizontal compressive stress in SW Iceland. These studies span 20 years. Thus, the orientation of the stress field appears not to have varied, even locally at the individ- ual recording sites. The alternative interpretation is that the anisotropy is not strongly affected by stress-aligned microcracks, but rather controlled by the long-term stable ex- tensional stress field through long-lived and relatively large-scale fractures.

Visual techniques for measuring shear-wave splitting are often considered to be sub- jective, since they involve visual estimation of the polarizations and visual correlation be- tween the waveforms of fast and slow shear waves to determine the time delay (Crampin and Gao, 2006). The two visual approaches used here are similar, but involve a difference in how the particle-motion pattern is viewed. The initial approach focuses on detail - the behavior of the initial motion at the beginning of the polarization diagram, while the

second approach is based on the over-all motion. Neither method excludes the possibility of wrong interpretation, but the consistency of the fast polarization direction obtained with the second approach gives confidence in those results.

As much as the NE-SW alignment of the polarization directions is predominant, the time delays measured for the events do not give a solid base for drawing conclusions. Weak trends of increase of the time delay as a function of depth are visible for stations BJA, LHA and HVD, but BLF and KRI show, if anything, the opposite. Comparable results for time delays at KRI are reported by Jing Liu, MIT. This may indicate that the scatter in the results is not simply caused by measurement difficulty and reflects strongly laterally- varying anisotropy or strong common scattering effects. The amount of data for BLF is rather small and they cannot be argued to give a robust base for conclusions. An increase of the time delays with depth becomes visible when data from all three stations in the Hengill region are superimposed (figure 16). This would be consistent with an anisotropic layer starting at about 3 km depth and extending to at least 9 km. However, we have no reason to expect the strength of anisotropy to be the same or similar beneath all three stations contributing to this figure. Furthermore, it would be difficult to reconcile such a layer with increasing pressure successively closing more of the fractures with increasing depth and increasing temperature shortening the life time of stress-induced microcracks.

The only robust conclusion that can be drawn from our measurements of time delay is that they are smaller at station BLF than at other stations with events at the same depths, indicating that the strength of anisotropy beneath this site is smaller than at the other sites.

In previous studies of shear-wave splitting the problem of earthquake prediction is widely discussed. It is claimed that variation in measured splitting times before the occurrence of an earthquake reflect changes in the stress field and therefore could help to forecast the event in the short term. Based on our own measurements and the difficulties met in particular in measuring the time delays, we find it difficult to believe that such changes can be resolved. Our observations indicate that the difference between the split shear waves is hard to measure with confidence. The complex character of the data is

an obstacle for both visual and automatic estimations. It is often possible to identify the initial polarization direction, but the time delays are much more difficult to measure because the slow shear wave can be obscured by scattered energy in this highly fractured rift zone. If the available catalogue is limited to a small number of earthquakes then no firm conclusions can be drawn from the uncertainly-measured time delays.

Acknowledgements

I would like to thank Prof. Ólafur Gudmundsson, my supervisor, for his guidance, en- couragement and dedicated time throughout the entire work. I am very grateful that he shared his knowledge and valuable experience with me. By giving me this topic he gave me the chance to work on something different and expand my horizons in seismology.

My report would not be the same if it was not for the insightful comments and im- pressive editing of Karin Berglund, my precious friend, always giving me her strong moral support, and also my opponent on the thesis defense.

I acknowledge the help I received from: Jing Liu and Michael Fehler from MIT, who shared their results for comparison, Ari Tryggvason who gave detailed data about the used velocity model, and Zeynab Jeddi who provided the earthquake records.

I am very happy that I had the chance while writing my thesis to work with the wonderful people in the Swedish National Seismological Network. My special thanks goes to: Björn Lund for his support and shared seismological wisdom, Conny Holmqvist - my coach in picking, Hossein Shomali - for helping me to find my way with SAC, my boss Reynir Böðvarsson - for giving me the opportunity to study and work at the same time, and also Arnaud Pharasyn - for making a better programmer out of me.

I can hardly express my gratitude to my parents and incredible sister for always being there for me, supporting all my decisions and giving me their love. I am thankful to all my friends for giving me strength and keeping my spirits high. I am especially grateful to Fatima Afsar and Plamena Raykova for their invaluable friendship, professional understanding and encouragement during the past few months. This journey would not be the same without the love and patience of my amazing Fartash, thank you for believing in me and giving me your shoulder to lean on!

References

Booth, D. C. and Crampin, S. (1985). Shear-wave polarizations on a curved wavefront at an isotropic free-surface. Geophysical Journal of the Royal Astronomical Society, 83:31–45.

Chen, T. C., Booth, D. C., and Crampin, S. (1987). Shear-wave polarizations near the North Anatolian Fault - III. Observations of temporal changes. Geophysical Journal of the Royal Astronomical Society, 91:287–311.

Crampin, S. (1990). The Scattering of Shear-waves in the Crust. Pure and Applied Geophysics, 132.

Crampin, S. (1994). The fracture criticality of crustal rocks. Geophysical Journal Inter- national, 118:428–438.

Crampin, S., Evans, R., and Atkinson, B. K. (1984). Earthquake prediction: a new physical basis. Geophysical Journal of the Royal Astronomical Society, 76:147–156.

Crampin, S. and Gao, Y. (2006). A review of techniques for measuring shear-wave split- ting above small earthquakes. Physics of the Earth and Planetary Interiors, 159:365–

385.

Crampin, S. and Lovell, J. H. (1991). A decade of shear-wave splitting in the Earth’s crust: what does it mean? What use can we make of it? And what should we do next?

Geophysical Journal International, 107:387–407.

Crampin, S. and Zatsepin, S. V. (1997). Modelling the compliance of crustal rock: II - response to temporal changes before earthquakes. Geophysical Journal International, 129:495–506.

Gao, Y., Zheng, S. H., and Sun, Y. (1995). Crack-induced anisotropy in the crust from shear-wave splitting observed in Tangshan region, North China. Acta Seismologica (English Edn), 8:351–363.

Goldstein, P., Dodge, D., Firpo, M., and Minner, L. (2003). SAC2000: Signal processing and analysis tools for seismologists and engineers, Invited contribution to “The IASPEI International Handbook of Earthquake and Engineering Seismology. Academic Press,

"London". Edited by WHK Lee, H. Kanamori, P.C. Jennings, and C. Kisslinger.

Jousset, P., Haberland, C., Bauer, K., and Arnason, K. (2011). Hengill geothermal vol- canic complex (Iceland) characterized by integrated geophysical observations. Geother- mics, 40(1):1–24.

Lay, T. and Wallace, T. C. (1995). Modern Global Seismology. Academic Press.

Menke, W., Brandsdottir, B., Jakobsdottir, S., and Stefansson, R. (1994). Seismic Anisotropy In the Crust At the Mid-Atlantic Plate Boundary In South-West Iceland.

Geophysical Journal International, 119(3):783–790.

Rial, J. A., Elkibbi, M., and Yang, M. (2005). Shear-wave splitting as a tool for the charac- terization of geothermal fractured reservoirs: lessons learned. Geothermics, 34:365–385.

Rutter, E. H. and Brodie, K. H. (1991). Lithosphere rheology - a note of caution. Journal of Structural Geology, 13:363–367.

Sato, H. and Fehler, M. C. (2012). Seismic Wave Propagation and Scattering in the Heterogeneous Earth. Springer.

Tryggvason, A., Rögnvaldsson, S. T., and Flovenz, O. G. (2002). Three-dimensional imaging of the P- and S-wave velocity structure and earthquake locations beneath Southwest Iceland. Geophysical Journal International, 151:848–866.

Volti, T. and Crampin, S. (2003). A four-year study of shear-wave splitting in iceland: 1.

background and preliminary analysis. Geological Society, London, Special Publications, 212(1):117–133.

Zatsepin, S. V. and Crampin, S. (1997). Modelling the compliance of crustal rock: I - response of shear-wave splitting to differential stress.Geophysical Journal International, 129:477–494.

Appendices

A Free surface correction

The recordings of the signals at the seismic station depend on the reflection and the conversion at the free surface, but also are affected by the topography and the anisotropy along the ray path. The free surface is the strongest disruption that interrupts a seismic signal. Its importance is due to the placement of the stations at the free surface. In the calculations of the incidence angle, a correction for the free surface has been included. It is derived, based on the following theory, equations and assumptions.

Figure 18: Geometry of the incidence and resulting rays.

In this case, a plane wave has been considered, incident on a boundary between a medium and a free surface as shown in figure 18. The particle motion is assumed to have no lateral variations and to be just in the xz-plane. The scalar potentials of the incidence and reflected P- waves (φ_i and φ_r respectively) and the vector potential of the S-wave (ψ_r) are written as follows:

φ_i = A exp(iω(t − px − ηz)), (9)

φ_r = B exp(iω(t − px + ηz)), (10)

ψ_r = C exp(iω(t − px + ξz)), (11)

where A, B,and C are the respective amplitudes of the waves, p is the horizontal slowness, η and ξ are the vertical slownesses of the P- and S-waves, respectively. α and β are the P- and the S-wave, respectively. Using the velocities and the incidence angle i, the three parameters p, η and ξ can be expressed as:

p = sin i

α , (12)

η = r 1

α² − p² = cos i

α , (13)

ξ =r 1

β² − p² = cos j

β . (14)

According to the Helmholtz’s theorem, every vector field u can be presented in terms of a scalar potential φ and a vector potential ψ by

u = ∇φ + ∇ × ψ, (15)

if φ is curl free (∇ × φ = 0) and ψ is divergence free (∇ · ψ = 0) (Lay and Wallace, 1995).

The displacement components of the three waves are

(u, 0, w) = ∇(φ_i+ φ_r) + ∇ × (0, ψ, 0), (16)

from where the displacements u and w can be written as:

u = ∂φ_i

∂x + ∂φ_r

∂x − ∂ψ_i

∂z w = ∂φ_i

∂z +∂φ_r

∂z + ∂ψ_i

∂x

(17)

which by substituting the potentials and after some algebra become:

u = −iω[pA exp(−iωηz) + pB exp(iωηz) + ξC exp(iωξz)] exp(iω(t − px)) w = −iω[ηA exp(−iωηz) − ηB exp(iωηz) + pC exp(iωξz)] exp(iω(t − px)).

(18)

Tractions at the free surface must vanish, i.e. σ_zz and σ_xz must be equal to 0 at z = 0.

In general, these stress components are expressed as:

σ_zz = λ ∂u

∂x +∂w

∂z
+ 2µ∂w

∂z (19)

σ_xz = µ ∂u

∂z +∂w

∂x
. (20)

Applying the boundary conditions on (19) and (20) for z = 0 and using the expressions for u and w in (18), we get the following equations for the amplitudes and the slownesses:

(3η²+ p²)(A + B) − 2pξC = 0 2pη(A − B) + (p² − ξ²)C = 0,

(21)

assuming a Poisson solid λ = µ. From here, the sum and difference between the amplitudes of the incidence and reflected P-waves can be expressed, such as:

A + B = 2pξC 3η²+ p² A − B = −C(p²− ξ²)

2pη .

(22)

For the displacements at the surface (z = 0) we the expressions:

U = −iω[p(A + B) + ξC] exp(iω(t − px)) W = −iω[η(A − B) + pC] exp(iω(t − px)).

(23)

A parameter, R, is defined as the ratio between the displacements at the free surface of the horizontal and vertical components:

R = U

W = p(A + B) + ξC

η(A − B) + pC. (24)

Substituting (22) in (24) gives the following simplified relation between R and the three slownesses:

R = 2pξ

3η²+ p², (25)

again using the assumption of a Poisson solid. This expression can also be written in

terms of the incidence angle after substitution using (12), (13) and (14):

R = 2 sin ip

3 − sin²i

3 − 2 sin²i . (26)

A parameter s is defined as s = sin i, which can be isolated using (26) gives an expression for s(R):

s = r3

2 s

1 − 1

√R²+ 1. (27)

This is valid for a Poisson’s ratio σ = 0.25 (α/β = √

3,λ = µ). This expression is used to estimate the incidence angle i based on the apparent incident angle according to R. Thus, the incidence angle is corrected for the effect of the free surface.

In case the ratio between the seismic velocities is different, i.e. α/β = γ, then:

s

λ + 2µ

ρ = γr µ

ρ ⇒ λ = (γ²− 2)µ, (28)

where λ and µ are the Lamé parameters, and ρ is the density of the medium. Substituting in the first boundary condition, where (19) is set to 0, this expression for λ and the displacements in (18), after some calculations, we get the relation for the sum of the amplitudes of the incident and reflected P-waves as follows:

A + B = 2pξC

γ²η²+ p²(γ²− 2). (29)

The equation for the difference of A and B will stay unchanged.

Using the same steps as before by substituting the expressions for the new A + B from (29) and A − B from (22), the general expression for R is derived:

R = 2 sin ip

γ² − sin²i

γ²− 2 sin²i . (30)

Correspondingly, s(R) is:

s = γ

√2 s

1 − 1

√R²+ 1. (31)

B Synthetic data

The following set of synthetic data (figures 19 to 30) presents the evolution of the particle motion of two waveforms (with a period of 0.1 s) for different time delays. The waveform is a multilobed signal with a varying number of lobes - from 1 to 4, and every lobe is a sinusoidal wave. The difference between the onsets of the waves is increased from 0.0 s to 0.12 s, or from no time delay up to a delay bigger than the wave period or the length of the waveform included (for the one-lobed signal). All seismograms and particle motion diagrams are shown within 0.4 s - window and the first 5 samples of the sinusoidal waveform are colored in red in the particle-motion diagrams.

In aadition to the number of lobes and the varying time delay, the SNR is also changed.

Three SNR-cases are displayed: without any added noise, with some added noise with a SNR ≈ 3.75, and with added noise of a higher amplitude corresponding to SNR ≈ 1.7.

The SNR intervals are chosen to reflect the typical quality of the seismic signals used in this study - above and below the SNR = 2.0 - threshold. The artificial noise is a low- pass filtered white Gaussian noise, band-pass filtered to simulate the noise characteristics encountered in this study. The noise added to the slow component is independent of the noise to the noise added to the fast component.

The synthetic data gives an opportunity to trace how the particle motion develops with increasing time delay between the two split shear waves. The plots show that there is a clear fast wave which on the figures is the initial horizontal particle motion to the right.

This evolves to a round-shaped motion and ends with a vertical motion representing the slow wave after the leading wave has passed. It is clear that the simulated, higher-noise levels obscure the fast and the slow polarization directions and in cases of small time delays and short waveform duration they may be completely hidden in the noise.

Figure 19: Synthetic data for one-lobed waveform without added noise.

Figure 20: Synthetic data for one-lobed waveform with added noise (SNR ≈ 3.71).