The 2013 September–October seismic sequence offshore Spain: a case of seismicity triggered by gas injection?

Abstract

A spatially localized seismic sequence originated few tens of kilometres offshore the Mediterranean coast of Spain, close to the Ebro river delta, starting on 2013 September 5, and lasting at least until 2013 October. The sequence culminated in a maximal moment magnitude M_w 4.3 earthquake, on 2013 October 1. The most relevant seismogenic feature in the area is the Fosa de Amposta fault system, which includes different strands mapped at different distances to the coast, with a general NE–SW orientation, roughly parallel to the coastline. However, no significant known historical seismicity has involved this fault system in the past. The epicentral region is also located near the offshore platform of the Castor project, where gas is conducted through a pipeline from mainland and where it was recently injected in a depleted oil reservoir, at about 2 km depth. We analyse the temporal evolution of the seismic sequence and use full waveform techniques to derive absolute and relative locations, estimate depths and focal mechanisms for the largest events in the sequence (with magnitude mbLg larger than 3), and compare them to a previous event (2012 April 8, mbLg 3.3) taking place in the same region prior to the gas injection. Moment tensor inversion results show that the overall seismicity in this sequence is characterized by oblique mechanisms with a normal fault component, with a 30° low-dip angle plane oriented NNE–SSW and a subvertical plane oriented NW–SE. The combined analysis of hypocentral location and focal mechanisms could indicate that the seismic sequence corresponds to rupture processes along shallow low-dip surfaces, which could have been triggered by the gas injection in the reservoir, and excludes the activation of the Amposta fault, as its known orientation is inconsistent with focal mechanism results. An alternative scenario includes the iterated triggering of a system of steep faults oriented NW–SE, which were identified by prior marine seismics investigations.

Earthquake dynamics, Earthquake source observations

1 INTRODUCTION

On 2013 September 5, the seismicity rate suddenly increased at a localized region located offshore the Mediterranean coast of Spain, at about Lat 40.3°N, Lon 0.7°E. The seismic sequence remained active for more than one month and was still ongoing in 2013 October. Several institutes run seismic stations and monitoring systems in the region, including the Spanish National Geographic Institute (IGN), the Catalan Geological Institute (IGC) and the Ebro Observatory (EBR). By 2013 October 15, 511 events had been detected and located by IGN, with magnitudes between 0.7 and 4.2 (mbLg), and 1002 by EBR, with local magnitudes between 0.0 and 4.0. The largest event (magnitudes 4.0–4.2) struck on 2013 October 1, at 03:32:45 UT, with preliminary locations at Lat 40.32°N, Lon 0.79°E according to IGN (Lat 40.37°N, Lon 0.66°E in the EBR catalogue). Source depths are typically shallow, with most hypocentres located at less than 10 km depth. According to these locations, earthquake foci are spatially distributed along an elongated region (Fig. 1), roughly extending along a NW–SE direction, almost perpendicular to the coastline. Such spatial distribution could suggest the geometry of the involved rupture region, but could also reflect the spatial extension of location uncertainties, given that most of the closest stations are located inland, towards NW. The location quality is strongly limited by the absence of close stations and the strong network asymmetry, with the average major semi-axis of the location error ellipse of about 7 km (EBR catalogue). The azimuthal gap is partially closed towards NE, thanks to the ocean bottom station COBS (IGC), and towards SW, with stations EIBI and ETOS (IGN) located on Ibiza and Mallorca islands. These stations are, however, at larger distance, at more than 150 km from the epicentral region. The latter (here limited within Lon 40.00–40.80°N, and Lat 0.25–1.25°E) has been affected in the past by sporadic, low-magnitude earthquakes. The largest one in the IGN catalogue has a magnitude mbLg 3.3 and was recorded on 2012 April 8. Considering that the magnitude of completeness (M_c) for the epicentral region has been estimated in M_c ∼ 2 since 2002 March (González, in preparation), the IGN catalogue includes only seven events in the 12-year period 2002 March–2013 February (less than 1 events yr^–1) with magnitude above completeness. On the contrary 58 events have been included in the same catalogue in the first 32 d starting on 2013 September 5 (>660 events yr^–1). The seismic activity of 2013 September and October was unusual in terms of maximal magnitude and seismicity rate, if compared to instrumentally recorded seismicity occurring in the last two decades. According to the review by Martínez Solares (2003), no relevant historical seismicity is known in the area. However, a comparison with longer-term past seismicity is not possible as instrumental and historical catalogues are incomplete and have poor locations and magnitudes estimations for moderate earthquakes with offshore epicentres.

Figure 1.

Spatial distribution of epicentres and operating stations in the area. Circles: epicentres of the studied seismic sequence offshore Spain (according to the catalogue of the Ebro Observatory). Green stars: epicentres of historical seismicity (with maximum intensity I > V in the IGN catalogue). Red stars: epicentres of previous instrumental seismicity (with M ≥ 4 in the IGN catalogue). White square: location of the Castor platform of gas injection. Triangles: seismic stations (either with broad-band or short period sensors, as well as the COBS ocean bottom seismometer) operated by IGN, IGC and EBR in the region.

The epicentral area is located in the Gulf of Valencia, within the Catalan-Valencian extensional domain. The extensional tectonics is accommodated through complex normal fault systems (Perea et al.2012). The main faults are oriented ENE–WSW in the north and NNE–SSW to N–S in the south (Roca & Guimerà 1992), with most faults oriented parallel to the coast. Three strands (Western, Central and Eastern strands) of the Fosa de Amposta fault system are located in proximity of the epicentral region (Roca 1992, 1996; Roca & Guimerà 1992; Perea 2006; Perea et al.2012). The three strands are building up a system of normal faulting of different dip orientation. The strike of these geological features would support the interpretation that the NW–SE distribution of seismicity is rather due to locations uncertainties. The offshore Eastern strand, which is the closest to the preliminary located epicentres of the largest events in the sequence has a variable NE–SW strike (172–200°) and dip towards the coastline (dip ca. 60° towards NW). The central and Western strands are also oriented parallel to the coast (strike 25° and 27°, respectively) but dip towards SE (dip ca. 60° for both strands).

