Three dimensional velocity structure and relocated aftershocks for the 1985 Constantine, Algeria (Ms = 5.9) earthquake

Local earthquake travel-time data were inverted to obtain a three dimensional tomographic image of the region centered on the 1985 Constantine earthquake. The resulting velocity model was then used to relocate the events. The tomographic data set consisted of P and S waves travel-times from 653 carefully selected aftershocks of this moderate size earthquake, recorded at 10 temporary stations. A three-dimensional P-wave velocity image to a depth of 12 km was obtained by Thurber's method. At shallower depth, the velocity contrasts reflected the differences in tectonic units. Velocities lower than 4 km/s corresponded to recent deposits, velocities higher than 5 kmls to the Constantine Neritic and the Tellian nappes. The relocation of the aftershocks indicates that most of the seismicity occured where the velocity exceeded 5.5 M s . The aftershock distribution accurately defined the three segments involved in the main shock and led to a better understanding of the rupture process.


Introduction
On October 27, 1985, a moderate size earthquake occurred near Constantine (Algeria) (36.46"N, 6.76"E, depth = 10 km, M, = 5.9, from NEIC).It was the strongest event recorded since 1947 in the eastern part of the Tellian Atlas.The focal mechanism of the main shock obtained by P and SH wave modelling shows a left lateral strike-slip motion on 217" striking fault (Aki's convention), dipping Mailing address: Dr. Catherine Dorbath, Ecole et Obsewatoire des Sciences de la Terre, 5 rue RenC Descartes, 67084 Strasbourg Cedex, France; e-rnail: cath@mapu.u-strasbg.fr at 84"NW (Deschamps et al., 1991), which is in agreement with the solution determined from direct reading of first motions as well as fault traces observed on the surface (Bounif et al., 1987).An examination of the quaternary faults in the region indicated that the surface rupture may correspond to the northeastern continuation of the Ai'n Smara fault (fig.1).
Three days after the main shock a network of 10 MEQ 800 stations was deployed in the epicentral area in order to record the aftershock sequence.These stations operated for more than three weeks, from October 30 to November 23.Some of them were moved during the experiment, so that 15 distinct sites were occupied (fig.1).The geometry of the array was designed for hypocentral determinations, with a good azimuthal coverage and short epicentral distances relative to the aftershock sequence.
During this field experiment about 1500 events were recorder.A first study concerning 150 aftershocks has been published by Bounif et al. (1987).A report on the complete data set has been presented by Bounif (1990) and was summarized in Deschamps et al. (1991).The spatial distribution of the epicenters defines a 30 km long and 2 km wide, NE-SW striking zone.The foci are distributed on nearly vertical fault planes.Examination of the aftershock dis-tribution shows evidence of three ruptured segments, characterized by differences in their general orientation and the depth of the corresponding events.As no crustal velocity model was available for this region, simple half-space or two-layer models were assumed for hypocentral locations; this leads to large uncertainties on final positions, particularly on depths.
The geology of the region (fig.2) is domi-affect middle to upper Quaternay deposits, nated by nappe-structures which developed which are aligned in a NE-SW direction, as is mainly during the Late Cretaceous and Paleo-the main surface faulting associated with the gene and resulting from Alpine orogeny.The earthquake.These deposits lie unconformably surface ruptures observed after the earthquake on a substratum of mainly Eocene flysch Fig. 2. Simplified structural sketch of the El Ana region, after Coiffait et al. (1997).1) Mio-Pliocene and Quaternary deposits; 2) Numidian nappe; 3) Tellian nappes; 4) Constantine Neritic nappe.A = main faults; B = thrust faults; C = thrust sheets.The tomography grid is superimposed in grey, toghether with the position of some stations.(Coiffait et al., 1997;Vila, 1980), which is affected by several faults trending in the same direction (Bounif et al., 1987).
In order to obtain an estimation of the velocities in the Constantine earthquake area, we performed an inversion of aftershock traveltime data.The velocity model obtained was then used to relocate the hypocenters with increased precision, leading to an improved image of the aftershock distribution and a better understanding of the rupture process.

