The model for the earthquake sequence at the SISZ

During PRENLAB-2, this model was further improved. The main features of this model are given again to ease comparison with the new results.

The method

The forward modelling of stress fields is done by applying static dislocation theory to geodetic data and data obtained through the seismic moments from seismograms. It allows to calculate displacements, strain and stresses due to double-couple and extensional sources in layered elastic and inelastic earth structures. Besides the change in displacement during the event, the changes caused by the movement of plates are included (for further details see e.g. Roth 1989).

Usually, for earthquake hazard estimation, the location, the magnitude and the statistically estimated recurrence period of former events is used. To improve this, here the rupture length and width as well as the tectonic setting and the crustal deformation rates are considered while calculating the space time development of the stress field.

The targets

In general, with these models, Subpart 7B aims:

$\bullet$

to achieve a better understanding of the distribution of seismicity in space and time, its clustering and migration in Iceland.

$\bullet$

to provide models for the joint interpretation of the data gathered in the whole research programme, of which this is one part.

$\bullet$

to compare models of stress fields at SISZ to those for stress fields in other regions, e.g. the North Anatolian fault zone.

$\bullet$

to make a contribution to the intermediate-term earthquake prediction in this populated and economically important region of Iceland.

The tectonic setting

The SISZ is situated between two sections of the mid-Atlantic ridge, the Reykjanes ridge (RR) and the eastern volcanic zone (EVZ). Even though the angle between the SISZ and the neighbouring ridges is far from 90$^\circ$, it is considered as a transform fault. Following the transform fault hypothesis, left-lateral shear stress is expected along the E-W striking zone. This is equivalent to right-lateral shear stress on N-S oriented rupture planes. In fact, earthquakes seem to occur on N-S trending en-echelon faults (cf. Einarsson et al. 1981; Hackman et al. 1990), They are located side by side between the Hengill-Ölfus triple junction, where the RR meets the low activity western volcanic zone (WVZ) and Hekla volcano, in the EVZ (cf. Einarsson et al. 1981) (Figure 50). As we further know from Subprojects 4 and 5, the orientation of the larger horizontal principal stress is NE-SW, i.e. fits to an active N-S or E-W trending fault, which is (at least for the period of those investigations) not a weak fault like the San Andreas fault (cf. Zoback et al. 1987). Moreover, the stress orientation seems to have been constant since Pliocene time.

  

Figure 50: The South Iceland seismic zone showing mapped surface breaks and regions in which over half of the buildings were destroyed in historic seismic events (after Einarsson et al. 1981). The north-south dashed line near Vatnafjöll indicates the estimated location of the fault on which the May 25, 1987, earthquake occurred (after Bjarnason and Einarsson 1991). The structural features and the coastline are after Einarsson and Sæmundsson (1987).

In detail, the questions to be solved are:

$\bullet$

Do these events, placed on parallel faults, release all the energy stored in the 3-D volume of the SISZ?

$\bullet$

Do the earthquakes always take place in areas of high stress?

$\bullet$

What is the critical stress level? How large is its variability?

$\bullet$

Where are the highest stresses nowadays?

  

Figure 51: Map of Iceland and surrounding area. Thick red lines indicate mid-Atlantic ridge segments, as used in the models. The smaller box shows the region of the model on the SISZ. The SISZ extends approximately between 21.4$^\circ$W, 63.95$^\circ$N to 18.8$^\circ$W, 64$^\circ$N. The large box gives the region for the Iceland rift model. RR: Reykjanes ridge, KR: Kolbeinsey ridge, WVZ/EVZ: western/eastern volcanic zone, SISZ: South Iceland seismic zone, TFZ: Tjörnes fracture zone.

The area investigated extends from 18$^\circ$ to 24$^\circ$W and from 63$^\circ$ to 65$^\circ$N. The origin is set to 24$^\circ$W, 64$^\circ$N (cf. Figure 51) it includes the SISZ, $\pm 1^\circ$north and south of 64$^\circ$N, the SW edge of the EVZ, and the northeasternmost part of the RR.

