Robust features of the source process for the 2004 Parkﬁeld, California, earthquake from strong-motion seismograms

We explore a recently developed procedure for kinematic inversion based on an elliptical subfault approximation. In this method, the slip is modelled by a small set of elliptical patches, each ellipse having a Gaussian distribution of slip. We invert near-ﬁeld strong ground motion for the 2004 September 28 M w 6.0 Parkﬁeld, California, earthquake. The data set consists of 10 digital three-component 18-s long displacement seismograms. The best model gives a moment of 1.21 × 10 18 N m, with slip on two distinct ellipses, one with a high-slip amplitude of 0.91 m located 20 km northwest of the hypocentre. The average rupture speed of the rupture process is ∼ 2.7 km s − 1 . We ﬁnd no slip in the top 5 km. At this depth, a lineation of small aftershocks marks the transition from creeping above to locked below, in the interseismic period. The high-slip patch coincides spatially with the hypocentre of the 1966 M w 6.0 Parkﬁeld, California, earthquake. The larger earthquakes prior to the 2004 Parkﬁeld earthquake and the aftershocks of the 2004 earthquake ( M w > 3) also lie around this high-slip patch, where our model images a sharp slip gradient. This observation suggests the presence of a permanent asperity that breaks during large earthquakes, and has important implications for the slip deﬁcit observed on the Parkﬁeld segment, which is necessary for reliable seismic hazard assessment.


I N T RO D U C T I O N
The aim of inversions for the earthquake source process is to find the slip distribution history that produces the best fit to the recorded ground motion. The usual approach in kinematic inversions is to subdivide the fault plane into rectangular subfaults, an approach first taken by Trifunac (1974). In each of these subfaults, source parameters are then retrieved. However, this approach has the disadvantage of requiring a large number of subfaults to adequately model the fault plane, introducing non-uniqueness into the problem (Das & Kostrov 1990, 1994, which is dealt with by using additional constraints. In this paper, we investigate a recently developed procedure for kinematic inversion involving elliptical subfaults, which, in addition to creating an intrinsically smooth slip distribution inside the slipping region, has the advantage of reducing the number of parameters to be inverted for. We use only a small number of ellipses, each of which is described by seven parameters. To test this approach with near-field records, we analyse the strong-ground motion data for the 2004 September 28 M w 6.0 Parkfield, California earthquake to obtain its robust source features. The Parkfield segment is part of the San Andreas transform fault system which accommodates right-lateral tectonic motion between the Pacific (PAC) and the North-American (NAM) plates ( Fig. 1). Bounded by a creeping section to the northwest and a locked section to the southeast (the Cholame segment), this region experienced five earthquakes of about M w 6.0 since 1881 (e.g. Bakun & McEvilly 1979;Toppozada et al. 2002;Smith & Sandwell 2006). Knowledge of the occurrence of the 1934 and the 1966 Parkfield earthquakes and the suggestion that the 1966 earthquake was an almost exact repeat of the 1934 earthquake, led to the setting up of the Parkfield Prediction Experiment (Bakun & Lindh 1985). As a consequence, the 2004 earthquake at Parkfield was widely recorded by seismic as well as other types of instruments.
The data set used here consists of 10 digital three-component strong-motion displacement seismograms with a duration of 18 s. Six stations are located on the northeastern side of the fault and four stations are on the southwestern side (Fig. 1). In addition, there are also 33 analogue stations that recorded this earthquake. The use of digital stations has two advantages: first, they have absolute timing, and secondly, the first P-wave arrival is recorded. Though we do not use data from analogue stations for the inversions, we shall use them later for a cross-check on our model by calculating the displacements at those stations and examining how well the waveforms are matched. All the displacement records were bandpass  Borcherdt et al. 1985) are shown as red triangles, whereas the blue triangles represent available analogue stations (CGS network). The fault trace used for the inversions is shown by the thick black line; the grey lines show the surface expression of the San Andreas Fault. The 'beach ball' is connected by a thin line to the earthquake epicentre. Start time of each trace is the origin time of the earthquake. filtered between 0.16 and 1 Hz. The lower frequency limit is chosen based on the instrument capability of analogue stations. Though the lower frequency limit for the digital stations could be taken down to 0.1 Hz, we used 0.16 Hz to keep consistency between different data sets. The higher frequency limit is determined by the accuracy of the velocity structure and Green's functions (Liu et al. 2006).