Data and inversion
Among the 1500 events in the Constantine sequence recorded by the temporary network, 1200 were located by Bounif (1990) using only the P-wave arrival times.To perform the hypocentral determination with the HYP~INVERSE program (Klein, 1978).Bounif used a simple velocity model consisting of a 12 km thick layer (V, = 5 M S ) above a half space (V, = = 6.5 M S ) , which minimizes the RMS residuals and the standard errors.
As we had only vertical seismometers, S wave were often difficult to read.However, for the purpose of this study and in order to increase the control on the hypocentral depths, we included reliable S readings in the data set.The reading errors are estimated to be less than 0.1 and 0.5 S for P and S waves respectively.Hypocenters were relocated with the HYPOIN-VERSE program, using the same initial velocity model, by scanning the trial depth.We kept only the very best located aftershocks and we rejected hypocentral locations with RMS greater than 0.2 S, horizontal error greater than 1.5 km or vertical error greater than 3.0 km.So we selected 653 hypocentres with mean RMS = 0.07 S, ERH = 0.65 km and ERZ = 1.17 km.They provided 5288 arrival times, 90% of which correspond to P waves.
The inversion program we used is the latest version of the program adapted by Eberhart-Phillips (1990) from the method developed by Thurber (1983).The method has been applied to numerous earthquake source regions and is fully described by Thurber (1993) and Eberhart-Phillips (1993).This method combines pa-rameter separation (Pavlis and Booker, 1980) and damped least square inversion.The iterative process inverts P and S-wave arrival-time data for a 3D velocity structure and hypocentral parameters.The parametrization of the region under study is achieved by assigning velocity values at fixed points on a 3D grid.An approximate ray tracing algorithm (ART, Thurber, 1983) is used to calculate travel-time between the source and the station.
For the initial P-wave velocities, we chose a smoothed version of the velocity model used for hypocentral locations; S-wave velocities were set by using a V,/V, ratio of 1.73.The velocity grid was rotated in order to follow the mean direction of the aftershocks distribution which is N34"E; the horizontal grid spacing is 6 km, the vertical grid spacing is 4 km from the surface down to the maximum depth of hypocenters.
To select the optimal value of the damping factor, we constructed <<trade-off>> curves of resolution against solution error.The optimal value substantiallv increases the resolution without introducing large velocity errors.We performed ten inversions with damping factor values ranging from 2 to 200, and finally selected the value of 10 which allows a good resolution with errors lower than 0.1 MS for P-wave velocities.
The stability of the solutions was been tested by standard methods: variation of the horizontal and vertical spacing of the grid points, rotation of the grid, translation of the coordinates of the origin of the grid, variation of initial velocities.In any case, the velocity changes in the P velocity models obtained are lower than the errors and therefore not significant.The images of the velocity structures depend on the source-station distribution and not on the parametrization of the space.
The initial P data variance was 0.03 1 s2.After 5 iterations, the reduction of the initial data variance was 60%, to 0.013 s2, which is close to the value expected from the magnitude of reading errors.
The reliability of the inversion is essentially indicated by the value of the diagonal elements of the resolution matrix (fig.3), while the value of off-diagonal elements indicates the coupling between neighbouring nodes.The largest values of the diagonal elements of the resolution matrix reach 0.6.At the surface, the highest values are clearly correlated with the position of the stations.The best resolved depths are 4 and 8 km, where the ray paths are numerous.Below 12 km, the resolution becomes too low because the nodes are not well sampled by rays.For the final P model, the standard deviation within each layer is about 0.25 M S , whereas the average standard error is about 0.06 M S , therefore the P velocity variations are significant with respect to the uncertainties.The S velocity model has a very low resolution, due to the small quantity of S data, and it will not be considered here; however the S travel-time data were used in the relocation process during the inversion because they constrain hypocentral depths.

. Velocity model
The 3D velocity model is presented in map view in fig. 4 for the four layers, together with the location of the aftershocks at corresponding depths.At the surface, we observe some clear coincidences with the geological units.The lowest velocities, down to 3.5 krnls, correspond to Mio-Pliocene and Quaternary deposits which have been slightly or not tectonised (around stations KES or ABY, fig.2).The highest velocities, up to 5.5 M S , correspond to the Constantine Neritic nappe (station BCD) and the Tellian nappes (stations CNE, GRB), which are Jurassic to Paleogene in age.The Miocene Numedian nappe corresponds to intermediate velocities.The same general velocity contrasts but weaker (from 4.5 to 5.5 M S ) are observed in the second layer, at a 4 km depth.From this extension in depth of the velocity contrasts we infer that the nappes are a few kilometers thick.In the third layer, the velocity patterns are more scattered.The highest velocities, more than 6.5 M S , are observed to the south east of the place where the activity is maximum.The velocity contrasts in the deepest layer are not very strong, which may result from the poor resolution.