The seismic sequence raised great interest among the scientific community and civil society, given its temporal coincidence with a test for cushion gas injection in the proximity of the epicentral region. This test was performed by the Castor project, which aims to store natural gas (Batchelor et al.2007) in the depleted Amposta oil reservoir, exploited in the period 1973–1989 (Seeman et al.1990; Batchelor et al.2007). The injection phase took place from 2 to 16 September at 1750 m under the sea level. Prior test injections of cushion gas were performed since 2013 June, not accompanied by seismicity. The total injected volume of gas was ≈1.02 × 10⁸ m³, measured at the standard conditions of 25 °C and 1 bar according to the official protocols for the Spanish underground gas repositories (Ministerio de Industria, Turismo y Comercio 2007). The seismicity rate increased first on September 5, consistently with the operation of the last injection phase. However, the largest events in the seismic sequence took place from September 24 on, after the injection was stopped. The spatial vicinity among the injection well and the epicentres, the shallow seismicity depth, and the beginning of the activity during the injection phase raise the question of whether the gas storage operations could have triggered or induced the increased seismic activity.

Seismicity can be triggered (meaning that its occurrence is anticipated because of a stress perturbation) or induced (in this case the seismicity is controlled in nucleation and size by the perturbing stress, and would not take place without it) by a variety of human related operations (Cesca et al.2013a; Dahm et al.2013). It has been shown that seismicity can be triggered or induced by mining (e.g. Feignier & Young 1992; McGarr 1992; Trifu et al.2000; Fletcher & McGarr 2005; Bischoff et al.2010; Sen et al.2013), water reservoirs impoundment (Gough & Gough 1970; Gupta & Rastogi 1976; Guha 2000; Do Nascimento et al.2005), wastewater or fluid injection also including brine from hydraulic fracturing of shale gas (Healy et al.1968; Cornet & Jianmin 1995; Sasaki 1998; Tadokoro et al.2000; Ake et al.2005; Sileny et al.2009; Ellsworth 2013; Kim 2013), geothermal fields operation (e.g. Ross et al.1999; Deichmann & Giardini 2009; Brodsky & Lajoie 2013) and oil and gas field exploitation (e.g. Grasso & Wittlinger 1990; Rutledge et al.1998; Dahm et al.2007; Bardainne et al.2008; Cesca et al.2011). McGarr (1976, 2014) and Klose (2013) have also proposed empirical relations to link the size of the human action (for example in terms of mass shift, volume change or injected volume) with the size of possible induced seismicity. The problem of discrimination among natural, triggered and induced seismicity has recently been debated (Cesca et al.2013a; Dahm et al.2013), following a broad range of questioned anthropogenic induced seismicity cases.

Gas storage systems are typically used to improve the handling of temporal market variations, e.g. injecting gas during periods of low demand and withdrawing it when the demand increases. Different types of underground reservoirs are currently in use for gas storage, including depleted gas and oil reservoirs, aquifers and salt caverns. Injection-induced seismicity (IIS) has been most often discussed for the cases of fluid or wastewater injection, and more recently in relation to hydraulic fracturing operations. Known cases of seismicity induced by gas injection and storage operations are rarer. Microseismicity has been observed accompanying few gas storage operations, e.g. at the Bergermeer gas field, in the Netherlands (Kraaijpoel et al.2012; Kwee 2012), so this sequence offshore Spain would be among the most important study cases of seismicity correlated to gas storage operations. Moreover, different microseismicity rates have been observed while monitoring the exploitation of different gas reservoirs, possibly in consequence of the variable storage type, size, shape and depth, the local seismotectonic setting, the presence and orientation of faults, and the volumes and rates on injected/extracted gas. Benetatos et al. (2013) investigated few microseismic events occurring in proximity to the (cavity) natural gas storage at Haje, Czech Republic, and reported a maximal magnitude of M_w 0.4. The ongoing monitoring of the Collalto, Italy, natural gas storage, where the gas is stored in a depleted gas reservoir (differently from the Castor project case, where an exploited oil reservoir is used), has not detected so far any relevant microseismicity. The Collalto network ensures a completeness magnitude of about 0.0, up to 5–8 km distance to the reservoir (Priolo et al.2013).

In our case, the distribution of seismic stations around the epicentral region (Fig. 1) is quite poor and strongly asymmetric. Stations are operated by different institutions (IGN, IGC and EBR); instrumentation includes broadband, short-period sensors and one broad-band OBS sensor. The closest station (ALCN) has a distance of more than 20 km from the Castor project platform. The inadequate configuration of the monitoring network limits the analysis of microseismicity. Fig. 1 illustrates the network configuration, EBR preliminary locations and the injection platform site, as well as the inferred geometry of the Fosa de Amposta fault system (from García-Mayordomo et al.2012).

This work aims to provide a first, quick analysis of the seismic sequence using a range of seismological techniques. We first analyse the catalogue information, to evaluate the temporal evolution of seismicity rate, the maximum magnitudes and the temporal clustering of seismic events. Then, we focus on the largest events, which mostly occurred after the injection stopped, until 2013 October 4. To tackle the limitation due to the monitoring network geometry, we use a specific full-waveform-based seismological technique and can satisfactorily characterize the largest events in the seismic sequence and compare them to a previous event, taking place prior to the gas storage exploitation. In particular, our work aims to precisely locate the largest events, both by means of absolute and relative locations, and to analyse the main focal mechanisms, through moment tensor inversion. This work provides detailed information on the seismic sources and their locations.

