# Shallow afterslip following the 2003 May 21, Mw = 6.9 Boumerdes earthquake, Algeria

A. Mahsas, K. Lammali, K. Yelles, E. Calais, A. M. Freed, P. Briole

► To cite this version:

A. Mahsas, K. Lammali, K. Yelles, E. Calais, A. M. Freed, et al. Shallow afterslip following the 2003 May 21, Mw = 6.9 Boumerdes earthquake, Algeria. Geophysical Journal International, 2008, 172, pp.155-166. 10.1111/j.1365-246X.2007.03594.x. insu-03603738

HAL Id: insu-03603738

https://insu.hal.science/insu-03603738

Submitted on 10 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Shallow afterslip following the 2003 May 21, M<sub>w</sub> = 6.9 Boumerdes earthquake, Algeria

A. Mahsas,<fnref n="1" /> K. Lammali,<fnref n="1" /> K. Yelles,<fnref n="1" /> E. Calais,<fnref n="2" /> A. M. Freed<fnref n="2" /> and P. Briole<fnref n="3" />

Accepted 2007 August 24. Received 2007 August 6; in original form 2007 March 26

## SUMMARY

We investigated post-seismic deformation following the 2003 May 21, M<sub>w</sub> = 6.9 Boumerdes, Algeria, earthquake using surface displacements from six continuous Global Positioning System sites that operated in the epicentral area for 2.5 yr following the event. We find up to 4 cm of cumulative horizontal displacement during that time period, with a time-dependence well fit by a logarithmic decay. Post-seismic deformation appears to continue at all sites after the 2.5-yr observation period, with rates on the order of 1 cm yr<sup>-1</sup> or less. The data is consistent with shallow afterslip (0–5 km) and shows no evidence for afterslip downdip of the coseismic rupture. The data is poorly explained by viscoelastic relaxation in the lower crust or upper mantle, or by poroelastic rebound. The concentration of afterslip adjacent to and updip of the coseismic rupture, at least in the western half of the fault, suggests that afterslip is driven by coseismic stresses. The correlation between the depth of afterslip and that of the sedimentary wedge along the Algerian margin, while coseismic slip occurs in deeper basement rocks, suggests (1) that post-seismic deformation may also involve folding and (2) that spatial variations in frictional properties along the fault correlate with the type of rocks involved.

**Key words:** Satellite geodesy; Seismic cycle; Transient deformation; Africa.

## 1 INTRODUCTION

Despite the well-established importance of post-seismic deformation during the earthquake cycle (e.g. Thatcher 1974, 1983), debate persists about the driving mechanisms (i.e. afterslip versus viscoelastic relaxation versus poroelastic rebound), where the deformation occurs (lower crust or upper mantle), and what are the relevant constitutive parameters of the candidate processes (e.g. Newtonian versus power-law viscosity). For instance, an early phase of rapidly decaying deformation following the 1992, M<sub>w</sub> = 7.3 Landers earthquake has been explained by afterslip below 10 km (Shen et al. 1994; Savage & Svarc 1997). Some of the early deformation following that event, especially the observed pattern of uplift and subsidence near the fault, suggests poroelastic adjustments to the earthquake stress changes (Peltzer et al. 1996). Finally, longer-term post-seismic deformation following that same event has been interpreted as the result of viscoelastic relaxation in the lower crust and upper mantle, possibly involving non-linear mantle rheology, either biviscous (Pollitz et al. 2001) or stress-dependent (Freed & Bürgmann 2004). Multiple mechanisms may also be acting over different spatial and temporal scales. Freed et al. (2006) showed that post-seismic deformation following the 2002, M<sub>w</sub> = 7.9 Denali, Alasaka, earthquake resulted from the combination of viscoelastic flow in the upper mantle and possibly the lower crust, shallow afterslip in the upper crust and poroelastic rebound in the immediate vicinity of the rupture.

The difficulty for identifying post-seismic processes stems, in part, from the lack of sufficient observational data sets to constrain models. Here we report on post-seismic deformation following the 2003 May 21, M<sub>w</sub> = 6.9 Boumerdes earthquake, Algeria (Fig. 1; Ayadi et al. 2003; Yelles et al. 2004). Within 3 weeks of the event, we deployed six semi-permanent Global Positioning System (GPS) stations in the epicentral area of the earthquake and operated them continuously for 2.5 yr to investigate the mechanism of post-seismic deformation that was likely to follow the event. In this paper, we describe the post-seismic signal recorded at the GPS sites and show that, over this 2.5-yr time period, post-seismic deformation is best modelled by shallow afterslip (0–5 km). Some of this shallow slip may contribute to the development of fault-propagation folds in the sedimentary wedge at the toe of the Algerian margin.

## 2 TECTONIC SETTING

On 2003 May 21, a M<sub>w</sub> = 6.9 earthquake struck northern Algeria, about 50 km east of its capital, Algiers, with about 2400 casualties and 10 000 injured, leaving 200 000 homeless, and causing extensive damage (intensity X) in the Algiers-Dellys area (Fig. 1; Ayadi et al. 2003; Yelles et al. 2003). This earthquake is among the largest well-monitored events to occur in the Western Mediterranean over the past 25 yr. Its magnitude is comparable to the Campania (Italy)