Earthquake locations
The relocated seisrnicity (fig. 1) shows the same global pattern as in the previous studies but the aftershock distribution is better defined and concentrates in a much narrower zone.We present in fig. 5 on a map view a comparison between the relocations made with our final velocity model (circles) and the initial epicentral locations (extremity of the relocation vectors).Overall, the shifts are modest but reduce the scattering.In particular, the segmentation in three segments characterised by specific trend is very clear.The most important shifts are observed at the junctions between the different segments.In the northern segment, the differences are very small, about 500 m, whereas they are about 2 km in the southern segment.
Cross sections through the aftershocks along profiles orthogonal to the three segments show that larger shifts are observed on the hypocentral depths, as expected (fig.5).The seismicity is mainly confined between 6 and 13 km on the northern and central segments, between 8 and 10 km on the southern one.Very few earthquakes locate in the upper part of the crust, none below 15 km.Moreover, the cross sections exhibit a decrease in the dip of the fault plane from north to south.The fault plane is vertical on the northern segment, dipping 80" to the northwest on the central segment and 70" to the northwest on the southern segment.This modest dip to the northwest is in good agreement with the dip of the focal mechanisms obtained for the main shock by Deschamps et al. (1991) (84"), or with the dip of the CMT solution (71").

Relation between the velocity structure and the relocated hypocentres
To select our initial data set, we imposed very restrictive criteria on the hypocentral locations in order to reject poorly constrained solutions.Nevertheless, the uncertainty on the depth drastically increases for the shallower events when no seismic station is situated just above the hypocentre.Thus, the criterion on the vertical error leads to the rejection mainly of the shallower events.We attempted to eliminate this bias by selecting a more comprehensive data set including shallow locations with a higher ERZ, the criteria on RMS and ERH remaining unchanged.Then we relocated all these events with our final 3D velocity model.
In fig.6 we show five velocity depth sections and the relocated hypocentres, oriented NW-SE, perpendicular to the global trend of the aftershocks (N34"E).As the seismicity is distributed along three segments with slightly different trends, the sections are not orthogonal to the segments and the fault planes are not as precisely defined as on the fig. 5.In fig.7 we present the y = 36 km section which cuts the seismicity lengthwise.The most interesting feature evident on these sections is the relation between the seismic activity and the velocity structures.Most hypocentres locate where the velocity exceeds 5.5 M s .The minimal depth of these clusters of events increases as the 5.5 km/s velocity isoline becomes deeper, from x = 24 to x = 42 km on fig. 7.Only at the junction between the southern and the central segments do we observe a cluster of events up to the 5 km/s isoline (X = 24 km).
Along the southern segment (fig.6, X = 18 km, and fig.7), we observe the lowest velocities at depth (5 km/s down to 7.5 km, 5.5 km/s down to 9 km).This segment is about 8 km a long and the seismicity is mainly concentrated between 9 and l 1 km in depth.The central segment is about 14 km long and the mean depth of the aftershocks increases toward the northeast (fig.7).The main activity is observed at its northern extremity (X = 30 to 36 km) between 8 and 12 km in depth, underneath the surface rupture zone.We can notice that the ruptures are located where the shallow velocity is the highest along the central segment.
The cluster located at its junction with the southern segment (fig.6, X = 24 km, and fig.7) is the most active with hypocentral depths between 6 and 10 km.It is situated below the lowest shallow velocities (4.5 km/s down to 2 km).This cluster was not seen on the data subset studied by Bounif et al. (1987), which included only the earlier aftershocks up to November 11.It fills the horizontal shift that they observed between the southern and central segments.The northern segment (fig.6, X = 36 to 42 km, on fig.7) is about 8 km long.It presents the weakest activity spread over the largest depth distance with a maximum between 9 and 12 km.This shallowing of the hypocentres depth corresponds to a steep shallowing of the 5.5 km/s isoline to the north of this segment.
The comparison between velocities and depths of the hypocenters shows that the structures with velocity lower than 5.5 km/s are almost not brittle.Only the junction between the central and the southern segments does not exactly follow the same relation.Moreover, the brittle part of the crust is thin, about 5 km.The mean uncertainty on the relocated hypocentres is weak, 0.62 km horizontally and 1.35 km on the depth, so the relation between aftershock locations and velocity structures is reliable.A similar relation between activity and velocity has been found in other parts of the world where local earthquake tomography has been performed, particularly in California (Eberhart-Phillips, 1990;Dorbath et al., 1996).In spite of the different initial velocity models used in these papers, the seismicity generally occurs at depths where P-wave velocity exceeds 5.5 M s .In those studies, the seismogenic layer is associated with the basement layer and it is consistent to make the same association in the Constantine region.
On a cross section parallel to the aftershocks activity, Bounif et al. (1987) observed that the central and the northern segments corresponded to a circular distribution of foci between the surface and 15 km depth around a quiet zone.They interpreted these two patches as asperities broken during the main shock, and then assumed that the principal segment of surface faulting (3.8 km) is only the visible part of the fault plane.Therefore they consider the entire fault plane to be 20 km long and 10 km wide with a 10 cm dislocation, in which case the moment released would be M, = 6 .1017 N .m, which is consistent with the value of 6.2 .1017 N .m estimated from the centroid solution.Such an interpretation is difficult to deduce from Our results on the distribution of the aftershocks (fig.7).We propose to attribute only the central segment, from the surface to the deepest aftershocks, to the main shock.With this hypothesis, the fault plane is then 14 km long and 12 km wide and dips 80" to the West.
Taking a seismic moment of 6.2 .10" N .m, we find a 12 cm dislocation, which is reasonable for a magnitude M, = 5.9 earthquake and consistent with surface observations.For example, the dislocation was up to 18 cm during the 1986 Kalamata (Greece) earthquake (M, = 5.8, Lyon-Caen et al., 1988), and up to 20 cm during the 1989 Chenoua (Algeria) earthquake = 6.0,Bounif et al., 1998, submitted).The surface faulting is seen at the sur-face in the part of the fault plane where the surface velocity is the highest.In this case, the low seismicity in the upper part of the fault plane (down to about 7 km) is not related to a broken asperity but to the mechanical properties of a layer with a velocity lower than 5.5 kmls.This would imply that the aftershocks are related to a very localized brittle layer, perhaps at the base of the mainshock rupture zone, and not to the mainshock rupture itself.
The northern and southern segments, which strike in slightly different directions, are fault planes activated by the main shock.The cluster observed at the junction between the central and southern segment could then be related to stress concentration at the southern end of the ruptured zone.

Conclusions
Seismic tomography of the Constantine earthquake source region reveals 3D P-wave velocity to depths of 12 km.The velocity contrasts in the upper layers can be roughly related to the different tectonic formations.The lowest velocities (down to 3.5 M S ) correspond to Mio-Pliocene and Quaternary deposits, the highest velocities (up to 5.5 M S ) to the Constantine Neritic and the Tellian Jurassic to Paleogene nappes.The extension in depth of the velocity contrasts infers that the nappes are a few kilometers thick.
A comparison of initial earthquake locations with those resulting from the tornographic inversion shows that the latter better define possible fault planes and a segmentation of the aftershocks into three ruptured segments.Despite modest shifts relative to initial hypocentres, the relocated hypocentres are less scattered and their concentration in depth is strengthened.We observe that the relocated seismicity occurs at depths where the velocity exceeds 5.5 M S , and assign the top of the basement to this isoline.Moreover.the thickness of the seismogenic layer does not exceed 5 km.

Fig. 1 .
Fig. 1.Distribution of the stations and events used for the local earthquake tomography.The shaded circles show the aftershocks of the 1985 Constantine earthquake relocated after the inversion, the darker the deeper.The triangles indicate the location of the stations with their code names.The grid spacing used in the tomographic study is presented in grey.Black solid lines represent surface ruptures generated by the main shock, dashed lines the A'in Smara fault.The insert shows the location of Constantine in Mediterranean region.

Fig. 3 .
Fig. 3. Diagonal elements of resolution matrix for the final model.Resolution reflects station spacing and dishbution of seismicity (see fig. 1 for location of earthquake and stations relatively to the grid).

Fig. 4 .Fig. 5 .
Fig. 4. Three dimensional P-wave velocity in the four layers of the model.The gridded results have been smoothed for the purpose of display.The location of the stations is reported on the first layer.A unique grey scale is used.

OFig. 7 .
Fig. 6.Cross sections through the relocated hypocentres (including shallow events, see the text) and the 3D velocity mode1 for grid lines x = 18 to 42 km (see fig. 1).The width of the earthquake sections is 6 km.