K I N E M AT I C I N V E R S I O N M E T H O D
We model the fault as a rectangular plane, 40 km long (30 km to the northwest and 10 km to the southeast of the hypocentre) and 16 km wide along depth, its surface projection being shown as a straight black line on Fig. 1. The size of the fault is based on the location of the aftershocks that occurred within 24 hr of the main shock; strike and dip are taken as N140E and 87 • , respectively; hypocentral location is 35.82 • N, 120.37 • W, at a depth of 8.3 km (Thurber et al. 2006). The velocity model used for the computation of Green's functions is the 1-D structure used by Liu et al. (2006).
The method of elliptical subfault approximation has been used for kinematic inversions by Vallée & Bouchon (2004), Peyrat et al. (2010) and Di Carli et al. (2010). In this approach, each elliptical patch is defined by the seven parameters shown in Fig. 2. In each ellipse, the slip distribution has a smooth Gaussian distribution from the maximum slip amplitude at the centre to zero slip amplitude at boundaries. We require that the first ellipse contains the hypocentre. The centre of the ellipse is calculated using two parameters: h r and α h , which control the position of the ellipse relative to the hypocentre. α h is the azimuth of the centre of the ellipse about the hypocentre; h r is the distance between the hypocentre and the centre of ellipse (this distance cannot be greater than the length of the semi-major axis). For the first ellipse, the rupture time is calculated assuming a circular rupture front, starting from the hypocentre at t = 0, and propagating at the rupture speed associated with the first ellipse (This rupture speed is one of the parameters we invert for). The rupture of the second ellipse is initiated at the point where the circular rupture front, from the hypocentre, makes its first contact with the ellipse. To determine the time when the second ellipse starts rupturing, we calculate the time needed to travel from the hypocentre to the initiation point with the rupture speed of the second ellipse. This rupture speed is also a parameter inverted for. The same process is repeated for other ellipses, when used.
The complete wavefield Green's functions, including near-field terms, are computed over the fault using 500 m × 500 m cell sizes, using the AXITRA code (Cotton & Coutant 1997), which combines the reflectivity method (Kennett & Kerry 1979) with the discrete wavenumber decomposition (Bouchon 1981). The same code is used to simulate the wave propagation and compute solution seismograms.
During an inversion, the aim is to minimize a cost function. This function measures the difference between observed data (u o ) and modelled seismograms (u s ). For our study, we choose the following cost function (Spudich & Miller 1990): (1) In eq. (1), W i is the weight given to each station. A higher weight is given to stations with a high signal-to-noise ratio, with values chosen as in Liu et al. (2008). N d is the number of records and (t b , t e ) gives the beginning and end times. The search algorithm for the best parameters used in this paper is the neighbourhood algorithm (Sambridge (1999a,b), which simultaneously searches for the best values of the parameters. This algorithm obeys the following reasoning: The first step is to uniformly sample n i models inside the parameter space. The seismograms are then computed for each model, and a misfit value is found for each, using the cost function ( ). n r models with the smallest misfit are then selected from the initial n i models. According to the distribution of those n r models, a Voronoi diagram (Voronoi 1908) is built, where each model is associated with a Voronoi cell. A new set of n i models in the regions defined by the Voronoi cells is then sampled. Initially, tests were performed using artificially constructed data to obtain better insight into the method of inversion, for various cases (see Appendix). In particular, in 'Test 3', we address the question of the number of ellipses to be used, when no previous study of the earthquake exists.
During each inversion, we choose to fix two source parameters: the rise-time (τ ) and the rake, constant over the entire fault. Due to the range of frequencies used in this study, the rise-time may not be resolvable (Liu et al. 2006). Based on a previous kinematic inversion and dynamic modelling, we choose a value of 0.5 s for τ (Liu et al. 2006); (Ma et al. 2008). Previous studies also show that the rake angle does not have large fluctuations, so it is reasonable to assume the rupture to be purely right lateral (Custodio et al. 2009).

C H O I C E O F T H E P R E F E R R E D M O D E L
We performed 12 inversions using different parameters to find the source process for the 2004 Parkfield earthquake (see Supporting Information Table S1). Fig. 3 shows the final slip distributions for the 12 inversions and Table S1 gives details of the different inversions, and the resulting source processes. For 11 inversions, two ellipses were used, and one case with three ellipses is also considered. The choice of two ellipses was based on the slip distributions inferred from other studies of this earthquake, using different methods (e.g. Murray & Langbein 2006;Allmann & Shearer 2007;Ma et al. 2008). For some of the inversions, we also required the second ellipse to be connected to the first one. In that case, (h r , α h ) relate to the centre of the first ellipse, rather than to the hypocentre. We also carry out some inversions in which we vary the rise-time (τ ) and the rake to examine their influence on the solutions. In addition to that, two inversions using only analogue stations were also carried out and are reported in Sections S14-S15 of the Supporting Information S2. We discuss next the choice of our preferred model among these inversions. Fig. 3 shows that several models have similar misfits, and this leads to the problem of selection of the preferred solution. In spite of the diversity, we find that there is consistency between the models. To highlight the robust features which are independent of the choice of a priori conditions (Table S1), a simple average of all the models, inversely weighted by the misfit value, is plotted in  Information Table S1 for details), using digital stations. The misfit E is given in the top right of each fault. The red star shows the position of the hypocentre. The preferred model is highlighted by the red rectangle.   Fig. 4, we shade regions of slip which are smaller than the standard deviation calculated over the 12 inversions. Therefore, the remaining area focuses on the robust features extracted from the 12 inversions. The misfit of 0.35 for this average model (calculated assuming a constant rupture velocity of 3.0 km s −1 ), falls within the range of values obtained for the 12 inversions, showing that the features highlighted by this averaged slip distribution does not lead to a non-realistic model. On average, two slip patches are necessary to explain the data, one close to the hypocentre and a second one located between 15 and 20 km northwest of the hypocentre, with no slip in the top 5 km. This suggests that the slip distribution of the preferred model should include these properties. It is also important to note that we find a strong resemblance with previous models of this earthquake, obtained from analysis of strong motion data but using different inversion methods (Ma et al. 2008); Custodio et al. (2009). This indicates that the main features are independent of the approach used to find the slip distribution.
As we also invert for the rise-time and the rake in some of the inversions, it is interesting to discuss their influence on the obtained solutions. The rake varies between 140 and 180 • in all the inversions, with a constant value over the entire fault. This is consistent with the results of Liu et al. (2006) showing that some parts of the fault exhibits a combined right-lateral and reverse motion. The rise-time shows a higher variability, making it difficult to discuss the reliability of the results. However, it is interesting to see that Inversions 6 and 8 show a similar final slip distribution with different values for the rise time (see Tables S6 and S8 in Supporting  Information S2).
This implies that this parameter cannot be resolved, and therefore its impact on the final solution may not be significant. Of the 12 inversions, eight have significant slip at the hypocentre (Inversions 2, 3, 5, 6, 8, 9, 10, 12-see Fig. 3). Among these eight, only Inversion 5 does not have a high-slip patch northwest of the hypocentre. If we also require the preferred model to have a moment value within ±15 per cent of the CMT value (i.e. between 0.96 × 10 18 and 1.30 × 10 18 N m), then of the remaining seven models, only Inversions 2, 6, 8 and 10 satisfy this condition. We may reject Inversion 2 from the list of acceptable models, due to the fact that the rupture speed in the first ellipse exceeds the shear wave speed at the hypocentral depth by ∼30 per cent, and no previous study found supershear rupture speed for this earthquake (e.g. Fletcher et al. 2006, and others referred earlier). Finally, by generating the seismograms at the analogue stations for the three remaining models, we calculate the values of the misfit for digital stations and analogue stations together. Since Inversion 6 has the lowest global misfit (0.82, with 0.26 and 0.56 for digital and analogue records, respectively) we take it as our preferred model. Fig. 5 shows the comparison of the waveform data generated by the preferred model with the seismograms recorded at the 10 digital stations. We see that the seismograms show an excellent match in the wave shape and an excellent timing for the main pulses. In general, the amplitude is very well retrieved in the early part of the record, though we note that the later arrivals are not well fitted. It has been known for some time (Li et al. 1990); Ben-Zion (1998) that waves from the fault zone gouge exist and could have a strong effect on the recorded strong-ground motion (Fukushima et al. 2010). The deterioration of the agreement between the data and the solution seismograms in the later portions of the record at some stations may stem from the fact that these fault-zone waves are not modelled in our study. The same reason could also explain the mismatch in amplitude in the early part of some stations (e.g. MFU), situated very close to the fault. Fig. 6 shows the comparison between waveform data and the preferred model solution seismograms at the 33 analogue stations. In the forward direction of rupture propagation (i.e. northwest from the hypocentre), a good agreement is observed only for the main pulses. In the backward direction, waveforms at some stations are not well matched. To have a more quantitative view of this mismatch, we plot the misfit values at all the analogue stations separately for our preferred model (Fig. 7a). The red dashed line represents the mean of the individual misfit value, so we can consider that stations above this line exhibit a significant mismatch. We identified three locations of higher mismatch (Fig. 7b), one at the extreme northwestern end of the fault (Station COAL), one near the region of the high-slip patch (Stations FZ12, FZ15, VC1W, VC2W, VC3W and VC5W) and one in the southeastern end of the fault (Stations FZ1, C1E, CH2W and CH3W). It is interesting to see that most of these stations are very close to the fault, which is a propitious area to be influenced by fault zone trapped waves (e.g. Li et al. 2004), and could explain the high amplitude and the highest mismatch between observed and calculated seismograms at those stations. For stations close to the high-slip patch, we believe that this area which released a significant amount of the total energy during the earthquake, could have resulted in enhanced site effects, which could explain the higher mismatch for stations at larger distances to the fault.