The initial stress field

The initial stress field is determined as follows: A tensional stress acting N103$^\circ$E (nearly parallel to the SISZ; cf. DeMets et al. 1990) is assumed, due to the opening of the mid-Atlantic ridge in the region adjacent to the transform fault. While this rifting induces mainly shear stresses in the SISZ with a small opening component, the rift segments (RR and EVZ) are modelled with tensional stress and a small shear stress contribution. Tensional stresses at both ridges are modelled as constantly being released to end up with zero values at the rifts. This induces additional stress in the transform zone. The stress magnitude, which is unknown, is set to a value that produces left-lateral shear stresses in E-W direction as large as the stress drop determined for the largest event (M=7.1) in the studied earthquake sequence.

On this initial field, the stress changes due to earthquakes are iteratively superposed as well as the stress changes due to further spreading at the ridge segments. From global geodetic measurements an opening of 2 cm/year is found, e.g. in DeMets et al. (1990). This was used as a zero-th order approach but was reduced to only 1 cm/year, as discussed later. Further, as the simplest assumption, lacking other data, the spreading rate is taken to be constant during the modelled time period, even though this can be questioned as for instance the present debate on the stress increase in the New Madrid seismic zone shows (cf. Schweig and Gomberg 1999; Newman and Stein 1999) The stress field before every event is thus the sum of the initial field, the stress drop of all preceding events, and the plate tectonic stress build-up since the starting time of the model, which is set to 1706, when the first event in the series occurred.

Results were calculated for 280x220 test-points covering 280 km in E-W direction and 220 km in N-S direction. Stresses were computed for a homogeneous half-space, as a starting model. Although surface stress changes are calculated, these should be representative for crustal stresses using values for the moduli, that are typical for oceanic crust (see Dziewonski et al. 1975) and not for sedimentary layers at the surface. Moreover, as the faults are introduced vertically into the unlayered environment, the stresses do not vary much with depth, besides at the lower end of the fault.

Changes and improvements in PRENLAB-2

In the first phase of PRENLAB-2, the models developed in PRENLAB-1 were improved:

A

The test-point density was increased from 56x44 (5 km distance) to 280x220 (1 km distance) to get more details of the stress field and to reduce interpolation errors.

B

At the western end of the SISZ, segments with aseismic oblique slip (mainly normal faulting with a smaller component of left-lateral strike-slip) were introduced, to better fit the Reykjanes ridge (RR) between the southwest tip of the Reykjanes peninsula to Hengill-Ölfus triple junction (Figure 51).

C

At the eastern and western tips of the SISZ two areas with steady stress release were introduced (see dashed lines in Figure 52). This could occur by a high rate of small events and maybe by creep. This is likely as these areas did never show strong events (M$\ge$6) in the seismic history of Iceland, with the exception of possible events in 1546 and 1632 at the western end (about 21.3$^\circ$W) and one event in 1311 (about 18.9$^\circ$W) (cf. Halldórsson 1991). See under Task 2 for a discussion of this.

D

To investigate the model resolution a set of different models is produced: Besides the main model, several extreme cases are assumed and the variation of the main results under these assumptions is observed.

E

The stress field was extrapolated to 1999 and was now updated to the situation after the June 21, 2000, event.

  

Figure 52: The box gives the area in SW-Iceland used in the modelling as indicated with the small box in Figure 51.

The earthquake data

All events with M $\ge$ 6 since 1706 were used (listed in Table 8., following Halldórsson et al. 1984; after Hackman et al. 1990; Stefánsson and Halldórsson 1988; and Stefánsson et al. 1993). The catalogue is supposed to be complete from 1706 for these earthquakes.

All ruptures were set to be oriented N-S, according to the isolines of damage intensity and surface ruptures shown in Figure 50. As only the events in 1912 and 2000 were instrumentally recorded, the source parameters are not very accurate - a problem to be discussed further below. The position of the epicenters and the rupture planes used in the models are given in Figure 53 and Figure 54, respectively.

 