<footnote n="1">¹Centre de Recherche en Astronomie, Astrophysique et Géophysique, Bouzaréah, Algeria</footnote>

<footnote n="2">²Purdue University, Department of Earth and Atmospheric Sciences, West Lafayette, IN 47907, USA. E-mail: ecalais@purdue.edu</footnote>

<footnote n="3">³Institut de Physique du Globe de Paris, 4 Place Jussieu, 75005 Paris, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Active tectonic framework of the western Mediterranean basin. The 2003 May 21, Boumerdes earthquake and the 1980 October 10 El Asnam events are shown in red. Focal mechanisms are shown in black for other $M_w > 6.0$ earthquakes (Harvard CMT database). Black circles show magnitude greater than 4.5 earthquakes (NEIC database). White arrows at the bottom of the figure show model velocities for the Nubia plate with respect to Eurasia (Calais *et al.* 2003).</figcaption></figure>

earthquake of 1980 November and only slightly less than that of the El Asnam earthquake that struck Algeria on 1980 October 10 (Yielding *et al.* 1989, Fig. 1).

Regional seismicity and GPS measurements show that most of the ~5 mm yr⁻¹ oblique convergence between the Nubian and European plates in the western Mediterranean is accommodated along the coastal regions of northern Africa (Nocquet & Calais 2004). Onshore active structures define a series of ENE–WSW trending folds and reverse faults affecting Neogene and Quaternary basins, with a right-stepping *en échelon* pattern and connected by NW–SE to E–W trending strike-slip faults (Meghraoui 1991; Meghraoui *et al.* 1996). Recent offshore investigations show that the margin is also tectonically active, with fault-propagation folds that develop above flat-ramp faults and affect the most recent sediments (Déverchère *et al.* 2005). The Boumerdes earthquake was a shallow thrust event on a ENE–WSW trending, south-dipping, reverse fault that project to the surface about 10 km offshore the Algerian coast. It is likely to have ruptured a segment of this active offshore thrust system (Yelles *et al.* 2004).

Estimates of source kinematics from teleseismic waveforms, GPS measurements, and coastal uplift data, show that most of the coseismic slip (up to 2.5 m) occurred in the western half of the rupture, between 5 and 10 km depth (Delouis *et al.* 2004; Yelles *et al.* 2004; Semanne *et al.* 2005). The source kinematics is less well determined in the eastern half of the rupture, where Delouis *et al.* (2004) find slip at shallow depth (0–5 km) while Semanne *et al.* (2005) find that slip concentrates between depths of 10 and 15 km.

The Boumerdes earthquake is an interesting analogue to reverse faulting events along the San Andreas fault system in the Los Angeles basin, with similar transpressional tectonic settings and shortening rates. Lin & Stein (1989) noted the similarity of the 1980 $M_s = 7.3$ El Asnam earthquake in Algeria with the 1983 $M_s = 6.5$ Coalinga, 1984 $M_s = 6.1$ Kettleman Hills and 1987 $M_l = 5.9$ Whittier Narrows earthquakes in California in that they occurred beneath actively growing folds. The same holds for the more recent $M_w = 6.6$ Northridge earthquake, a typical blind thrust event in the northern Los Angeles basin (Hauksson *et al.* 1995). Because of its relatively large moment, the accessibility to geophysical measurements, and our ability to immediately respond to the event, the Boumerdes earthquake is a particularly valuable target to investigate post-seismic deformation in a transpressional regime.

## 3 GPS DATA AND PROCESSING

Five days after the Boumerdes event, we initiated the deployment of GPS instruments in the epicentral area (hanging-wall block) in order to capture the post-seismic signal. We used dual frequency Ashtech Z-XII and ZXtreme receivers with Geodetic III antennas. The installation of six sites was completed on June 16, 2003. Most stations remained operational for 2.5 yr. In order to capture the early post-seismic signal on a tight schedule, antennas were first installed on tripods. In July 2003, monumentation was upgraded to concrete pillars. Several days of simultaneous observations were acquired at both monuments in order to establish an accurate tie at each site.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Post-seismic displacements measured at continuous GPS station in the epicentral area of the 2003 May 21, Boumerdes earthquake, shown with respect to site EMAP. Ellipses are 95 per cent confidence. Top panel: Cumulative displacements for the May 2003 to October 2005 period. Numbers by the arrow tails are the cumulative vertical displacement in millimetres. Circles show aftershocks from 2003 May 25–30 (Bounif *et al.* 2004) colour-coded as a function of depth. Focal mechanism and centroid location are from Delouis *et al.* (2004). Bottom panel: Displacements for four 7-month long successive time periods.</figcaption></figure>