D I S C U S S I O N O F T H E P R E F E R R E D M O D E L
Previous models using seismic data (e.g. Allmann & Shearer 2007;Custodio et al. 2009) as well as our study have some slip at the hypocentre. A study of the Parkfield earthquake using geodetic data (InSAR and GPS) by Johanson et al. (2006), has argued that slip at hypocentre may be due to rapid after-slip. Hypocentral slip may only be needed to explain the large amplitude observed at analogue stations FZ1, C1E and C2W (see Fig. 6), as also pointed out earlier (e.g. Shakal et al. 2006). Also, some of our inversions either have no slip or, no significant slip at hypocentre (Inversions 1, 4, 5, 7 and 11). Among these, Inversions 1, 5 and 9 explain the signal observed at analogue stations even better than our preferred model. To test if some slip at the hypocentre is a source signal or an artefact induced by fault zone waves affecting stations at the southeastern extremity of the fault, we run an inversion using all digital and analogue stations, except FZ1, C1E, CH2W and CH3W, which have significantly larger misfit, and the same set-up as in Inversion 6. The resulting slip distribution (Fig. 7c), shows that even without these high-misfit stations, which also have very high amplitudes, some slip at the hypocentre is needed to explain the data.   Fig. 5 but for the analogue stations. A comparison of the observed and the solution seismograms for all inversions carried out in this study can be found in Sections S1-S12 of the Supporting Information S2.