Table: Earthquakes M$\ge$6 since 1706 in the South Iceland seismic zone. 1 ) Data taken from Stefánsson et al. (1993); for the events in 1706, 1732, and 1734 no exact date is known. Data on the events of June 2000 are from Stefánsson, Guðmundsson and Halldórsson (pers. comm.), with magnitudes as moment magnitudes. 2 ) Position in the model coordinate system with origin at 64$^\circ$N, 24$^\circ$W. 3 ) Calculated via the magnitude moment relationship log M₀ [dyne cm] = 1.5M_S+(11.8 - log( $\sigma _a/\mu$)) with the apparent stress $\sigma _a$ = 1.5 MPa and the shear modulus $\mu$ = 39 GPa (after Kanamori and Anderson 1975), followed by using the values of $\mu$ above, the rupture length as given in the table as well as a vertical fault width of 14 km east of 21$^\circ$W and 7 km between 21$^\circ$W and 21.2$^\circ$W. Finally, the values for all events before those in 2000 were reduced by a factor of 2, following the discussion of Hackman et al. (1990). All this does not apply to the June 17 and 21, 2000, earthquakes, for which good instrumental data are available, and e.g. a maximum rupture depth of 9 and 7 km were given, respectively. 4 ) Besides for the June 2000 events, calculated using log L [km] = 0.5 M - 2 (after Qian 1986) which results in slightly lower values compared to e.g. Schick (1968).

Date1	Magnitude1	Epicenter1		South end of		Co-seismic	Rupture
				rupture2		slip3	length4
		Lat. $^\circ$N	Long. $^\circ$W	x [km]	y [km]	U0 [m]	L [km]
1706	6.0	64.0	21.2	131	-5	0.30	10
1732	6.7	64.0	20.1	183	-11	0.77	22
1734	6.8	63.9	20.8	150	-23	0.96	25
14.08.1784	7.1	64.0	20.5	164	-18	1.9	35
16.08.1784	6.7	63.9	20.9	145	-22	0.77	22
26.08.1896	6.9	64.0	20.2	178	-14	1.2	28
27.08.1896	6.7	64.0	20.1	183	-11	0.77	22
05.09.1896	6.0	63.9	21.0	140	-16	0.30	10
05.09.1896	6.5	64.0	20.6	159	-9	0.48	18
06.09.1896	6.0	63.9	21.2	131	-16	0.30	10
06.05.1912	7.0	63.9	20.0	187	-27	1.5	32
17.06.2000	6.5	64.0	20.4	169	-12	0.9	16
21.06.2000	6.4	64.0	20.7	154	-13	1.1	18

 

  

Figure 53: Location of the earthquakes in space and time (cf. Table 8). As the events have been located on N-S trending faults in an E-W trending fault zone, their location is very accurately displayed in this graph.

  

Figure 54: The position of the rupture planes of the earthquake sequence as used in the models (cf. Table 8).

The results

The stress fields at 20 dates were calculated: the pre- and post-seismic situation for all 13 events. The time before 6 events was too short to accumulate appreciable plate tectonic stresses since the preceding event. In these cases, the post-seismic stress field of the preceding event was set equal to the pre-seismic stress field of these earthquakes.

  

Figure 55: The stress field after the May 6, 1912 earthquake.

Originally, an extrapolation was done from the 1912 earthquake to spring 1999. After the two earthquakes happened this June, it was updated to June 17, 2000 (Figure 56) and the effect of both events were determined (see Figure 57 and Figure 58).

  

Figure 56: The stress field before the June 17, 2000, M=6.5 earthquake occurring at 20.4$^\circ$W, 64.0$^\circ$N, i.e. (169, -12 -- +4).

  

Figure 57: The stress field after the June 17, 2000, M=6.5 earthquake occurring at 20.4$^\circ$W, 64.0$^\circ$N), i.e. (169, -12 -- +4).

  