2 TEMPORAL EVOLUTION OF THE 2013 SEPTEMBER–OCTOBER SEISMIC SEQUENCE

The reference catalogue from the Ebro Observatory contains 1002 earthquakes, 982 of them with magnitude assigned (M_L, with three significant digits), from 2013 September 5 to October 15. For the following analysis, we will avoid mixing magnitude scales, using only M_L, because M_w is systematically higher than M_L for the events for which M₀ was calculated. The magnitude–frequency distribution (MFD) of this catalogue was fitted to a Gutenberg–Richter relationship. For this purpose, the entire-magnitude-range method of Woessner & Wiemer (2005) was used, with 5000 bootstraps. This enables an automatic determination of the magnitude of completeness (M_c) and its standard deviation, and a maximum likelihood fit of the b value (based on Aki 1965 and Utsu 1965) and its standard deviation (based on Shi and Bolt 1982). The b value for the whole sequence is not significantly different from unity, b = 0.96 ± 0.05 (Fig. 2). The relatively good instrumental coverage yielded a relatively low M_c = 1.32 ± 0.08. In comparison, the national Spanish earthquake catalogue at this location has M_c ∼ 2 (González, in preparation).

Figure 2.

Cumulative magnitude–frequency distributions. Coloured lines: Observed distributions for the whole sequence, during the injection (phase 1, until 2013 September 16) and after the injection (phase 2, since September 17). Black lines: fits to the Gutenber–Richter relation, indicating the magnitude of completeness (M_c) and the b value of each distribution.

The temporal evolution of the sequence is here discussed. As is commonly observed in other catalogues, earthquakes were predominantly detected during night time (19-6 hr UT, and especially 0-2 hr UT), because of the lower noise level. These strong systematic differences disappear when only the 424 earthquakes with M_L ≥ M_c are considered. The temporal evolution of the seismicity shows an increase of the seismic rate before September 16 and a further sharp increase from September 28 (Fig. 3). We therefore divide the catalogue into two periods, until (phase 1) and after (phase 2) September 16, that is the stop of the gas injection. For these two different periods, we studied the MFD separately as before. The results show that the b value in the first phase (b₁ = 1.41 ± 0.21) is significantly larger than in the second phase (b₂ = 0.87 ± 0.06). In contrast, the magnitude of completeness for these periods (M_c1 = 1.40 ± 0.16 and M_c2 = 1.31 ± 0.08) is consistent with the overall value when their uncertainty ranges are taken into account.

Figure 3.

The temporal evolution of the seismic sequence (first 44 days since the beginning of the injection) is described by means of the daily evolution of number of events with magnitude above completeness (top panel), the maximum daily magnitude (centre) and the increase of the cumulative seismic moment (bottom panel), considering only events above completeness, according to the Ebro Observatory catalogue.

The simplest way to characterize temporal correlations of events in a seismic catalogue is to consider the distribution of time intervals between subsequent events in the catalogue, which are named interevent-times. To describe the tendency of events to be in a temporal cluster, we use the coefficient of variation (CV = δ/μ) parameter, which is the ratio of the standard deviation (δ) over the mean value (μ) of the interevent-time distribution. CV > 1 indicates that the events are clustered and tend to be closer to each other in the short timescale (Kagan & Jackson 1991; Zöller et al.2006; Maghsoudi et al.2014), whereas CV < 1 is typical for swarm-type activity. For the two different phases defined before, only the events with magnitudes greater or equal to above the respective M_c are further considered to compute interevent-times and consequently CV values. The large values derived for both phases before and after September 16 (CV = 2.18 and 2.77, respectively) indicate that the events are similarly clustered in both periods.

Finally, the cumulative seismic moment release of the sequence is plotted in Fig. 3. We consider the M₀ values which we calculated (see below) for the largest earthquakes (which take up most of the total), and for the remaining ones we assume that M_L = M_w, calculating M₀ from the definition of M_w (Hanks & Kanamori 1979). The total cumulative moment release is 1.37 × 10¹⁶ Nm. The choice of the M_L/M_w ratio can partially affect the estimation of the cumulative moment release as well as the overall b value, but not the general pattern of their temporal variation.

The results of this analysis point out some differences before and after September 16. Whereas the magnitude of completeness and the interevent-time distribution are comparable, a significant change is seen in the b value, which decreases from 1.40 to 0.87. This may suggest a change in stresses but also a change in the dominant focal mechanism (Schorlemmer et al.2005). The change in b value is accompanied by an increase of maximum observed magnitudes, which is experienced in the period following September 16. In the following, we will analyse in detail the largest events in the sequence. However, given the previous considerations, the derived source parameters may remain indicative only for the second phase of the sequence.

3 ABSOLUTE AND RELATIVE LOCATION