T H E S O U RC E P RO C E S S
The slip distribution associated with our preferred model (Inversion 6) is shown in Fig. 8(a). The average displacement over the whole fault is of ∼0.07 m, the maximum displacement within the second ellipse being 0.91 m. This high-slip patch is located ∼17 km from the hypocentre in the northwesterly direction, and explains the high amplitude for stations located in the northwestern end of the fault. This is one of the most robust features of the 2004 Parkfield earthquake, and has also been found in other seismic (Allmann & Shearer 2007); (Custodio et al. 2009) and geodetic (Johanson et al. 2006); (Johnson et al. 2006); (Murray & Langbein 2006) studies. It is very important to note that this area of high slip coincides spatially with the hypocentre of the 1966 M w 6.0 Parkfield earthquake (white star in Fig. 9). Harris & Segall (1987) and Malin et al. (1989) had earlier identified this region of high slip as a locked zone or asperity. We can therefore interpret this high-slip patch as a permanent asperity, which has been activated during the 2004 Parkfield earthquake. Custodio & Archuleta (2007) suggested the presence of persistent asperities in the Parkfield region that can rupture together or individually during earthquakes. The unconstrained inversion of the 1966 earthquake gives a high-slip patch located in the same area as the high-slip patch observed for the 2004 Parkfield earthquake. Though the limited data for the 1966 earthquake are inadequate to resolve details, Custodio & Archuleta (2007)'s study gives an idea of where slip could have occurred, which suggests that some asperities ruptured both in 1966 and 2004, whereas others broke only in one of the events.
Our preferred model shows that the rupture of the 2004 Parkfield earthquake propagates to the northwest at an average speed of ∼2.7 km s −1 which is about 80 per cent of the local shear wave speed. Our propagation time agrees well for the first 3 s of the total ∼8 s rupture process, with Fletcher et al. (2006), and with Allmann & Shearer (2007) for the entire process. The hypocentral ellipse has a rupture speed of ∼2.2 km s −1 . The second ellipse, located between 15 and 25 km northwest of the hypocentre, has a higher   . Larger background seismicity (3 < M w ≤ 5) prior to the 2004 Parkfield earthquake, from 1984 January 3 to the day before this earthquake, plotted over our preferred slip distribution. The catalogue used is from Thurber et al. (2006). rupture speed of ∼3.1 km s −1 , which is within the range of velocities found in other studies (Borcherdt et al. 2006); (Liu et al. 2006); (Ma et al. 2008). After the first ellipse completes rupturing in 3 s, a short pause of ∼1 s is observed, before the rupture starts to break the second patch, this taking ∼5 s (Fig. 8b).
Most inversions of the 2004 Parkfield earthquake find the slip to be continuous over the fault (e.g. Ma et al. 2008;Custodio et al. 2009). However, as emphasized by Vallée & Bouchon (2004), the aim of the method used here is to focus on finding the major slip areas, which explain a large part of the waveform and are also the best resolvable features of the co-seismic slip. This is why we do not observe any transitional slip between patches while other models do.