Figure 58: The stress field after the June 21, 2000, M=6.4 earthquake occurring at 20.7$^\circ$W, 64.0$^\circ$N, i.e. (154, -13 -- +5).

As a simple assumption, one might expect, that earthquakes in a certain fault zone usually occur at about the same critical shear stress level. We examine here if such an expectation matches the known facts about the earthquakes, given above, and the stress field in the SISZ from knowledge about plate motion and the modelling here. Figure 59 summarizes the mean shear stress level before each of the earthquakes at the area of the impending event. The stress level is near the average (1.8 MPa) or higher for most of the events. The highest value (for 1896e) is mainly influenced by the fact that the rupture area for the event in 1706 was located completely north of that in 1896 (cf. Figure 54). The second 1784 earthquake (two days after the first, 0.4 in magnitude smaller, 19 km away) might have been an aftershock and therefore situated in a lower stress area (1.6 MPa). The same might apply to the second of the 1896 events. The forth event in 1896 and the first shock in 2000 are both influenced by the largest event in the series, i.e. the first one in 1784. Thus, the accuracy of the source parameter of this 1784 earthquake strongly influences the whole model. Concerning the first event in June 2000, it has to be noted that it took place in a very inhomogeneous stress field, i.e. there are low stresses in the north of the rupture plane and high ones in the south (see Figure 56). Calculating an average pre-seismic stress level might be especially misleading for this event, when all the test-points around the rupture plane are considered equally.

  

Figure 59: Cross plot of the pre-seismic shear stress level at the site of the impending earthquakes vs. occurrence time. The stress values at 2 to 14 test-points at 0.5 km distance to the surface trace of the rupture plane were averaged. Instead of a plate velocity of 2 cm/year only 1 cm/year was used (see text). Letters "a" through "e" denote the events in one year in temporal sequence. The dashed line gives the average stress, the solid line represents the trend found by a linear least squares fit.

Checking the performance of the model in a qualitative way, we examined if the earthquakes hit the high stress area and how large the high stress areas with no event were at the same time (the range in longitude with high stresses was summed when the N-S extension of the area was at least 5 km and the longitude range for the event that occurred was subtracted, usually 0.1 to 0.2 degrees in longitude). The results can be found in Table 9 and are quite satisfying with respect to the named question. Almost all events hit high stress areas and the size of high stress areas with no event was rather small after the earthquakes in the eighteenth century. Nevertheless, the problem remains, why some events did not occur earlier (at lower stress), just passing the "limit in pre-seismic stress" (i.e. here: the average pre-seismic stress).

 

Table 9: Qualitative evaluation of the model performance. 1 ) As the initial stress field is unknown, it was assumed to be homogeneous. So the first series of events that ruptures all across the fault zone is very strongly influenced by this assumption. 2 ) The stress field before the first earthquake in June 2000 is very inhomogeneous, i.e. there are low stresses in the north of the rupture plane and high ones in the south (see Figure 56). Looking at the development of the stress distribution, this situation is still influenced by the rupture position and magnitude of the first event in 1784 (cf. Figure 54). These historical data have some uncertainty (cf. Figure 50).

Event	Pre-seismic	Earth-	Hit high	Size of other
time:	stress high	quake	stress	areas with	Remarks:
	at $^\circ$W:	at $^\circ$W:	area?	high stress:
1706	19.8 - 21.2	21.2	yes	large	tuning phase1
1732	19.8 - 21.1	20.1	yes	large	tuning phase
1734	19.8 - 20.0,	20.8	yes	large	tuning phase
	20.2 - 21.1
	19.8 - 20.0,
1784a	20.2 - 20.7,	20.5	yes	medium	tuning phase
	20.8 - 21.1
1784b	19.8 - 20.0,	20.9	yes	small	tuning phase
	20.9 - 21.1
1896a	19.8 - 20.3,	20.2	yes	medium	other high stress
	20.8 - 21.2				area hit 10 days later
	19.8 - 20.0,				other high stress