We located the seismic events using the waveform stacking approach proposed by Grigoli et al. (2013a,b). In this procedure, the space of possible locations is scanned and the STA/LTA traces of selected characteristic functions are stacked along the P and S traveltime surfaces corresponding to the selected hypocentre. We compute the STA/LTA of the vertical energy trace, which enhance the signal of P wave, and STA/LTA of the squared maximal eigenvalue of the instantaneous covariance matrix obtained from horizontal components traces, which is more sensitive to S wave (for further details, see Grigoli et al.2013b). Iterating this procedure on a 3-D grid we retrieve a multidimensional coherence matrix, whose absolute maximum corresponds to the coordinates of the seismic event. We performed event location by direct search within a grid with size 175 × 175 × 40 km³ and with a spacing of 500 m. The P and S traveltimes were computed by using the eikonal solver developed by Podvin & Lecomte (1991). Since at least three different P-wave velocity models have been proposed for the study area, we tested those to evaluate their performance in locating the largest events. Two models were based on reflection seismic profiles (Torné et al.1992; Vidal et al.1998) and one was chosen from the CRUST 2.0 database (Bassin et al.2000). The best results (i.e. higher coherence values obtained during the location process) have been obtained using the average velocity model extracted from the CRUST 2.0 database corresponding to the Amposta platform location. We then applied the Wadati method to estimate an S-wave velocity model, under the assumption of a constant Vp/Vs ratio. For this dataset we retrieved a Vp/Vs of 1.69. Absolute locations are then finally performed using the CRUST2.0 model, and with a fixed Vp/Vs of 1.69. The location is performed for 73 events, including all events with magnitude above 2 (mbLg magnitude estimations from the IGN catalogue), plus a few additional smaller earthquakes (e.g. following the main shock on October 1), and also including the previous event of 2012 April 8 (mbLg magnitude estimations from the IGN catalogue). For smaller earthquakes the seismic signals are strongly contaminated by noise and main onsets only seen at the closest stations, posing a strong limitation to the location resolution, which are not considered here. Fig. 4 shows the coherence matrices related to the strongest studied event, which occurred on 2013 October 1. The position of the coherence matrices maximum corresponds to the hypocentral location of the considered seismic event. Fig. 5 shows the location of all events of the dataset. Strong variations are determined, when using or excluding the OBS station in the location procedure. The exclusion of station COBS improves the maximal coherences, and results in more spatially clustered locations. Such a result suggests a timing problem or anomalous performance of the OBS station. The figure highlights some important findings. First, the seismic event of 2012 April 8, which took place prior to the gas injection operations and is here considered as a reference for typical natural seismicity in the region, is located about 10 km away from the platform and the region struck by the recent seismic sequence. Its proximity to a mapped fault suggests its relation with that seismogenic structure. A second relevant result is that the location method is able to reduce the spatial scattering of epicentral locations, which affected other catalogues. All studied events of the new sequence cluster spatially and are located within a very limited region, with a maximal horizontal extension of less than 5 km. All events are very shallow, with depths between 0 and 3 km. It is also evident that the epicentral region is consistent with the location of the injection, a result which is very different from previous location results (NW–SE distributed epicentres, with largest events located at about 5–6 km SE of the platform). Locations uncertainties are calculated using a bootstrap approach, iteratively perturbing the STA/LTA parameters (i.e. the length of the short and the long-time window). After each perturbation an event is then located again (up to 20 times). For all events we retrieved hypocentral uncertainties within the range 0.5–2.5 km. These values confirm a localized seismogenic region, located very close to the injection point and at comparable shallow depths. Note that these uncertainties may be underestimated as they only account for the effects of the STA/LTA parametrization, but not, for example the choice of a simplified velocity model.

Figure 4.

Absolute location result for the 2013 October 1 M_w 4.3 earthquake, the largest of the seismic sequence. The three plots illustrate the horizontal and vertical (north–south, east–west) projections of the 3-D coherence matrix, according to the given colour scale. The preferred location is identified at a very shallow (<2 km) depth.

Figure 5.

Absolute location results for the largest events in the sequence, and faults in the epicentral area, including (a) or excluding (b) the OBS station COBS. Circles: epicentral locations, with colour corresponding to the timescale on the right side; the prior 2012 April 8 event is also shown. White square: location of the Castor platform. Coloured lines: faults in the proximity of the injection site (faults further away are omitted). In the largest map, the rough locations of the Amposta fault strands (red lines) and a different mapped fault (blue line) are plotted according to García-Mayordomo et al. (2012). The small panel shows a more detailed view of the epicentral region: more detailed digitalized faults are shown (green lines), according to Geostock (2010). These include the Eastern Amposta fault, striking NNE–SSW below the platform, different steep subparallel faults striking NW–SE on the NW side of the Amposta fault, and few subfaults on the opposite side of the Amposta fault striking NE–SW but with different dip angles.

Absolute locations provide clear evidence that seismicity occurred in the vicinity of the injection wells, but the resolution is not sufficient to infer the geometry of the hosting fault. We use a relative location technique to further investigate the spatial distribution of hypocentres of the largest events at a small scale. This method is based on waveform correlation, and could be adopted because the events show similar waveforms. Velocity waveforms are bandpass filtered between 1 and 4 Hz, and tapered to consider time windows of 4 s length, centred at the arrival of the first P and S onsets. Only similar waveform pairs are used for the relocation procedure, fixing the threshold for waveform similarity at a 0.6 correlation coefficient. The weighting of waveform pairs in the inversion procedure scales linearly with the correlation coefficient (a test using a weighted scheme upon the square of correlation coefficients produced almost identical results). The velocity model used for the synthetic traveltime is the modified CRUST2.0 model adopted for the absolute location. The inversion for the relative location is performed by solving the following equation:  

\begin{equation} D_{\it kij} = (x_j - x_i )s_k + (t_j - t_i ), \end{equation}

(1)