We processed phase and pseudo-range GPS data using the GAMIT software (version 10.2; King & Bock 2005), solving for station coordinates, satellite state vectors, 7 tropospheric delay parameters per site and day, and phase ambiguities using double-differenced GPS phase measurements. We used the final orbits from the International GNSS Service (IGS), earth orientation parameters from the International Earth Rotation Service (IERS), and applied azimuth and elevation dependent antenna phase centre models, following the tables recommended by the IGS. We included 10 global IGS stations in Europe and Africa to serve as ties with the terrestrial reference frame (ITRF2000; Altamimi *et al.* 2002). The least squares adjustment vector and its corresponding variance-covariance matrix for station positions and orbital elements estimated for each independent daily solution were then combined with global SINEX (Solution Independent Exchange format) files from the IGS daily processing routinely done at the Scripps Orbital and Permanent Array Center (http://sopac.ucsd.edu). The reference frame is implemented using this unconstrained combined solution by minimizing the position and velocity deviations of 41 IGS core stations with respect to the ITRF2000 while estimating an orientation, translation and scale transformation. In this process, height coordinates were downweighted using a variance scaling factor of 10 compared to the horizontal components.

We converted the resulting position time-series from ITRF2000 to a Nubia-fixed frame using the Nubia-ITRF2000 angular rotation from Calais *et al.* (2003). We find that EMAP, contrary to all other sites, shows a linear displacement as a function of time, with no indication of time-dependent transient. Velocity at EMAP for the 2 yr measurement time span is $1.5 \pm 4 \text{ mm yr}^{-1}$, with a relatively large uncertainty due to the short observation period. Assuming that the entire Nubia-Eurasia plate motion at the longitude of Boumerdes (5.5 mm yr$^{-1}$) accumulates as elastic strain on the offshore fault system dipping south at 42° and locked to a depth of 20 km results in $2 \text{ mm yr}^{-1}$ of expected velocity at EMAP with respect to Nubia. In other words, there is no evidence of accelerated displacement at site EMAP after the Boumerdes earthquake compared to expected secular rates. We therefore, mapped horizontal displacements with respect to local site EMAP (displacements shown on Fig. 2,
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Position time-series for sites CDJP, DLYP, SFNP, THNP and BJMP. Dots show daily positions corrected for an annual and a biannual periodic term and for instrument offsets. Red curves show best logarithmic fit.</figcaption></figure>

time-series shown on Fig. 3), which has the advantage of significantly increasing the signal-to-noise ratio of the position time-series over relatively short baselines compared to absolute position time-series expressed in a global reference frame.

In a second step, we cleaned the position time-series by estimating and removing annual and semi-annual terms, likely to result from unmodelled non-tectonic effects such as hydrological or atmospheric loading, and instrument offsets. To do so, we modelled site positions as the sum of (1) a linear term representing secular elastic strain accumulation, (2) a logarithmic term representing

post-seismic deformation (Langbein et al. 2006), (3) an annual and semi-annual periodic term representing seasonal effects not modelled in the GPS data analysis and (4) DC offsets due to equipment changes at the site. The model equation is:

$$y = at + b + \sum_{i=1}^{n} c_i H_i(t) + d \sin(2\pi t) + e \cos(2\pi t) + f \sin(4\pi t) + g \cos(4\pi t) + h \ln(1 + t/\tau), \quad (1)$$

where a, b, c, d, e, f, g and h are estimated by inverting the site position data y using a singular value decomposition scheme.
<table><caption>**Table 1.** Cumulative post-seismic displacements in millimetres (east, north, up) and associated 1-σ error.</caption><thead><tr><th>Site</th><th>Lon.</th><th>Lat.</th><th>East</th><th>North</th><th>Up</th><th>σ<sub>E</sub></th><th>σ<sub>N</sub></th><th>σ<sub>U</sub></th></tr></thead><tbody><tr><td>EMAP</td><td>3.253</td><td>36.810</td><td>0.0</td><td>0.0</td><td>0.0</td><td>1.0</td><td>1.0</td><td>4.0</td></tr><tr><td>BJMP</td><td>3.723</td><td>36.733</td><td>-12.7</td><td>25.6</td><td>21.4</td><td>0.9</td><td>0.9</td><td>4.3</td></tr><tr><td>CDJP</td><td>3.720</td><td>36.878</td><td>-19.5</td><td>36.7</td><td>26.1</td><td>0.4</td><td>0.4</td><td>2.0</td></tr><tr><td>DLYP</td><td>3.892</td><td>36.916</td><td>-10.0</td><td>31.0</td><td>-1.7</td><td>0.4</td><td>0.4</td><td>1.9</td></tr><tr><td>SFNP</td><td>3.556</td><td>36.786</td><td>-20.2</td><td>28.4</td><td>-27.0</td><td>0.4</td><td>0.4</td><td>1.8</td></tr><tr><td>THNP</td><td>3.560</td><td>36.724</td><td>0.5</td><td>27.4</td><td>-12.7</td><td>0.5</td><td>0.5</td><td>2.3</td></tr></tbody></table>

Parameter $τ$ is found by iteratively minimizing the model misfit over a range of values. $H_i(t)$ is a binary operator equal to 0 or 1 before or after offset $i$, respectively. The average amplitude of the annual signal of 3.4 mm on the vertical and 1.5 mm on horizontal components. The amplitude of the biannual signal is about half of the annual one. The logarithmic time constant $τ$ is on the order of 0.3 yr (~110 d for the horizontal components. It is smaller (~0.1 yr or 36 d), but also less well determined, for the noisier vertical component. For BJMP (north component) and THNP (east and up components), the scatter in position time-series together with the relatively short time span they cover does not allow us to reliably estimate a logarithmic decay term.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Post-seismic slip distribution estimated from cumulative displacements at GPS sites. Black arrows = observations; red arrows = model predictions. Top four panels: Slip distribution for four consecutive time periods of 7 months each. Bottom: Cumulative slip distribution for the entire observation time period (2.5 yr). Numbers by the site names show the vertical displacements in millimetres: model/observed.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Resolution test of the slip distribution for deep versus shallow slip. Top panel: input slip and synthetic GPS displacements used in the inversion. Bottom panel: estimated slip and predicted GPS displacements.</figcaption></figure>

Once estimated, we remove the DC offsets and annual and semi-annual terms, keeping the secular and logarithmic terms as the constraint set for our post-seismic study. One-sigma uncertainties on cumulative displacements are taken as the rms of the scatter of the position data about the model presented above. The resulting post-seismic horizontal displacements (Fig. 2, Table 1) show up to 40 mm of horizontal displacement in 2.5 yr in the same direction as coseismic displacements predicted using Delouis et al.'s (2004) rupture model. Horizontal displacements at the coastal sites are larger during the first 7 months following the event, after which they gradually slow down. Horizontal displacements at the inland sites appear more stable with time. Numerous data gaps and position scatter at site BJMP after mid-2004 lead to larger uncertainties on the estimated displacements after that time. The position time-series (Fig. 3) show that post-seismic deformation continues at a significant rate at all the sites after the GPS stations were removed in October 2005.

## 4 AFTERSLIP INVERSION

We first model post-seismic deformation as the result of continued slip on the coseismic rupture plane. We use the same fault geometry as in Delouis et al. (2004), with a N70°E strike and 40° dip. Model

fault strike is consistent with the azimuth of active reverse faults mapped offshore (Déverchère et al. 2005) and dip with the depth distribution of aftershocks (Bounif et al. 2004). We extend the fault from the surface (offshore) to a depth of 20 km, about 10 km below the area where coseismic slip has been reported, in order to allow for downdip afterslip. The modelled fault captures the entire area affected by aftershocks (Bounif et al. 2004).

We invert the GPS displacements (Table 1) for slip of the fault described above. We divide the fault into 20 (along-strike) by 8 (downdip) rectangular patches of dimension 2.8 × 3.7 km each and compute Green's functions G relating slip s on each patch to the 3-D displacement u at GPS sites assuming an elastic half-space with a Poisson ratio of 0.25 (Okada 1985). We solve for dip-slip motion only, which results in 160 estimated parameters for 36 observations. We regularize the inversion by applying smoothing via a finite difference approximation of the Laplacian operator (∇²) and an associated weighting factor κ. We apply positivity constraints in order to avoid implausible and overly rough slip distributions. The inversion minimizes the weighted residual sum of squares (WRSS) and the roughness of the slip model using the functional:

$$ \|W(G\bar{s} - \bar{u})\|^2 + \kappa^2 \|\nabla^2\bar{s}\|^2, \qquad (2) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Resolution test of the slip distribution for along-strike resolution at shallow depth. Top panel: input slip and synthetic GPS displacements used in the inversion. Bottom panel: estimated slip and predicted GPS displacements.</figcaption></figure>

where *W* is the weight matrix associated with the GPS data. The value of *κ* that optimizes smoothing versus data fitting is found iteratively using a cross-validation technique (Matthews & Segall 1993).

The inversion of the GPS data results in a reduced $\chi^2$ of 1.8, with a good fit to the horizontal displacements (weighted rms = 2.4 mm). The fit to the vertical displacement is not as good, with a weighted rms of 12.2 mm, due for a large part to site BJMP which also has the largest uncertainty (See Fig. 3 and Table 1) and the poorest time-series, with a large scatter and significant data gaps. The modelled slip distribution (Fig. 4) for the entire 2.5 yr period shows that slip is concentrated at shallow depths, with maximum amplitude (up to 25 cm) in the upper 5 km. Little to no slip is estimated below 10 km. The equivalent moment release estimated for that time period is $3.23 \times 10^{18}$ N m ($M_w$ = 6.3). The slip distribution for four successive 7-month periods (Fig. 4) shows a similar pattern and indicates

that most of the slip occurs during the first 7 months following the earthquake, with an equivalent moment release of $2.4 \times 10^{18}$ N m ($M_w$ = 6.2). Subsequent time periods show a significant decrease in slip and an associated progressive decrease of the equivalent moment release from $8.3 \times 10^{17}$ to $5.3 \times 10^{17}$ N m.

We performed resolution tests to determine how well the spatial distribution of slip is estimated from our data set. In particular, we want to assess whether the inversion could have missed slip downdip of the coseismic rupture and how well shallow slip is resolved. In a first test, we use an input slip model consisting of two fault-parallel stripes with uniform slip (0.1 m of pure dip-slip per patch). The shallow one extends between 0 and 5 km, the deeper one between 12.5 and 17.5 km (Fig. 5, top). This slip model was used in a forward run to produce synthetic surface displacements at the six GPS sites, to which we added random noise with a ±5 mm standard deviation. The inversion of this synthetic data set (Fig. 5, bottom)
shows that the two areas of slip are recovered, with maxima spatially consistent with the input model. The solution however results in a significant smearing of the input model, in particular at depth. Input slip amplitudes are well-recovered overall for the shallow area, but are underestimated by about 50 per cent at depth. In addition, estimated slip on the shallow patch shows an asymmetric pattern with larger amplitudes in its western half, which does not match the input model. In a second test, the input slip model consists of seven shallow slip patches with 0.1 m of pure dip-slip per patch (Fig. 6). The inversion recovers the overall amount of slip (input magnitude is 6.1, recovered magnitude is 6.1) and the depth range, but individual patches are poorly recovered. These tests highlight the poor along-strike resolution of the inversion due to the very sparse GPS data set.

In a third test, we solved for afterslip on a fault starting at 7.5 km depth (i.e. just below the area of maximum slip found above) and extending to 20 km. The fit to the data is significantly poorer than when allowing for shallow slip, and the inversion still assigns larger slip to the shallower parts of the fault. With the resolution limitations listed above in mind, our data set is, therefore, consistent

with shallow slip (with a poorly resolved along-strike variability) and little to no deep slip over the time period studied here.

## 5 VISCOELASTIC MODELLING

We next test whether viscoelastic relaxation in the lower crust and/or upper mantle could explain the observed post-seismic deformation. To do so, we developed a 3-D viscoelastic finite element model (FEM) of the rupture surface and surrounding region that simulated the inferred coseismic slip distribution (Fig. 7a), and allowed a viscoelastic lower crust and/or upper mantle to relax for the duration of GPS post-seismic observations. We utilized the finite element software Ideas (www.eds.com) that we have used in several previous post-seismic studies (e.g. Freed & Bürgmann 2004; Freed et al. 2006). Coseismic slip is simulated by enforcing the slip distribution of Delouis et al. (2004) on an explicit fault in the mesh. Fig. 7(b) shows a reasonable fit between FEM predicted and GPS observed coseismic surface displacements from Yelles et al. (2004).

Calculated maximum shear stress changes induced by coseismic slip (Fig. 7a), available to drive relaxation of a viscoelastic

<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Cutaway view of finite element mesh showing a portion of the footwall plate and calculated maximum shear stress changes induced by the Boumerdes quake. The rupture surface is indicated with the black/white dashed boundary. The full model is 500 km along strike by 300 km normal to strike, with a depth of 100 km. All side and bottom boundaries are fixed, which does not influence the stress state generated by coseismic slip (i.e. does not influence post-seismic results). The full model contains 30 000 elements. Coseismic slip is generated by enforcing displacements between coincident nodes along the rupture plane in accordance with the inferred slip distribution (Delouis et al. 2004). Enforced slip generates coseismic stress changes and deformation of the surface. (b) Comparison of observed GPS versus calculated coseismic surface displacements from the finite element model show a reasonable match.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of GPS observed and viscoelastic model predicted (a) horizontal surface displacements for lower-crustal and upper-mantle flow models and (b) vertical surface displacements for lower crustal flow model. All displacements are relative to station EMAP.</figcaption></figure>

lower crust and upper mantle, are significant in the upper crust (~1 MPa). They are much smaller at lower-crustal levels, with ~0.3 MPa at 22 km depth and ~0.1 MPa at the 32-km-deep base of the crust. These stress magnitudes are relatively small and, compared with other well-studied events (e.g. Freed & Bürgmann 2004; Freed et al. 2006), may not drive detectable viscoelastic flow in the lower crust and upper mantle in the time frame of the GPS observations.

To explore possible contributions of viscoelastic relaxation, we considered models of only lower crustal flow (22–32 km depth), only upper-mantle flow (below 32 km depth), and models of combined lower-crust and upper-mantle flow. Each calculation was conducted by applying the coseismic slip distribution, then allowing the viscoelastic region to relax for 2.5 yr. For each configuration we assumed a range of viscosities ($10^{17}$–$10^{19}$ Pa s) and searched for the viscosity that leads to the minimum misfit with respect to GPS observed horizontal surface displacements without overpredicting any

of the displacements. This last criteria is critical, as overpredicted horizontal displacements are unlikely to be countered by afterslip unless it is in the direction opposite to that of the coseismic slip, which is highly unlikely.

We find that models matching horizontal displacements at coastal GPS stations (CDJP, DLYP, SFNP) greatly overpredict displacements at the inland stations (BJMP, THNP). The best-fitting models are those that match horizontal displacements at the inland stations, although they underpredict displacements at the coastal stations (Fig. 8). To best match inland stations with a model of lower crustal flow only, we find a best-fitting viscosity of $6 \times 10^{17}$ Pa s. For a model of upper-mantle flow only, the best-fitting viscosity is $3 \times 10^{17}$ Pa s. The magnitude of displacements at inland stations can also be matched by considering relaxation in both the lower crust and upper mantle with a uniform viscosity of $1.4 \times 10^{18}$ Pa s. We also considered a model in which the lower crust had a viscosity five times that of the upper mantle (consistent with temperature
<figure><img src="image_10.png" /><figcaption>**Figure 9.** North–south cross-section of the Algerian margin offshore Boumerdes illustrating the spatial relationship between the coseismic rupture (western half of the rupture plane), the area of afterslip estimated here, and mapped faults along the margin. No vertical exaggeration. Modified from Dèverchère *et al.* (2004).</figcaption></figure>

increase with depth), and found that inland GPS displacements could be matched with a $3.5 \times 10^{18}$ Pa s lower crust and a $7.0 \times 10^{17}$ Pa s upper mantle.

While any of the viscoelastic structures considered can be used to match the amplitude of displacements at inland sites, the lower crustal flow model provides the closest match to their azimuths while the upper-mantle flow model leads to the largest errors (Fig. 8a). In addition, the upper-mantle flow model shows significant displacements at coastal stations CDJP and DLYP with large azimuth differences to those observed (though the azimuth match at station SFNP is very good). While this difference could potentially be made up by afterslip, it would require much more shear (strike-slip motion) than was observed in the coseismic slip distribution, which is counter-intuitive. Vertical displacements from all viscoelastic relaxation models predict subsidence over an area encompassing the rupture zone (Fig. 8b), becoming modestly broader with the deeper flow models. None of these models explain the observed uplift at sites CDJP and BJMP, and all greatly overpredict observed subsidence at THNP, DLYP and SFNP relative to EMAP. Overall, post-seismic displacements are better fit by flow in the lower crust than in the upper mantle, although none of the models tested here match the entire data set.

## 6 DISCUSSION

As shown above, none of the viscoelastic flow models tested provide a satisfactory match to the observed horizontal or vertical GPS displacements. In addition, viscosities in the best-fitting flow models are an order of magnitude lower than found in other studies where post-seismic crustal or upper-mantle flow is well documented (e.g. Pollitz *et al.* 2001; Freed & Bürgmann 2004; Freed *et al.* 2006).

We also tested whether coseismic pressure changes in the shallow crust (i.e. poroelastic rebound) could explain the GPS observations. We modelled this process using a variation in Poisson's ratio from an undrained condition to a drained condition in which fluid pressure

equilibrium is re-established, while using the same shear modulus (Roeloffs 1996; Jónsson *et al.* 2003). Using a similar reduction of Poisson's ratio as proposed for the Landers event (from 0.25 to 0.21; Peltzer *et al.* 1998), we find post-seismic surface displacements up to a maximum of 5 mm at any of the stations, about an order of magnitude smaller than the observations. In addition, the azimuth of calculated poroelastic rebound displacements are orthogonal to the direction of observed post-seismic displacements. Finally, poroelastic rebound produces primarily uplift, with up to 10 mm at station THNP, where subsidence was observed. Thus, although poroelastic rebound may be occurring since fluids in the upper crust are ubiquitous, it provides at most a small contribution to observed post-seismic surface deformation.

Afterslip, therefore, appears to be the primary post-seismic process at work in the 2.5 yr following the Boumerdes event. It concentrates updip of the coseismic slip patch reported by Delouis *et al.* (2004) and Semanne *et al.* (2005) in the western half of the rupture. In the eastern half of the rupture, afterslip also concentrates updip of the coseismic slip patch of Semanne *et al.* (2005) but coincides with the coseismic patch found by Delouis *et al.* (2004). The lack of afterslip downdip of the coseismic rupture is an interesting feature of the Boumerdes event. For instance, post-seismic deformation following the 1999, $M_w = 7.3$ Chi-Chi earthquake in Taiwan, another example of thrust event in a convergent setting, shows that afterslip is concentrated south and downdip of the areas of largest coseismic slip (Hsu *et al.* 2002). Hutton *et al.* (2001) find a similar result for the 1995 $M_w = 8.0$ Jalisco earthquake, with afterslip concentrated below 20 km and migrating downward with time along the plate interface. Afterslip following the 2003, $M_w = 8.0$, Tokachi-Oki is also localized outside of the area of large coseismic slip, although mostly along-track of it (Miyazaki *et al.* 2004).

In the framework of rate-and-state friction theory (Dieterich 1979; Scholz 1998), the absence of afterslip downdip of the coseismic rupture may indicate a velocity-weakening behaviour in the lower crust, which responds to coseismic slip by a sudden decrease in
frictional stress but remains strongly coupled between earthquakes. Reciprocally, the concentration of afterslip outside of the coseismic rupture area may be indicative of a velocity-strengthening behaviour, where the area updip of the rupture responds to coseismic slip by short-term increase in frictional stress that, in turn, induces afterslip. Since this behaviour prevents accelerated slip, stresses in velocity-strengthening areas have to be released by processes other than earthquake rupture, such as afterslip or interseismic creep (e.g. Heki et al. 1997).

Although there is no evidence for shallow creep along the Boumerdes fault, there is ample data showing soft deformation through folding in the sediments at the toe of the Algerian margin slope. Seismic reflection surveys coupled with detailed bathymetry have shown that most of the Algerian margin, including offshore the Boumerdes area, is characterized by a series of active thrust faults and fault-propagation folds (Déverchère et al. 2005). The depth range at which afterslip concentrates (0–5 km) actually corresponds to the estimated depth of margin sediments (Fig. 9; Domzig et al. 2006). Post-seismic deformation (and potentially interseismic creep) may, therefore, involve folding of sedimentary layers and interbedding slip, a hypothesis previously mentioned by Donnellan & Lyzenga (1998) to explain shallow afterslip following the $M_w = 6.7$, 1994, Northridge earthquake in California.

## 7 CONCLUSIONS

Continuous GPS data collected at six sites in the Boumerdes area for the 2.5 yr following a $M_w = 6.9$ thrust event show clear evidence for post-seismic deformation with up to 4 cm of cumulative horizontal displacement and a time-dependence well fit by a logarithmic decay. We find that the data are consistent with shallow afterslip (0–5 km) but show no evidence for afterslip downdip of the coseismic rupture or for a significant contribution of viscoelastic relaxation in the lower crust or upper mantle, or poroelastic rebound. Afterslip is rapid during the first seven months or so, then decays significantly. Post-seismic deformation appears to continue at all sites after the 2.5-yr observation period, with rates less than $1 \text{ cm yr}^{-1}$.

The fact that afterslip concentrates adjacent to and updip of the coseismic rupture in the western half of the fault suggests that afterslip is driven by coseismic stresses, a mechanism proposed for several other thrust events (e.g. Hutton et al. 2001; Hsu et al. 2002; Miyazaki et al. 2004). We find that the area of afterslip corresponds to the depth range of the folded sediments of the margin while coseismic slip occurred below that depth and affected basement rocks. This suggests that spatial variations in frictional properties along the fault correlate with the type of rocks involved.

## ACKNOWLEDGMENTS

This work benefited from technical and field support from W. Bacha, A. Bellik (CRAG) and O. Charade (IPGP). We thank the Institut National des Sciences de l'Univers (INSU) for a long-term loan of GPS instruments for the experiment and the Algerian institutions and individuals that agreed to host the semi-continuous GPS stations. We are grateful to J. Déverchère for his encouragements and for discussions on the seismotectonics of the Algerian margin. We thank F. Cotton and an anonymous reviewer for their helpful comments on the first version of this paper. This work was supported by CMEP (Comité Mixte d'Evaluation et de Prospective de coopération interuniversitaire franco-algérienne) - Programme TASSILI - N. 041 MDU 619 and the Algerian PNR III Programme.

## REFERENCES

Altamimi, Z., Sillard, P. & Boucher, C., 2002. ITRF2000: a new release of the International terrestrial reference frame for earth science applications, *J. geophys. Res.*, **107**, 2214, doi:10.1029/2001JB000561.

Ayadi, A. et al., 2003. Strong Algerian earthquake strikes near capital city, *EOS, Trans. Am. geophys. Un.*, **84**, 50, 561–568.

Bounif, A. et al., 2004. The 21 May 2003 Zemmouri (Algeria) earthquake Mw 6.8: relocation and aftershock sequence analysis, *Geophys. Res. Lett.*, **31**, L19606, doi:10.1029/2004GL020586.

Calais, E., DeMets, C. & Nocquet, J.M., 2003. Evidence for a post-3.16 Ma change in Nubia-Eurasia plate motion, *Earth planet. Sci. Lett.*, **216**(81-92), doi:10.1016/S0012-821X(03)00482-5.

Delouis, B. et al., 2004. Slip distribution of the 2003 Boumerdes-Zemmouri earthquake, Algeria, from teleseismic, GPS, and coastal uplift data, *Geo-phys. Res. Lett.*, **31**, L18607, doi:10.1029/2004GL020687.

Déverchère, J. et al., 2005. Active thrust faulting offshore Boumerdes, Algeria, and its relations to the 2003 Mw 6.9 earthquake, *Geophys. Res. Lett.*, **32**, L04311, doi:10.1029/2004GL021646.

Dieterich, J.H., 1979. Modelling of rock friction, 1: experimental results and constitutive equations, *J. geophys. Res.*, **84**, 2161–2168.

Domzig, A., Yelles, K., Le Roy, E., Déverchère, J., Bouillin, P., Bracène, R. et al., 2006. Searching for the Africa-Eurasia Miocene boundary offshore western Algeria (MARADJA'03 cruise), *C. R. Géoscience*, **338**(1-2), 80-95.

Donnellan, A. & Lyzenga, G., 1998. GPS observations of fault afterslip and upper crustal deformation following the Northridge earthquake, *J. geophys. Res.*, **103**(9), 21 285–21 297.

Freed, A.M. & Bürgmann, R., 2004. Evidence of power-law flow in the Mojave Desert mantle, *Nature*, **430**, 548–551.

Freed, A.M., Bürgmann, R., Calais, E., Freymueller, J. & Hreinsdóttir, S., 2006. Implications of deformation following the 2002 Denali, Alaska earthquake for postseismic relaxation processes and lithospheric rheology, *J. geophys. Res.*, **111**, doi:10.1029/2005JB003894.

Hauksson, E., Jones, L.M. & Hutton, K., 1995. The 1994 Northridge earthquake sequence in California: seismological and tectonic aspects, *J. geophys. Res.*, **100**, 12 335–12 356.

Heki, K., Myazaki, S. & Tsujii, H., 1997. Silent fault slip following an interplate thrust earthquake at the Japan trench, *Nature*, **386**, 595–598.

Hsu, Y.J., Bechor, N., Segall, P., Yu, S.B., Kuo, L.C. & Ma, K.F., 2002. Rapid afterslip following the 1999 Chi-Chi Taiwan Eqrthquake, *J. geophys. Lett.*, **29**, 16, doi:10.1029/2002GL01467.

Hutton, W., DeMets, C., Sànchez, O., Suárez, G. & Stock, J., 2001. Slip dynamics during and after the 9 October 1995 Mw = 8.0 Colima-Jalisco, Mexico, *Geophys. J. Int.*, **146**, 637–658.

Jónsson, S., Segall, P., Pedersen, R. & Björnsson, G., 2003. Post-earthquake ground movements correlated to pore-pressure transients, *Nature*, **424**, 179–183.

King, R.W. & Bock, Y., 2005. Documentation for GAMIT GPS processing software Release 10.2, Mass. Inst. of Technol., Cambridge, MA.

Langbein, J., Murray, J.R. & Snyder, H.A., 2006. Coseismic and initial post-seismic deformation from the 2004 Parkfield, California, earthquake, observed by global positioning system, electronic distance meter, creepmeters, and borehole strainmeters, *Bull. seism. Soc. Am.*, **96**(4B), S304-S320.

Lin, J. & Stein, R.S., 1989. Coseismic folding, earthquake recurrence, and the 1987 source mechanism at Whittier Narrows, Los Angeles basin, California, *J. geophys. Res.*, **94**, 9614–9632.

Matthews, M.V. & Segall, P., 1993. Statistical inversion of crustal deformation data and estimation of the depth distribution of slip in the 1906 earthquake, *J. geophys. Res.*, **98**, 12 153–12 163.

Meghraoui, M., 1991. Blind reverse faulting system associated with the Mont Chenoua-Tipaza earthquake of 29 October 1989 (north-central Algeria), *Terra Nova*, **3**, 84–93.

Meghraoui, M., Morel, J.L., Andrieu, J. & Dahmani, M., 1996. Tectonique plio-quaternaire de la charne tello-rifaine et de la mer d'Alboran. Une
zone complexe de convergence continent-continent, *Bull. Soc. Géol. Fr.*, **167**, 141–157.

Miyazaki, S. et al., 2004. Modeling the rupture process of the 2003 September 25 Tokachi-Oki (Hokkaido) earthquake using 1-Hz GPS data, *Geophys. Res. Lett.*, **31**, doi:10.1029/2004GL021457.

Nocquet, J.M. & Calais, E., 2004. Geodetic measurements of crustal deformation in the western Mediterranean and Europe, *Pure appl. Geophys.*, **161**, 661–681, doi:10.1007/s00024-003-2468-z.

Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space, *Bull. seism. Soc. Am.*, **75**, 1135–1154.

Peltzer, G., Rosen, P., Rogez, F. & Hudnut, K., 1998. Poroelastic rebound along the Landers 1992 earthquake surface rupture, *J. geophys. Res.*, **103**, 30131–30145.

Peltzer, G., Rosen, P., Rogez, F. & Hudnut, K., 1996. Postseismic rebound in fault step-overs caused by pore fluid flow, *Science*, **273**, 1202–1204.

Pollitz, F.F., Wicks, C. & Thatcher, W., 2001. Mantle Flow Beneath a Continental Strike-Slip Fault: postseismic Deformation After the 1999 Hector Mine Earthquake, *Science*, **293**, 1814–1818.

Roeloffs, E., 1996. Poroelastic techniques in the study of earthquake-related hydrological phenomena, *Adv. Geophys.*, **37**, 135–195.

Savage, J.C. & Svarc, J.L., 1997. Postseismic deformation associated with the 1992 Mw = 7.3 Landers earthquake, southern California, *J. geophys. Res.*, **102**, 7565–7577.

Semanne, F., Campillo, M. & Cotton, F., 2005. Fault location and source process of the Boumerdes, Algeria, earthquake inferred from geodetic and strong motion data, *Geophys. Res. Lett.*, **32**, L01305, doi:10.1029/2004GRL021268.

Shen, Z.K., Jackson, D., Feng, Y., Cline, M., Kim, M., Fang, P. & Bock, Y., 1994. Postseismic deformation following the Landers earthquake, California, 28 June 1992, *Bull. seism. Soc. Am.*, **84**, 780–791.

Scholz, C.H., 1998. Earthquakes and friction laws, *Nature*, **39**, 37–42.

Thatcher, W., 1974. Strain release mechanism of the 1906 San Francisco earthquake, *Science*, **184**, 1283–1285.

Thatcher, W., 1983. Non-linear strain buildup and the earthquake cycle on the San Andreas fault, *J. geophys. Res.*, **88**, 5893–5902.

Yelles, K., Lammali, K., Mahsas, A., Calais, E. & Briole, P., 2004. Coseismic deformation of the May 21st, 2003 earthquake, Algeria, from GPS measurements, *Geophys. Res. Lett.*, **31**, L13610, doi:10.1029/2004GL019884.

Yielding, G., Ouyed, M., King, G.C.P. & Hatzfeld, D., 1989. Active tectonics of Algerian Atlas Mountains—evidence from aftershocks of the 1980 El-Asnam earthquake, *Geophys. J. Int.*, **99**(3), 761–778.