1896b	$\pm$ 20.6,	20.1	no	medium	area hit 9 days later
	20.8 - 21.2
	19.8 - 20.0,				other high stress
1896c	$\pm$ 20.6,	21.0	yes	medium	area hit next day
	20.8 - 21.2
	19.8 - 20.0,				other high stress
1896d	$\pm$ 20.6,	20.6	yes	medium	area hit at the same
	20.8 - 21.2				day
1896e	19.8 - 20.0,	21.2	yes	small
	20.8 - 21.2
1912	19.8 - 20.0,	20.0	yes	small
	20.8 - 21.2
2000a	$\pm$ 19.9,	20.4	no/	small	hit short (N-S)
	20.5 - 21.2		yes		high stress area2
2000b	$\pm$ 19.9,	20.7	yes	small
	20.6 - 21.2

 

As stated earlier, the plate velocity used in the model was set to 1 cm/year, only half of that what is measured. The reason can be seen in Figure 59 in comparison with Figure 60. The higher plate velocity leads to a steady and strong increase in the pre-seismic stress level, which is very unlikely, while the reduced rate entails an almost constant level. As the model cannot change the geodetic results, the reason for the discrepancy might be that only 50% of the stress build-up by plate tectonics is released seismically. For the rest, stable sliding or aseismic creep could be responsible.

  

Figure 60: Cross plot of the pre-seismic shear stress level at the site of the impending earthquakes vs. occurrence time. A plate velocity of 2 cm/year was used. For further explanations see Figure 59.

Two problems were addressed next: (i) To see how sensitive the results depend on the model parameters, it was begun to check extreme cases and their outcome. (ii) The average stress level before some events is considerably lower than for others (cf. the discussion on Figure 59). This is, among other reasons, due to the fact that the rupture planes, used until now, extend rather far to the north and south of the SISZ. The damage areas from historical records are not gathered by scientists and are usually biased by uneven population density. So the magnitudes and locations are not very accurate, as stated earlier. And as mentioned in the footnotes of Table 9, there are doubts on the correct rupture size from global relations between magnitude and rupture length.

From both reasons given here, a model was calculated that uses the same seismic moment of the events, but cuts the fault length to 50% while doubling the co-seismic displacement. It will be termed "short rupture model". One side-effect of this change is an increase of the stress level, as the moment release is concentrated to a smaller area. The initial stress field amplitude was increased accordingly, because - as described above - this field is adjusted to the average stress change of the strongest event. It is important to note that the increase in stress level does not change the stress pattern of the initial stress field; as we are not looking for specific stress amplitudes but for stress concentrations, the change in level is not important. The resulting variation of the pre-seismic stress level is - as expected - smoother than before due to the concentration of stress release to high stress areas. However, this approach could only be used, once the historical events are re-evaluated with the result of shorter ruptures.

The processing of the recent earthquakes last June led to maximum rupture depth of less than ten kilometres. Even though there are smaller events located down to 13 km (cf. Stefánsson et al. 1993), the assumption that all ruptures extend to no more than a depth of 10 km seems to be reasonable. Such a model was calculated too, replacing the maximum depth of 14 km for most events (cf. Table 9) by 10 km. In this case too, the stress release by the events is higher, as it is concentrated to an area closer to the rupture plane. Moreover, the interaction of the events is lower due to this concentration in space. However, the stress level before the main events remains in a similar range as before (average now 2.1 MPa instead of 1.8 MPa), and the variation in the pre-seismic stress does not differ much (the standard deviation is 0.87 MPa instead of 0.88 MPa) from the model with deeper ruptures (cf. Figure 61).

  

Figure 61: Cross plot of the pre-seismic shear stress level at the site of the impending earthquakes vs. occurrence time. A maximum rupture depth of 10 km was used. For further explanations see Figure 59.

In general, the variation of the model parameters shows that the models are stable, i.e. small changes in the parameters do not lead to totally different results. Therefore, they fulfill this important quality criterium.