where D_kij is the observed differential time between P phase arrival times of event i and j at station k, s_k is the slowness take-off vector at master event (for instance event j) to station k (at 2 km depth), x_i is the spatial coordinates of event i, and t is the origin time of event i. Results (Fig. 6) were obtained for a subset of 51 earthquakes, after automatically removing those showing poor average correlation coefficients with the others. In the plot, a possible distribution of absolute epicentral locations is shown, which was chosen to have the average location consistent with the platform location. However, it is important to remind that the performed analysis can only provide relative locations. These show mainly that epicentres are distributed in a small region, elongated NNE–SSW, with a maximal extension of about 5 km. This result confirms that epicentres are confined to a small region, which location close to the platform is supported by the absolute locations. The limited extension of hypocentres in depth confirms the absolute location findings, in terms of a very narrow depth range at about 0–4 km. Uncertainties are estimated by repeatedly solving the linear inverse problem for perturbed time lags; the amount of perturbation has been estimated from the variances between modelled and observed time lags. The results indicate that the uncertainties are below 2 km in depth, and 400 m for the epicentral location. A rough spatio-temporal evolution of the seismic sequence can be seen. The earliest earthquakes were located at the northern side of the affected region, followed by others to the south after injection stopped and during 2013 September, and finally the latest ones originated again at the northern side in 2013 October, when the seismicity is more spatially distributed.

Figure 6.

Relative location results, as map view and vertical cross-section. Epicentres and faults are plotted according to conventions in Fig. 5, for the region defined in the small panel. Note that the absolute locations are not known and were chosen to be centred at the Castor platform location, and only the relative location can be discussed (i.e. they could either fit a distribution along the Eastern Amposta fault, along the subparallel faults on its SE edge or the activation of different faults on its NW edge).

4 SEISMIC SOURCE INVERSION

Once defined an improved P- and S-velocity model, and relying on previous absolute hypocentral locations, we focus on the analysis of the focal mechanisms. Moment tensor inversion is performed using the Kiwi tools (Cesca et al.2010; Heimann 2011), following the procedure described in Cesca et al. (2013b). The inversion is carried out in two steps. In the first step we compare observed and synthetic amplitude spectra of the whole waveforms to derive the best fitting focal planes, using a pure double couple (DC) and a full moment tensor (MT) point source model. In the second step, we compare full waveforms in the time domain to define the focal mechanism polarity and to obtain the centroid location and time. A similar approach was previously successfully adopted to invert moment tensors at different scales, including natural seismicity with a very poor azimuthal coverage (Custodio et al.2012; Domingues et al.2012), the analysis of a mixed natural/induced seismicity dataset at regional distances (Cesca et al.2013b), and the investigation of mining induced microseismicity sources at a local scale (Sen et al.2013; Cesca et al.2014). Synthetic waveforms and spectra are built using a single 1-D velocity model, upon the CRUST 2.0 database and 1.69 V_P/V_S velocity ratio discussed in the previous section. Synthetic seismograms include near-field terms. Displacement waveforms are filtered between 0.05 and 0.10 Hz, for events with mbLg > 3.5, and between 0.067 and 0.125 Hz when mbLg < = 3.5. These conditions reduce the available dataset to broad-band stations only (the OBS station is also excluded, because the 1-D model cannot account for the water layer). We further limit the used data to stations located at less than 200 km from each epicentre, which show the best quality signals, and remove single traces for specific earthquakes, whenever the signals are contaminated by strong seismic noise. The second inversion step is limited to stations with less than 150 km epicentral distance.

The DC inversion results are illustrated in Fig. 7, whose bottom panel shows the spectral and waveform fit for the largest event in the series (2013 October 1, M_w 4.3). The spectral match is good (L² norm misfits below 0.35 for all events with M_w larger than 3.3) for all studied events in the sequence, and support the quality of the derived focal mechanisms. The waveform fit is also very good at the closest stations; at further distance, at stations EIBI and ETOS located on Balearic Islands, we observe a time offset between data and synthetics, suggesting a velocity anomaly with respect to the adopted velocity model. The centroid locations are always found very close to the assumed hypocentral locations, and thus confirm the locations results obtained in the previous section. Centroid depths are inverted independently with the moment tensor method. Centroids are found at very shallow depths, typically of 2 km, a result consistent with location findings. The uncertainty on the depth estimation has been assessed through a bootstrap test (see Cesca et al.2013b), whose results (Fig. 7) confirm the finding of very shallow sources.

Figure 7.

Moment tensor inversion results. Top: source parameters (strike, dip, rake, best depth value and its range of estimate according to the bootstrap test, scalar moment and magnitude), best and mean pure DC focal mechanisms (black and blue focal spheres, respectively) as obtained in this study. DC focal mechanisms by IGN (grey focal spheres) have been included, whenever available from the IGN webpage. Bottom: comparison of amplitude spectra and waveform fit after different inversion steps, for the case of the 2013 October 1 M_w 4.3 earthquake (red is used for observations, black for synthetics).

The largest events in the sequence repetitively occur with a similar oblique focal mechanism, which suggest the activation of the same fault or a system of parallel faults. The possible orientation of the hosting fault (or faults) is either NNE–SSW, dipping with a low angle (15–42°) towards SE, or NW–SE, with a steep plane likely dipping towards SW. The robustness of focal mechanism solutions has been investigated using a jack-knife approach, repeating the inversions after excluding each station. The focal mechanism distributions confirm previous mechanisms, and indicate that the less resolved fault plane parameters are the dip angle of the NNE–SSW oriented plane, and the rake along the NW–SE plane. This effect is likely due to the combination of focal mechanism and source–receivers geometry. Single focal mechanisms range from oblique to strike-slip mechanisms, the latter associated with a steeper NNE–SSW plane, still confirming a dipping towards SE. In Fig. 7, in addition to the best solutions (black focal spheres), obtained using all stations, we also report the mean mechanisms (blue focal sphere). These mean mechanisms have been obtained from the distribution of focal mechanism results of the jack-knife analysis. They are determined by a direct search of the focal mechanism which minimize the differences from all mechanisms in the distribution, using the approach described in Cesca et al. (2014) and the Kagan angle (Kagan 1991, 1992) to measure the distance between DC solutions. Non-DC components are typically low, and the misfit improvement from the DC to the full MT solution is minor (e.g. misfits of 0.327 and 0.324 respectively for DC and full MT models for the largest event). These results exclude the presence of strong tensile components or volumetric changes, which could suggest, for example the occurrence of collapse processes (Cesca et al.2013b). Our focal mechanisms are in general agreement with solutions proposed by IGN for few events (Fig. 8), although our solutions are more homogeneous among them, and have lower dip angles for the NNE–SSW plane (further referred as plane 1, with strike 22–49°, dip 15–42°, rake –10° to 0°), and a steeper one for the NW–SE plane (plane 2, with strike 119–140°, dip 83–90°, rake –131° to 105°). As a result of these minor discrepancies, IGN solutions tend to have larger strike-slip components. Minor differences to IGN solutions may be explained in terms of the moment tensor inversion methodology as well as data used for the inversion. For example, for the 2013 October 1, we have used six stations, with azimuthal coverage on the first, second and fourth quadrant, whereas the IGN solution is based on three station with an azimuthal gap of 270°. Interestingly, when excluding the closest broad-band station and only using IGN stations, the moment tensor inversion provide focal mechanisms with steeper NNE–SSW planes of about 60° dip, fitting well with reference IGN solutions.