R E L AT I O N B E T W E E N H I G H -S L I P PAT C H A N D S E I S M I C I T Y P R I O R T H E 2 0 0 4 PA R K F I E L D E A RT H Q UA K E
We used the catalogue of Thurber et al. (2006) from 1984 January 3 to 2004 September 28, to examine the relationship between previous seismicity and the northwest asperity. Fig. 9 shows the locations of the main earthquakes (M w > 3). This figure shows that earthquakes, especially larger ones, surround the region of highest slip. This observation is in good agreement with the hypothesis of an asperity in this area, which produces stress accumulation at its edges, that is partly released during smaller earthquakes. This correlation between prior seismicity and asperities has been observed in previous studies (e.g. Hsu et al. 1985). Ben-Zion & Rice (1993) showed that this 'Parkfield asperity' has a major influence on prediction attempts, and must play a role in the slip budget of the Parkfield segment (the part of the San Andreas Fault which experienced the five last M w 6.0 Parkfield earthquakes). Toké & Arrowsmith (2006) show that the northwestern part of the Parkfield segment has a slip deficit close to zero, since this is adjacent to the creeping section of the San Andreas Fault. However, the southeastern part of the Parkfield segment (the Cholame segment) has a slip deficit of ∼5 m since 1857. The Parkfield segment also has a slip deficit but slightly lower than the Cholame segment (∼3.5 m), due to the release of energy by the recurrent Parkfield earthquakes. So, a hypothetical earthquake, which breaks the Cholame segment would have a higher magnitude if it also breaks the Parkfield segment at the same time.

R E L AT I O N B E T W E E N M A I N S H O C K S L I P A N D A F T E R S H O C K S
We used the relocated aftershock sequence of the 2004 Parkfield earthquake of Thurber et al. (2006) from 2004 September 28 to 2005 June 30, to compare their locations with our slip distribution. One salient feature is the horizontal level delineated by small aftershocks around 5-6 km depth (Fig. 10a). Waldhauser et al. (2004) suggest that this level may represent a change from creeping (above), to locked (below) in the interseismic period. Though the depth of the ellipses in our inversion was not fixed to be below this level, the seismic data require the ellipses to lie below it. The larger aftershocks (M w > 3, Fig. 10b) mainly lie where no slip has been observed, though there is a cluster ∼19 km northwest of the hypocentre, located at the edge of the northwestern high-slip patch, where there is a sharp change from high to low slip. This behaviour of aftershocks has been observed previously (see Das & Henry 2003, for other examples) and provides independent support for the reliability of the position of this high-slip patch.

CONCLUSION
We performed a kinematic inversion of the 2004 September 28 Parkfield, California, earthquake using a recently developed method, which defines the slip distribution as an aggregate of ellipses. The method was tested using artificial data (Appendix). We fit well the early portions of most seismograms and suggest that modelling of fault zone waves is required to explain some of the later waveforms. Our preferred slip distribution (Inversion 6) is composed of two distinct ellipses, and shows no slip in the top 5 km. A horizontal lineation of small aftershocks at this depth of 5 km has been suggested as marking a sharp transition in the interseismic slip rate. The highest slip occurs in a region located between 15 and 20 km from the hypocentre, in the northwestern direction. This patch can be interpreted as a permanent asperity, which is activated during large earthquakes. The presence of this asperity has important implications for seismic hazard assessment since it may be a characteristic feature of the Parkfield earthquakes.