Figure 8.

Sketches of the possible rupture scenarios, along (a) a EW section (after the Shell Spain seismic profile in Seeman et al.1990, with vertical scale in two-way traveltime) and (b) a horizontal projection, at the Castor platform. The sketch includes the rough location of the reservoir (yellow region) and rough depth at which the map view is plotted (dashed line). The blue and green lines are stratigraphic markers, with the blue being the contact between the reservoir rocks (karstified Mesozoic carbonates) and the impermeable Miocene overlying it, and the green being the Messinian unconformity, which is the level at which the faults in Fig. 8(b) are drawn. The first possible scenario involves the rupture of a NE–SW striking low-angle fault, dipping towards SE, which could match the fault identified by Seeman et al. (1990), as shown in the cross-section sketch. The alternative scenario requires the activation of different subparallel faults known as Montsia faults (Geostock 2010), which are better seen in the map view panel. The scenario of the activation of the main fault (Eastern Amposta fault) is excluded by the joint interpretation of earthquake locations and focal mechanisms.

The striking of plane 1 is well consistent with the average orientation of the Eastern strand of the Amposta fault system, but the dip is not (the plane dips towards the open sea whereas the Eastern Amposta fault dips towards the Spanish coast), and actually makes this plane perpendicular to its known fault geometry. Instead, the low-dip plane 1 results parallel to the sedimentary stratification above and below the reservoir, and could fit a small fault recognized by Seeman et al. (1990). The alternative plane 2, which is oriented NW–SE, is consistent in strike and dip with a system of small subparallel faults located on the NW side of the Eastern Amposta fault (Geostock 2010). Finally, the inversion result for the M_w 3.5 event occurring on 2012 April 8, shows a normal fault mechanism with different orientation. However, the reliability of this solution is limited, because of the few reliable traces and the large azimuthal gap (260°). There are no focal mechanism solutions available in the study area, for events occurring previously to this seismic crisis. However, there are few solutions for earthquakes occurring offshore towards NE. The nearest solutions to the Castor zone correspond to strike-slip and reverse solutions, both with horizontal pressure axis changing form N–S to NE–SW direction (Olivera et al.1992; Buforn & Udías 2003).

5 DISCUSSION

Can the rupture geometry associated to the largest events of the sequence be inferred from hypocentral locations and focal mechanisms? This work suggests that the question can be partially positively answered. Three scenarios are foreseen: (1) the activation of the Eastern Amposta fault, (2) the activation of a small fault with strike parallel to the latter, but dipping towards SE, which could fit to the seismic interpretation by Shell Spain (Seeman et al.1990) and (3) the repeated activation of one or more faults perpendicular to the Eastern Amposta fault and situated on its NW side, known as Montsia faults, which were identified by Geostock (2010). Location results, both using absolute and relative methods, suggest a spatial distribution of epicentres extended along a SSW–NNE direction. This finding supports the first two hypotheses. However, given the uncertainties of absolute locations, the third case cannot be excluded, if we assume that the seismicity is distributed along multiple parallel faults and not a single one. The analysis of focal mechanisms identifies two possible orientations for the rupture plane, one steep plane oriented NW–SE and one low-angle dip plane oriented NNE–SSW, dipping towards SE. This result excludes the case of the Eastern Amposta fault, which dips in the opposite direction towards the Spanish coast, with a dip angle inconsistent with focal mechanism solutions. The remaining two possible scenarios are illustrated in the Fig. 8. The model predicting the activation of a low-angle fault dipping SE may be problematic, as the only fault with a consistent orientation at the Castor location seems only about 1 km in depth in extent and is possibly too small to host a M_w 4.3 event (which have, e.g. 3–5 km rupture length for shallow earthquakes, Dahm et al.2007; Cesca et al.2011). An alternative hypothesis according to this scenario would be by means of interlayer slip, with the slip occurring along a bedding plane corresponding to a stratigraphic interface; the orientation of these interfaces at the injection point is consistent with the orientation of the fault plane of scenario 1, and bedding planes have a much larger surface area, which could accommodate the observed earthquakes with magnitudes above M_w 4. On the other hand, the alternative scenario 2 requires the reactivation of a complex system of parallel faults (Montsia faults). With the current results, we cannot fully discriminate among these two rupture scenarios but the location results would favour the first one. All results are consistent in finding a spatially limited seismogenic region, not exceeding a length of 5 km, and being limited to the first 3–4 km of depth.

Given the lack of significant seismic source studies in the region, it is difficult to judge how the September/October events are different, in terms of locations and focal mechanisms, from previous seismicity. We could show that the seismic sequence is highly localized in space, and very repetitive in terms of focal mechanisms, and that both locations and focal mechanism differ from the reference event of 2012 April 8. Focal mechanisms also differ from those previously known in the area, although only few focal mechanisms solutions are available, and at locations far from the epicentral region of the new seismic sequence.

Should the earthquakes be considered as cases of natural, triggered or induced seismicity? The spatiotemporal correlation between the last injection test and the seismic sequence, together with the size and rate of seismicity, seem compelling arguments to infer a correlation between the injection and the seismicity. All events occurred in a similar depth as the injection point, within only few kilometres or less around it. The first event occurred only few days after the beginning of the injection. After injection shut-in, the frequency magnitude distribution changed from b values of about 1.4 to about 0.8.

Different mechanisms of reactivation of existing faults due to fluid injection or fluid withdrawal have been proposed (e.g. McGarr et al.2002; Ellsworth 2013). For instance, oil or gas withdrawal from a porous reservoir formation, which is sealed from bottom and top, leads to a depletion of the reservoir layer and induces stress perturbations in the surrounding rock mass. Faults in the region of the perturbed stress field may be re-activated (model 1). This type of induced and triggered seismicity has been observed for conventional gas and oil fields during and after production, and might be a candidate mechanism for the Castor oil field exploited between 1973 and 1989. However, the reservoir where the injection of the Castor project took place is a karstic system with an active water drive, and is sealed only from top (Seeman et al.1990). Therefore, stress perturbation during injection and extraction should be small and transient. Moreover, model 1 does not explain the temporal correlation to the injection operation in 2013. Another mechanism for induced seismicity explains the instantaneous failure of faults by means of pore pressure changes, or lubrication, on the fault itself (model 2). Pore pressure increase will increase the effective Coulomb stresses acting on the fault, and will bring the fault closer to failure. Faults under critical pre-stress may rupture. The mechanism requires a hydraulic connection between the fluid injection point and the affected faults (e.g. McGarr 2014). Typically, this type of seismicity controlled by pressure diffusion begins in close proximity to the injection point and migrates with time to larger distances. Pressure diffusion is not very fast, and the triggering of earthquakes several kilometres from the injection point after only few hours is unusual. Also, an outward migration of seismic events was not observed in this case. On the other hand, as discussed before, the locations uncertainties are possibly too large to resolve a small-scale migration. Interesting for model 2 is that the pore pressure change does not affect the shear stress acting on the fault. This means that rupture, triggered by pore pressure change, is expected to slip in the direction of the pre-existing stresss (e.g. the regional tectonic stress). A third mechanism of fluid-induced seismicity is given by the formation of hydrofractures, that is growing, fluid-filled tensile cracks. Hydrofracture experiments are able to induce micro-earthquakes. If hydrofractures are growing large, the magnitudes of the events may become larger. For the Ekofisk oil field 2002 M_w 4.3 shallow event we have proposed that a hydrofracture triggered an earthquake rupture on a subhorizontal plane at the border of the oil reservoir, where the depletion-induced stresses were large (e.g. Cesca et al.2011). A similar mechanism may be suggested for the M_w 4.3 earthquake in this study. However, the formation of a hydrofracture needs an injection overpressure large enough to overcome the tensile strength of the sediments. Hydraulic fracturing was unlikely in this case, because during the injection test care was taken not to disrupt the reservoir sealing, so the injection overpressure was kept relatively small, not exceeding 0.8 MPa, according to the company (IIE 2013). This is more than 20 times smaller than the hydrostatic pressure at the injection depth. Given the geological information on the existing local faults, the source mechanism may give hints whether the fault and rupture have been favourable oriented with respect to the regional tectonic stress, and to possibly distinguish model 1 from model 2.

The world stress map project (WSM, Heidbach et al.2008) provides some information on the regional stress, indicating that the direction of maximal compression (SHmax) is about 10° NNE. Unfortunately, the WSM does not provide clear information on the stress regime at the location of the earthquakes, and both strike slip and normal faulting regimes are indicated at some distance to the study area, with a slight predominance to strike-slip cases. A similar variability of fault regimes was confirmed by Schindler et al. (1998), who reported a strike slip/normal faulting regime and a maximal horizontal compression (Shmax) direction varying from 8° to 36° NNE at about 10 km distance from Castor (wells Delta C-3, Delta E-3 and San Carlos III-I). We assume the strike φ of SHmax of 23° +−14° striking, and a magnitude of the vertical principal stress of similar to SHmax. Fig. 9 shows the relative size of the Coulomb stress acting on pre-existing faults with different strike and dip angles, assuming a frictional angle of 0.2. Additionally, the direction of shear stress is indicated in terms of a rake angle. The observed rake angle should be fit to the expected shear stress direction on the fault, if the rupture was driven by the background stress (model 2). The comparison with the source mechanism of the largest event shows, that plane 1 with strike 37° (strike relative to SHmax is 37°–23° = 14°), dip 30° and rake –3° was more favourable oriented and ruptured in more favourable direction of the pre-existing stress than plane 2 with strike 129°, dip 89° and rake –129°. If the earthquakes were triggered by pore pressure increase, the observed slip direction on plane 1 was more likely to occur. It is also likely the shallow SE dipping fault in Fig. 8 is connected to the reservoir layer.

Figure 9.

Relative positive Coulomb stress (size of circles) and direction of shear stress (colours, see colour bar) calculated for pre-existing faults with different strike and dip angles. A friction coefficient of 0.5 was assumed. The assumed regional stress field has a maximal compressive stress striking 10° NE, a least compressive stress striking 100° and an intermediate vertical stress of 80 per cent of the least compressive stress. The strike-dip-rake points of the two nodal plans of the largest studied events are indicated by coloured triangles.