A C K N O W L E D G M E N T S
Research by CT was supported by QUEST, funded under the EC Marie Curie Initial Training Network Grant Agreement No. 238007. We thank the Center for Engineering Strong Motion Data (CESMD) and the California Geological Survey (CGS) for the use of their data. We also thank Javier Ruiz who provided the k −2 -type slip distributions. Finally, we would like to thank two anonymous reviewers for very thoughtful comments which helped to improve and clarify the paper.

A P P E N D I X A : T E S T O F T H E I N V E R S I O N M E T H O D , U S I N G A RT I F I C I A L DATA
We construct artificial data, using the same configuration as presented in the main paper, for the digital stations. The rupture is initiated at the same location than the 2004 Parkfield earthquake hypocentre, rupturing with a prescribed rupture velocity of 3 km s −1 . Three tests, using three different slip distributions, were carried out. The rupture velocity and the moment of the artificial earthquake is retrieved almost perfectly for each of the tests.
For this earthquake, we have used results from previous studies to choose the number of ellipses to use. However, if we want to apply this method blindly, some tests on how to choose the number of ellipses need to be carried out. This will be considered in 'Test 3'.

A1 Test 1
Artificial displacement seismograms were generated at digital station locations for a model with two distinct ellipses. This test was carried out to see how the programmes performed in a very simple case. Fig. A1 shows the artificial model (a) and the inverted model (b). The artificial model is almost perfectly retrieved. The major difference is seen in the amplitude of the ellipses, with the right side ellipse being 7 per cent higher than the initial model and the left ellipse being 13 per cent lower. The fit to the data is very good ( = 0.03), and there is no difference between the artificial data and the solution seismograms.

A2 Test 2
The artificial seismograms were generated for a k −2 -type slip distribution (Fig. A2a), following the method described in Ruiz et al. (2007). In the inversion, we allow rupture and slip to occur only on a single asperity located near the hypocentre. The result is shown  in Fig. A2(b). The blue line on Fig. A2(a) shows the contour of the ellipse inferred from the inversion. We can see that the main asperity of the artificial model and the ellipse of the inversion are spatially close, though the ellipse is bigger than the main asperity. This is because the inversion is constrained to include all the slip of the artificial model within a single patch, while it aims to obtain the proper moment and fit the seismograms. This test shows that despite the heterogeneous slip distribution of the artificial model, the inversion method retrieves the main asperity. The seismograms show a good fit ( = 0.16), with small mismatch at some stations (Fig. A3). This can be explained by the small amplitude asperities, which we do not try to extract.

A3 Test 3
This test is performed using a k −2 -type, double asperity, slip distribution. One goal here is to gain insight into how to choose the  number of ellipses to be used in the inversion, when no other studies of an earthquake exists. It was suggested by Vallée & Bouchon (2004) that the number of ellipses should be increased when part of the signal, believed to be caused by source, is not fitted, a procedure that they applied to the 1995 Jalisco (Mexico) earthquake. To explicitly illustrate this, we first invert the artificial model using only one ellipse. The slip distribution obtained from the inversion and the comparison of the waveforms are shown in Figs A4 and A5, respectively. We can see that in the case where two major asperities are present, the inversion tends to produce a final slip distribution with an ellipse located at the centroid of the distribution. However, examination of the displacement waveforms shows non-negligible mismatch between the solution seismograms and the artificial ones ( = 0.19). In a real problem, this would have alerted us to the fact that a higher number of ellipses are needed. We next used two ellipses to see if the fit is improved. The results are shown in Fig. A6 (comparison of final slip distribution). The two-ellipse inversion correctly retrieves the two major asperities. The inverted model also shows a small patch with high amplitude between the two ellipses. This does coincides with a small peak in amplitude in the artificial model, but is more likely an artefact of the inversion due to the proximity of the two ellipses. The waveform fits are also significantly improved ( = 0.05).

S U P P O RT I N G I N F O R M AT I O N
Additional Supporting Information may be found in the online version of this article: Table S1. Details of the different inversions carried out in this study and the resulting source processes. Supplementary Information S2 (Section S1-S15). Comparison between the observed and solution seismograms for the two different sets of stations (digital and analogue). Labels at the top of each page give the name of the Inversion as described in the main text. A summary of the parameters describing each Inversion is also attached with this file.
Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.