Upon these results it seems realistic that the gas injection can have triggered seismicity by reactivating the shallowest part of the SE dipping fault (fault 1) or similarly oriented bedding planes. However, our approach cannot formally exclude the possibility that a large stress perturbation is still existing from the former exploitation of the oil reservoir, nor the hypothesis that the earthquakes were purely natural, although less likely. This would require further detailed information, not available at the moment: a detailed knowledge of the geological structure, information on the gas injection rate, further information on the previous oil exploitation of the reservoir, and a complete modelling of stress perturbation.

It has recently been proposed that the maximum seismic moment release, and the maximum magnitude, in an earthquake sequence triggered by the pore pressure increase due to fluid injection is proportional to total fluid volume injected (McGarr 2014). This model assumes that the fluid is injected in a fluid-saturated rock mass, with faults well oriented for slip in the ambient stress field and stressed to within a seismic stress drop of failure, and with a friction coefficient with the typically observed value of 0.6. This model does not apply to the Amposta field, because this is a karstic formation in drained conditions, where water can easily flow upon the application of an overpressure and thus the overpressure should be smaller than assumed in this model. We can check the resulting estimates as follows. The company reported a pore pressure in the reservoir after injection of ≈18 MPa (IIE 2013), with a static overpressure of 0.6 MPa. This is similar to the static pressure of 17.6 MPa (hydrostatic plus atmospheric) at the injection depth, considering the density of sea water (≈1.025 kg m⁻³, e.g. Talley et al.2011). The temperature at the location and depth of injection is ≈80 °C (Fernández et al.1990, 1998). Assuming thermal equilibrium of the gas, its volume in the reservoir can be calculated considering the ideal gas equation (e.g. Lyons & Plisga 2005). The reservoir rocks are almost pure, dense limestone with very low porosity (<2 per cent), according to the sonic logs (Batchelor et al.2007), consistent with the shear modulus of 30 GPa (Sayers 2008) used as a generic value by McGarr (2014). The maximum cumulative moment release, according to McGarr (2014) would then be ∑M₀ ≈ 2.01 × 10¹⁶, and the maximum magnitude of M_w ≈ 4.8. Both the maximum magnitude observed in the sequence (M_w = 4.3) and the cumulative seismic moment (∑M₀ ≈ 1.37 × 10¹⁶ Nm) are well below this upper bound. Finally, we also have to note a potential caveat of estimating maximum magnitude bounds. For triggered seismicity the induced stress perturbation is only responsible for the nucleation process and the rupture is driven by the tectonic background stress. Thus, the maximum magnitude may be only constrained by the size of the activated fault.

6 CONCLUSIONS

The seismic sequence which struck offshore Spain, in the Gulf of Valencia, is an interesting case of seismicity correlated to gas injection operations. The seismicity reached a magnitude M_w 4.3, and more than 1000 events were recorded during a period of about 40 d, starting on 2013 September 5. The seismic sequence shows a temporal variation, correlated with the beginning and the end of the injection process. Earthquake activity started with the beginning of fluid injection, and b values changed from the co-injection to the post-injection phase. Both indicate that the events could have been triggered by pore pressure changes on pre-existing faults.

The largest events mostly occurred in the post-injection phase, which is a common observation in other cases of seismicity related to fluid injection. We successfully applied full-waveform techniques to determine absolute and relative locations and to perform a moment tensor inversion for the larger events in the post-injection phase, despite of the poor network configuration. The results indicate that seismicity is confined in a very small region of less than 5 km size, in proximity to the gas injection wells. The combined interpretation of seismological analysis results (location, focal mechanisms) and local seismic surveys and small-scale fault mapping suggest two possible rupture scenarios, but exclude the reactivation of the largest fault in the vicinity, the Eastern Amposta Fault. The two possible cases either involve a low-dip failure plane striking roughly parallel to the Eastern Amposta Fault but dipping perpendicular to it, or a system of subvertical faults oriented NW–SE. Both cases would be consistent with small mapped faults interconnected to the reservoir layer. However, contrary to the rupture on the subhorizontal plane which is consistent to the trigger model by pore pressure change, rupture on the subvertical planes could only be explained if a large stress perturbation would be present, since these planes are not favourably oriented to the regional stress field and rupture direction differs from the resolved shear stress direction.

The total seismic moment of the sequence was 1.37 × 10¹⁶ Nm and the observed maximum magnitude of M_w = 4.3. At the time of this writing, further injection operations are halted. It is difficult to foresee, whether the activity would continue if the gas injection resumes in the future, and what the maximum magnitudes could be. It would depend on the available stress on pre-existing faults and the size of the subfaults hydraulically connected to the reservoir formation.

We cannot completely exclude a natural cause for the earthquakes, although this seems unlikely, since magnitudes and seismicity rates have no precedent in the region in the last two decades. However, earthquakes of similar magnitude could have occurred in the past, without being felt and/or reported in historical catalogues, given their moderate magnitudes and the offshore locations. Future monitoring of the area, possibly with a denser network and at closer distance to the epicentral region could be helpful to confirm our findings through the monitoring of microseismicity.

We are thankful to Dr J. Wassermann and two anonymous reviewers for their comments. This work has been funded by the German BMBF ‘Geotechnologien’ project MINE (Grant of project BMBF03G0737) and Ministerio de Ciencia e Innovación project ALERT-ES (CGL2010-19803C03-01). ÁG was supported by the Caja Madrid Foundation (Spain) and the Spanish Government grant FIS2010-19773. We are thankful to the Seismic Network of the National Geographical Institute of Spain (IGN), to the Catalan Seismic Network of the Geological Institute of Catalunya (IGC) and to the Ebro Observatory for providing seismic waveform and catalogue data.