Earthquake static stress transfer in the 2013 Gulf of Valencia (Spain) seismic sequence

On 24 September 2013, an Ml 3.6 earthquake struck in the Gulf of Valencia (Spain) near the Mediterranean coast of Castelló, roughly 1 week after gas injections conducted in the area to develop underground gas storage had been halted. The event, felt by the nearby population, led to a sequence build-up of felt events which reached a maximum of Ml 4.3 on 2 October. Here, we study the role of static stress transfer as an earthquake-triggering mechanism during the main phase of the sequence, as expressed by the eight felt events. By means of the Coulomb failure function, cumulative static stress changes are quantified on fault planes derived from focal mechanism solutions (which act as both source and receiver faults) and on the previously mapped structures in the area (acting only as stress receivers in our modeling). Results suggest that static stress transfer played a destabilizing role and point towards an SE-dipping structure underlying the reservoir (or various with analogous geometry) that was most likely activated during the sequence. One of the previously mapped faults could be geometrically compatible, yet our study supports deeper sources. Based on this approach, the influence of the main events in the occurrence of future and potentially damaging earthquakes in the area would not be significant.

Abstract.On 24 September 2013, an M l 3.6 earthquake struck in the Gulf of Valencia (Spain) near the Mediterranean coast of Castelló, roughly 1 week after gas injections conducted in the area to develop underground gas storage had been halted.The event, felt by the nearby population, led to a sequence build-up of felt events which reached a maximum of M l 4.3 on 2 October.
Here, we study the role of static stress transfer as an earthquake-triggering mechanism during the main phase of the sequence, as expressed by the eight felt events.By means of the Coulomb failure function, cumulative static stress changes are quantified on fault planes derived from focal mechanism solutions (which act as both source and receiver faults) and on the previously mapped structures in the area (acting only as stress receivers in our modeling).Results suggest that static stress transfer played a destabilizing role and point towards an SE-dipping structure underlying the reservoir (or various with analogous geometry) that was most likely activated during the sequence.One of the previously mapped faults could be geometrically compatible, yet our study supports deeper sources.Based on this approach, the influence of the main events in the occurrence of future and potentially damaging earthquakes in the area would not be significant.

Introduction
With an increasing demand for energy resources, the Spanish government launched the Castor Underground Gas Storage (UGS) project in 1996 (Fig. 1a).Its development should have provided more autonomy to the Spanish gas system, which is highly dependent on incoming gas from northern Africa and Europe.The geological structure selected to store gas (Fig. 1b) was the depleted Amposta oil field, which had been operated by Shell in the period 1970-1990 by exploiting its naturally contained heavy oil of 17 • API (Seeman et al., 1990;Batchelor et al., 2007;Escal, 2014).
Natural seismicity in the area is very low (Instituto Geológico y Minero de España, IGME, 2015a, b) despite the fact that several faulting structures had been described or inferred by previous studies in the surroundings of the UGS, mainly using geological approaches (Fonboté et al., 1990;Roca and Guimerà, 1992;Gallart et al., 1995;Vergés and Sàbat, 1999;Perea, 2006;IGME, 2015a, b).Escal, the operating society of the UGS, explored a 3-D cube by using seismic reflection profiles, which enhanced the existing knowledge of the geological structure surrounding the reservoir (Geostock, 2010).The reservoir itself is located roughly 1.7 to 2.5 km beneath the seabed.The fault system derived from Geostock (2010), as in Cesca et al. (2014), is considered here (Fig. 1c).Throughout the paper, we refer to these faults as the "previously mapped faults" or the "mapped faults" and model them along with the faults derived from focal mechanism solutions.
During the third injection phase of cushion gas, carried out from 5 to 17 September 2013, seismic activity built up to more than 30 events per day (Fig. 2).Although the immediate response after injections were halted on 17 September was a decrease in seismic activity, on 24 September the first felt event took place (M l 3.6), and on 1 and 2 October, while nearly 100 events were being recorded each day, 3 events of M l > 4 occurred, raising public concern.The last day with felt events was 4 October, and seismicity returned to basal Figure 2. Histogram of Castor's seismic sequence showing the number of earthquakes per day and their magnitudes.Two separate phases can be distinguished in this sequence.The first lasted until 19 September, just 2 days after injections were stopped, and maximum magnitudes did not surpass M 3.After 4 days of almost no seismicity, the first felt earthquake took place on 24 September (M 3.7), and high levels of seismicity with three M 4 earthquakes were recorded during the 2 following weeks.Modified from ICGC (2013).levels at the end of the month.After the earthquakes felt by the local population, the project was paused and is at the time of the study in a hibernation phase.The Spanish government has to decide whether the facility should be dismantled or storage operations can be resumed (e.g., González, 2014).
Nonconfidential literature in relation to the Castor case has already addressed problems in earthquake locations and focal mechanism computation, as well as the frequencymagnitude relationship (Institut Cartogràfic i Geològic de Catalunya, ICGC, 2013;IGME, 2013;Instituto Geografico Nacional, IGN, 2013;Cesca et al., 2014;González, 2014;Gaite et al., 2016).Cesca et al. (2014) included some discussion in relation to the mapped faults that would have been more likely to slip based on background stress.However, quantification of the physical mechanisms that may generate earthquakes and their evolution in the sequence had not been done beforehand.
Fluid injection-induced earthquakes are typically initiated close to the injection point and then expand outwards under the assumption that fluid pressure diffusion controls the process.This is a common behavior often referred to as the Kaiser effect (e.g., Baisch et al., 2009), which can be explained by the stress memory of the system (i.e., slip is only produced in regions where maximum excess pore pressures are exceeded) and requires hydraulic connection between the injection area and the fractures (e.g., McGarr, 2014).The Castor UGS was placed in a karstic reservoir sealed only at the top with the faults limiting the lateral ends of the reservoir acting as structural traps (Seeman et al., 1990).This means that hydraulic junction with them should be minimum, and could support downward migration as the preferential fluid path.If that had been the case, fluid pressure diffusion, which takes time, could have triggered events on underlying faults after some days (the earthquakes under analysis in this paper took place between 1 and 2 weeks after the well shut-in).However, it should be pointed out that no evidence of outward migration of seismicity has been found (Cesca et al., 2014;Gaite et al., 2016), so this mechanism is somehow discredited.In addition, when a pressurized fluid is injected into a rock matrix, there is an instantaneous stress generation on the whole volume due to mechanical deformation of the matrix itself.Furthermore, coseismic stress transfer as studied herein can trigger events at a distance (e.g., King et al., 1994).Static earthquake triggering can be understood as the result of a fault slip (introducing strain) that, accounting for a confined medium, translates into a stress perturbation that can destabilize other faults.Although other physical mechanisms have been described for induced seismicity (e.g., Ellsworth, 2013), they were most likely less significant for the timescale that is under analysis here.
The performed modeling in this study considers static stress redistribution on its own, and therefore no quantitative account for fluid flow is given herein.Although limited in nature, this is a possible approach when the focus is exclusively on static stress redistribution (Baisch et al., 2009;Schoenball et al., 2012;Catalli et al., 2013).Specifically, we aim to determine, by means of a quantitative assessment, whether the aforementioned mechanism was a destabilizing factor during the strongest phase of the series (which started on 24 September; refer to Fig. 2).It should be noted that the seismic series near the Castor UGS spanned a period of roughly 2 months and started at the beginning of September 2013.While the destabilizing role of static stress transfer concerning the main events is under analysis here, the origin of the series cannot be explained by our modeling and is thus beyond the scope of the paper.
As a result of the latest earthquake locations (ICGC, 2015; Gaite et al., 2016) and obtained focal mechanism (FM) solutions in this work, experienced seismicity and the previously mapped structures cannot be linked in a straightforward way.Therefore, we compute cumulative Coulomb stress changes ( CS) on both the previously mapped and FM-derived faults as a result of slip on the FM faults only.Because these FMderived faults are earthquake sources but also receive stress along the sequence, the analysis allows for the assessment of whether each event could have been caused by the stress perturbation introduced solely by the previous tremors.We also evaluate the evolution and final stress state in terms of CS on the mapped faults and consider whether the experienced events would have shortened the occurrence of a characteristic earthquake in the Main Fault.

Data and resources
FM solutions for the eight strongest events in the sequence had already been found in previous literature (see Table A1, Appendix A).Solutions by both Frontera et al. (2013) and IGN (2013) were obtained shortly after the occurrences with unrevised locations.The Cesca et al. (2014) solutions were found to be shallow and mainly supported the (repetitive) rupture on a fault plane approximately striking NNE-SSW and dipping gently to the SE.However, later findings by Gaite et al. (2016), who used a new 3-D velocity model for the region, placed most of the events deeper underground (around 6 km).This later finding accords with reported depths by ICGC (2015).Owing to the fact that earthquake location is a necessary parameter for FM computation, and thus influences the result, FMs are again computed here for the main events in the series.
We use the Delouis (2014) method, FMNEAR, which uses waveform modeling to determine FM solutions (Table A1; summary provided in Table 1) via its online service (FMNEAR webservice, Delouis, Gerakis, Deschamps, Geoazur/Observatoire de la Côte d'Azur; Delouis et al., 2016).For FM computation, we select the input locations in ICGC (2015) and the recorded waveforms at stations shown in Fig. 1a (see ICGC, 2000).As an example, Fig. 3 shows the FM solution and waveform fit at stations ERTA and CMAS for the M 3.9 event of 30 September obtained by applying  FMNEAR.The velocity model selected to perform calculations corresponds to the default model used in the FMNEAR webservice unless otherwise specified by the user.The stability of the computation proved positive for all eight M l ≥ 3.5 events, but none below that threshold value could be computed with sufficient reliability.FM depth (see Table A1) is in good agreement with location results reported in Gaite et al. (2016).Throughout the analysis, we keep an eye on hypocentral locations reported by Cesca et al. (2014) and Gaite et al. (2016) in addition to those in ICGC (2015); when it comes to the discussion of the results as a whole, the distribution of the earthquake cloud (including events with M < 3.5) should not be ignored.

Coulomb failure function
The Mohr-Coulomb failure criterion has been widely used to determine failure conditions on rocks (e.g., Jaeger and Cook, 1979;Hudson and Harrison, 1997).Here, we use it in a similar way to that in Reasenberg and Simpson (1992), King et al. (1994), Sumy et al. (2014), and many others in order to model static stress changes.Starting with the Mohr-Coulomb failure criterion, we can define the threshold shear stress for fault failure as shown in Eq. ( 1): where τ is the failure shear stress that is related to a particular value of the total normal stress σ and pore pressure u acting on the plane given a coefficient of friction µ and an effective cohesion c.The resultant value of (σ − u) is known as the effective normal stress and has long been known to govern the behavior of soil and rock layers in the crust (Terzaghi, 1943).Omitting the cohesive term, the Coulomb failure function (CFF), which allows us to obtain the Coulomb stress change CS as shown in Eq. ( 2), can be understood from Eq. ( 1) by moving the frictional term to the left of the equal sign.In a situation of critical stability, the result of the resultant addition would equal 0. Two important remarks are made here regarding Eq. ( 2).First, while the shear stress change is always positive and thus destabilizing, the frictional term is positive when the fault is unclamped and negative otherwise (e.g., Lin and Stein, 2004;Toda et al., 2005).Second, as we are focusing on the stress change, the cohesive term of the Mohr-Coulomb failure criterion is omitted and the equation is formulated using incremental terms .
The fault zones undergoing analysis in this paper are located below the seabed and near the injection wells of Castor UGS.Therefore, the effect of pore pressure cannot be neglected and should be accounted for in the modeling of CS (e.g., Sumy et al., 2014).After a stress redistribution, the instantaneous increase in pore pressure can be written as in Eq. (3) for the undrained condition (no fluid flow).The contribution to pore pressure is made solely by the compressional terms of the stress tensor; assuming fault zone lithologies that are more ductile than those in the neighboring region, σ xx = σ yy = σ zz (Rice and Cleary, 1976;Rice, 1992;Simpson and Reasenberg, 1994;Sumy et al., 2014) and thus where B k is the Skempton's coefficient.Equation (3) allows us to rewrite Eq. (2) as in Eq. (4) after rearranging terms: Now, by introducing the effective coefficient of friction µ = µ (1 − B k ), Eq. ( 5) is obtained: (5) Hence, it is the effective coefficient of friction (µ') that incorporates the fluid pressure effect.In addition to the fact that the hydrogeologic parameters of the faults have not been published at the time of this study and are therefore not known to the authors, knowing the exact value of µ' on a particular fault is a matter of ongoing research and was not attempted here.However, it has long been acknowledged to range from 0.0 to 0.8 (e.g., Parsons et al., 1999;Sumy et al., 2014).Its recommended value in calculations for strike-slip or faults with unknown orientation is 0.4 (friction µ = 0.75; Skempton's coefficient β κ = 0.47), as it minimizes the maximum associated error if its previously defined range (0.0 to 0.8) is taken into account (Stein et al., 1992;Toda et al., 2011a).The fault model is built into the Coulomb 3.3 software (Toda et al., 2011b), which presupposes a homogeneous elastic halfspace and implements Okada's (1992) solutions to compute strains.Bearing in mind that the entire analysis is performed at depths of 11 km or shallower (compare FM depths in Table A1), the assumption of a homogeneous crust density model is sufficiently realistic in the area of study (Díaz and Gallart, 2009).

Shortening of the seismic cycle
The characteristic earthquake theory accepts that faults will slip with a series of identical (characteristic) earthquakes (Scholz, 2002).The concept of a characteristic earthquake was first defined by Schwartz et al. (1981), who proposed that faults have an essentially limited range of magnitudes close to the maximum.According to Schwartz and Coppersmith (1984), it follows that earthquake occurrences on individual faults and fault segments do not adhere to Gutenberg and Richter's (1954) log linear frequency-magnitude relationship (log N = a − bM).This strong hypothesis was first bolstered by paleoseismicity data, although the implicit idea of equal-magnitude occurrences can be found in earlier studies (Allen, 1968;Wallace, 1970).Since that time, several publications based on seismological and geological observations have reinforced the idea (e.g., Utsu, 1972;Bath, 1981;Sieh and Jans, 1984).
Not surprisingly, such a hypothesis was a common underlying assumption in earlier studies addressing the effect of CS on earthquake recurrence times, for instance, in studies focused on the influence of moderate-to-strong events occurring on traced faults (Working Group on California Earthquake Probabilities, Dieterich et al., 1990;Parsons, 2005) and those analyzing the impact of anthropogenic seismicity on known geological structures (Baisch et al., 2009).The selection of a particular model to describe the mechanical behavior of faults is of factual meaning and affects the calculations carried out.In this paper, fault mechanics consistent with a characteristic earthquake as introduced previously are assumed.In a more global context, the extensive discussions in Wesnousky et al. (1983) and Schwartz and Copper-smith (1984) regarding which model best describes fault behavior, support this choice (see also earlier references cited in this section).Previous studies in the area of interest, particularly those concerning earthquake recurrence in the Main Fault (Fig. 1c), also used the characteristic earthquake model (Perea, 2006;Fernandez et al., 2014).
Therefore, we can define the process of strain accumulation and release on a fault as a cyclical process that takes place with a recurrence time T r , namely the seismic cycle.The introduction of an external perturbation, caused for example by fluid injection or previous earthquakes, may accelerate the process (Fig. 4).Therefore, the shortening of the seismic cycle δ cyc can be determined as in Eq. ( 6) for an associated earthquake stress drop σ and T r .However, results can be given in terms of the relative shortening (e.g., as a percentage) without introducing T r and increasing uncertainty.Both alternatives are considered here.The variable δ cyc can be incorporated to improve probability estimates of earthquake occurrences (see Parsons, 2005, and references therein).In this paper, however, we limit our contribution to the quantification of its value and assess its significance in relation to the computed stress drops and recurrence times.
In Eq. ( 6), the ratio σ / T r represents the regional stressing rate, so the average σ is released by the characteristic earthquake on its own.It involves the fact that such earthquakes produced by a particular fault are always identical in terms of the associated rupture length and dislocation.Following Lay and Wallace (1995), the earthquake stress drop can be calculated as in Eq. ( 7) by relating it to the average strain change ( u L ) using Hooke's law.In Eq. ( 7), L refers to a characteristic rupture dimension (either length or width).It is worth mentioning that the stress drop depends on a scale usually difficult to constrain precisely [L 3 ]; in our case, the shortening of the seismic cycle is only computed for the Main Fault, for which we model both rupture area and L from available geological and geophysical studies.Then, the characteristic magnitude (M w ∼ 6.0) is determined from its dimensions in accordance with Wells and Coppersmith (1994).A detailed description of fault geometry is provided in Sect. 4. .Representation of the shortening of the seismic cycle due to a perturbation that could be man-made.In this case the perturbation is of a positive nature but could be of a negative nature as well.The latter is clearly not concerning.Modified from Baisch et al. (2009).
In the context of elasticity, G is the shear modulus or Lamé's second parameter, and λ is Lamé's first parameter.Under usual crust conditions, both are roughly on the order of 3E+10 Pa.Here we compute both strike-slip and dip-slip stress drops based on the assumed fault dimensions, but maximum and minimum probable values (30 and 10 bar, respectively) are taken into account as well (Beeler et al., 2000;Baisch et al., 2009).Following Eq. ( 6), if larger stress drops were considered (e.g., Kanamori and Anderson, 1975), the influence of the computed CS value would be smaller (the ratio CS / σ diminishes).A greater value of the ratio CS / σ means increased shortening, i.e., the eventual earthquake occurrence being brought closer in time.Thus, using larger stress drops, which has the opposite effect, would be a non-conservative hypothesis from a seismic hazard standpoint.The T r value can be computed as the ratio between the expected seismic moment M e o and the geologically assessed moment rate M g o (Wesnousky, 1986;Perea, 2006).M e o can be obtained from the moment magnitude M w by using Hanks and Kanamori's (1979) 8) to (10).In Eq. ( 9), A is the area of the fault and SR is the slip rate.SR values are only reported in the literature for the Main Fault in Fig. 1c (Perea, 2006;García-Mayordomo et al., 2015), which is the only structure for which we assess δ cyc .
4 Fault model assumptions In the study presented here, cumulative Coulomb stress changes are quantified on both the FM and the mapped faults.
The model and calculation implementation is performed as follows.First, the nodal plane (NP) for each of the eight FM events is selected (see Sect. 4.1 below).Then, each FM (Sect.4.2) and mapped fault (Sect.4.3) is built in according to its geometry.Thus, the total number of modeled faults in Coulomb 3.3 is 15, resulting from 8 FM faults and 7 mapped faults (Fig. 5).The computation procedure is divided into eight steps by order of occurrence, each corresponding to slips in FM faults one after the other.After each earthquake occurrence (simulated by introducing the appropriate slip value onto the corresponding FM fault), CS values are obtained for each of the 15 modeled faults.From the second to the last event, stress changes are added to the previous ones, and thus our analysis produces cumulative CS values.Sections 4.1 to 4.3 below describe the NP selection criteria and modeling of FMs and mapped faults.

Nodal plane selection
Various criteria can be employed to select the causative NP from the two provided by each FM solution.Direct options involve comparing geological (known faults in the area) and seismological (cloud of hypocentral locations nearby) observations to both NPs to decide which one is feasible and which one is not.In our study, a straight relation with mapped faults was dismissed due to depth differences (Fig. 5); cloud density, together with location uncertainty of small magnitude events, made it impractical to select the causative faults based on adjacent seismicity (see hypocentral location discussions in Cesca et al., 2014;ICGC, 2015;Gaite et al., 2016).We use the critical pore pressure (CPP) criterion (e.g., Pine and Batchelor, 1984) in order to select each NP on the grounds that the NPs with lower CPP are better oriented for fault slip in relation to the background stress.Using the Mohr-Coulomb failure criterion as described in Eq. ( 1), we determine the necessary excess pore pressure for shear slip on each NP pair.The selected NP is that with a lower CPP value (i.e., closer to shear failure based on background stress).To determine CPP on a plane, the necessary parameters are the orientation of the principal stresses S 1 , S 2 , S 3 , their magnitude (for which the rock specific weight γ R , depth Z, and faulting regime are needed), µ, and hydrostatic pressure, P H ). We perform a statistical selection based on results from 500 realizations with likely ranges of parameters and find that the NP2 family of FM solutions (see Table 1) is consistently better oriented for failure and is thus the chosen NP for all FM solutions.An approximate and simpler approach by directly applying Coulomb 3.3 to determine the optimally oriented fault planes (OOFPs) based on a strikeslip stress regime and vertical faults also yielded the NP2 family as closer to OOFPs.A detailed explanation concerning the determination of regional stress magnitudes and the CPP approach is provided in Appendix B, while a summary of the parameters is provided in Appendix C, Table C1.
L. Saló et al.: Earthquake static stress transfer

FM-derived faults
After the NP is selected for each FM solution, geometry is modeled according to that NP, and fault plane size is determined using Wells and Coppersmith's relations (1994).
Faults are centered at the location of each earthquake.Due to the small area A of the planes (roughly 1 × 0.5 km for an M w 4 rupture), we model FM faults using only one patch.As CS values are resolved on each of the eight causative planes after each event, these planes act as source faults when they slip but also as receiver faults along the sequence.Eq. ( 11) is used to select the appropriate slip value d for each event (e.g., Lay and Wallace, 1995).From the M w , the well-known Hanks and Kanamori (1979) formula is employed to find the seismic moment (M o ); then, d can be found via isolation: (11)

Mapped faults
Taking into account the reasons given above, mapped faults near the Castor UGS only receive stress variations due to the FM faults' slip.In order to generate the input file for Coulomb 3.Then, the listric morphology is simulated in Coulomb 3.3 using a three-patch assembly along the vertical dimension of the fault, each one with a different dip angle and according to descriptions in the literature mentioned in the paragraph above.Divisions along the base of the fault depend on its length (2 to 8 patches), resulting in mapped faults being divided into a total number of patches ranging from 6 to 24.Fault traces, which if looked at in detail appear curved, were digitized into a straight line maintaining the average strike direction.Rake of the mapped faults has been inferred from FM information (as a potential rupture should be consistent with the background stress regime), the decided value being the average of the rakes from the eight FM solutions.Of the two NPs provided by each solution, the nodal plane selected to find the average is the one closest to each fault strike (the NP1 of 2 October M w 4.0 event is excluded to compute this average as it is a particular case; refer to Table 1).However, regarding the Main Fault, both were tested, and the most adverse (i.e., higher values of resulting CS) was used to provide a conservative approach from a seismic hazard standpoint.The described model for calculations is shown in Fig. 5 and Table 2.

Involved assumptions and uncertainties
Both the underlying assumptions and computation procedure of the analysis are fundamental when it comes to interpretation of the results and thus require further comments which are provided as follows.First, we used the CPP criterion for NP selection.Its application is limited to cases in which the main introduced perturbation is the fluid overpressure and stress information is not derived from FMs (i.e., when no inversion of FM solutions to infer stress magnitudes is performed).For this case, a scenario in which the pore pressure generation by injection was the main perturbation is the most likely, though no information regarding fluid overpressures has been made public (other than a supposed maximum of 7-8 bar briefly shown in Escal, 2014).In addition, what can be stated with certainty is that maximum fluid overpressures would have been considerably smaller than the tensile strength of rock (e.g., Escal, 2014;Cesca et al., 2014), and thus failure must have occurred on preexisting faults.Provided that were the case, it follows that even the main introduced perturbation would be considerably smaller in magnitude than the dominant stress state in the area; hence it can be assumed that faults better oriented for shear failure in relation to background stress were the causative ones.Lower values of CPP are obtained for faults that are better oriented regarding the regional stress, so the selection of NPs based on this line of reasoning should be solid.
As required by the approach used, we determined the regional stress magnitudes, in this case by using frictional equations on a critically stressed crust (see Appendix B).Jaeger and Cook (1979) already showed that S 1 and S 3 can be explicitly related for a fault at its frictional limit, and Chapter 4 in Zoback (2007) provides an extensive discussion on why the hypothesis of a critically stressed crust is also reasonable in intraplate regions.Briefly, the evidence is the following: (i) induced seismicity by reservoir impoundment or energy technologies involving fluid injection (e.g., Healy et al., 1968), (ii) earthquakes triggered by small stress transfer due to other earthquakes (e.g., Stein et al., 1992), and (iii) well measurements of stress (e.g., Townend and Zoback, 2000), which provide proof that stress magnitudes as obtained from the Coulomb failure theory with frictional coefficients, µ, between 0.6 and 1.0 are coherent with the measured values.
It cannot be denied, however, that uncertainty exists in the stress regime, stress orientation, and friction coefficient in addition to stress magnitudes, which are all required parameters for NP selection based on CPP.Therefore, a statistical approach as explained in Appendix B was used (see also Table C1 in Appendix C) with a likely range of variations derived from the cited references in this section and obtained FM solutions in this and previous studies.Obtained results allow the selection of each NP with confidence and, with this result, we adopt a deterministic approach (i.e., NPs are not varied) for Coulomb stress transfer computations.The trustworthiness of the obtained FM solutions via FMNEAR is From the results in Díaz and Gallart (2009), it can be seen that rock density in the first 10 km of crust (roughly where all calculations are performed) does not experience significant changes in the area of interest.Thus, the supposition of a homogeneous isotropic crust, which is intrinsic in Coulomb 3.3 calculations, can also be justified.A reliable value of Young's modulus E and Poisson coefficient ν has been established for studies of this kind (e.g., Toda et al., 2011b).Here we use 8E+10 Pa and 0.25, respectively.
Because the geometry (strike, dip, and rake) of faults, the effective friction coefficient µ', and resolved depths play a key role in the performed Coulomb analysis, we provide a sensitivity test of the previously mentioned parameters (see results in Sect.6 and Tables C2 and C3 for our best estimate (BE) and parameter variations, respectively).Other possible sources of variation from our BE result are horizontal error in the FM solution location, the width / length ratio of both FM and mapped faults, and slip level.In this study we do not consider these sources as they are either better constrained or have a lower range of deviation.The change in CS levels after parameter variation is studied for the mapped faults, mainly due to their strike and dip variation (they are modeled as straight lines, but as shown in the map in Figs.1c and 5 this is a simplification; they are also listric, so dip angle changes with depth) and rake uncertainty (unknown).On the question of stress transfer, specifically between the FM faults, an evaluation of outcome variation if different NPs had been picked is briefly discussed in Sect.7. As indicated before, after the trusted NP selection from each pair, the performed analysis is deterministic in essence, but we aim to provide insights into how variations in crucial parameters that have different levels of uncertainty would affect our results.
Finally, the assumption of a characteristic earthquake was presented in Sect.3.2.Much evidence has been provided in the literature to support the fact that the mechanical behavior on individual faults and fault segments is best described by a characteristic earthquake model.Examples from different tectonic environments derived from geological and seismological studies seem to coincide in showing that rupture occurs in specific segments of fault zones and that ruptures occur cyclically (see Wesnousky et al., 1986, and references therein).Stress drops were calculated as shown in Eq. ( 7) for both strike-slip and dip-slip regimes, and values were obtained inside the recommended range of 10-30 bar that we consider.Studies of laboratory samples and small mininginduced earthquakes have reported values as low as 2 bar, although most of them are in the aforementioned range (Mc-Garr, 1994).In contrast, those are low-end values for shallow natural earthquakes (Aki, 1972;Kanamori and Anderson, 1985), which relate better to the scenario for which we consider δ cyc .On this basis, we provide a conservative approach from a seismic hazard point of view.The calculation of T r as in Eq. ( 10), which is consistent with the elastic rebound concept (e.g., Imamura, 1930;Kanamori, 1973), requires estimating the geologically assessed moment rate M g o , which is a function of the slip rate (SR).Values that could range from a lower limit of 0.04 to 0.63 mm yr −1 are reported in the literature for the Main Fault (Perea et al., 2006;Fernández et al., 2014), and these values are a highly significant source of uncertainty in the calculation of T r .In the results, the worst assumption influence in δ cyc resulting from parameter variations in the sensitivity test is also taken into account (i.e., increased CS).A summary of the specific parameters needed for δ cyc computation is provided in Table C4.
6 Analysis results

Static stress changes resolved onto the modeled faults
Initially, the analysis checks the evolution of CS on the FM faults derived from selected NPs after each event introduces stress changes in the area, which are added to those caused by the previous events.Thus, the computed values represent cumulative stress changes.The evolution of CS in each FM fault (as modeled from the selected NP) leading to their stress state just before slipping is shown in Table 3.
Values after slipping are not shown as the relative Coulomb unloading due to slip is much larger in magnitude than any positive stress change resulting from failure in neighboring faults (they are referred to as "unloaded" in the table).Appendix D provides in Fig. D1 the CS levels immediately before failure as well.It should be noted that the first event introduces changes that are taken into account for the rest of the modeled sequence, but the CS state on its corresponding plane before failure cannot be assessed in our analysis.It can be seen that CS values on all NPs except those of the third and fourth increased as the sequence advanced and that over 0.1 bar of positive CS values are found for both the sixth and eighth events.More generally, just before slip, five of seven events (71 %; events 2 to 8 considered) had positive values of CS as a result of the previous events.The result indicates that most of the earthquakes occurred due to slip on faults that would have been brought closer to failure ( CS > 0) based solely on static stress transfer.
Cumulative static stress changes resulting from slip on FM faults one to eight were also resolved onto the mapped fault planes.The final stress state is reported in Fig. 6.For further insight, Table D1 and Fig. D2 provide a detailed time series analysis on each fault and fault patch.Stress levels increase on nearly all of the mapped faults as the sequence evolves, though the magnitude of values is modest.The East 4 fault experiences a general decrease after the occurrence of the last event (FM 8).As shown in Fig. 6, at the end of the sequence, that fault is globally inhibited for failure in terms of CS.On the other hand, resolved CS values on two patches of the Main Fault achieve 0.1 bar of positive variation, and a remarkable number of patches are left with positive values at the end of the sequence.Nevertheless, the mean CS value is close to zero.The main variations on all faults are introduced by the fifth (the largest) and eighth (which has the shallowest FM solution) events.All calculations in this part were also carried out for a standard value of µ' = 0.4.The stability of this result is addressed next.
Owing to the fact that fault geometry, µ', and depth are dominant parameters in the model, we perform a sensitivity test to explore how robust our results on the mapped faults are.Results of the analysis are shown in Fig. 7 for the East 4 and Main Fault only, where more change was registered (variations introduced are ±15 • in strike, dip, and rake, 0.2 to 0.6 for µ', and up to +3 km in depth; see Table C3 in Appendix C).The figure depicts the mean CS (as obtained from averaging all patches) resolved on both faults after the eight main events.Excluding depth, average values practically remain unchanged on the Main Fault, but average negative values lower than −0.1 bar were attained on the East 4.The nature of CS did not change.Although not shown in Fig. 7, minima of around −0.3 bar were obtained on individual patches of the East 4, while maximum values near 0.2 bar could be obtained on the Main Fault.Switching source depths to shallower points made resolved CS of greater magnitude at localized patches, which influences the average result.A switch of 2 or 3 km was seen to change the nature of CS on most faults.
In this case, we plotted the results on the mapped faults for strike, dip, rake, and depth changes on the FM faults (the ones that slip).A strike direction change on the mapped faults Table 4. Shortening of the seismic cycle on the Main Fault.The result is given in % for two values of σ ( σ ss , σ ds ) derived from the calculation as shown in Eqs.(6, 7).See Table C2 in Appendix C for the best estimate (BE) of parameters.For all variations in strike, dip, rake, µ', and depth, the case in which the computed CS was higher is selected to provide a conservative approach.See Table C3 for considered ranges and C4 for σ computation.themselves is already covered by the variation in strikes in the FM faults (the relative variation in maintaining the source fault's strike as unchanged or moving the strike in any mapped fault is analogous).The same occurs for depth if the two systems are brought closer.Dips and rakes of the mapped faults were also varied within the ranges shown in Table C3 with no greater differences than those previously presented.From this analysis we conclude that the model is robust to moderate variations in geometrical parameters of faults (those that may result from modeling imprecision when digitizing a fault, for example) and µ'.For the considered variations, the model was found to be slightly more sensitive to dip and rake variations, while µ' was of less importance.On the other hand, depth may change the results substantially regarding the magnitude of localized CS.This is reflected in the result shown in Fig. 7.

Relating the studied sequence to the occurrence of future earthquakes
Table 4 shows the shortening of the seismic cycle δ cyc computed as shown in Eq. ( 6) regarding the Main Fault.The considered magnitudes of the earthquakes with seismic cycles that may be influenced by the studied sequence are 6.0, 5.0, and 4.5.The first corresponds to an upper-bound value of the characteristic magnitude (Schwartz et al., 1981;Schwartz and Coppersmith, 1984) taking into account the total modeled area.The last two correspond to moderate shakes that could be expected in the area with relatively common frequency (i.e., every few hundred years).The mean value of the www.solid-earth.net/8/857/2017/Solid Earth, 8, 857-882, 2017  C3.Except for panel (d), which directly reports the µ' value in the x axis, only the introduced variation (e.g., −15 to +15 • in the strike direction) is reported.Note that vertical scale changes for the depth variation (e).
CS as obtained from averaging values found on all patches is used when taking into consideration the M w 6.0 earthquake (for which the entire fault area is expected to slip; assumed rupture dimensions following Wells and Coppersmith, 1994), while the maximum positive value is selected for both the M w 5.0 and 4.5 events.
In Table 4, results are provided for the best estimate (BE) of input parameters (geometry of FM faults according to selected NPs from each pair and effective friction coefficient µ' = 0.4 owing to the explanation in Sect.3.1).Results from the previously described sensitivity test are also incorporated.We report the maximum shortening that could be attained for a given value of each parameter within the considered ranges.To provide the results in percentage terms of shortening, we assume two different values of the stress drop σ within the 10-30 bar range according to the reasons given in Sects.3.2 and 5.Because stress drops for dip-slip faults are greater, the shortening when σ ds is used is lower.The computation of T r for the Main Fault yielded values with high uncertainty that discouraged their usage, mainly because of the great range of possible slip rates reported by the literature (obtained standard deviation was greater than the T r mean value).For this reason, we only report relative shortening (in percentage terms) in this paper.The specific set of parameters for computing δ cyc following Eq.( 6) is reported in Table C4.
Regarding the characteristic earthquake, δ cyc for a BE of parameters does not reach 0.05 %, and thus the advancement in the occurrence can be neglected.Shortening values around 1 % could be attained for moderate events of lower magnitudes.On the other hand, by using the values obtained in the sensitivity test for variations in geometry, depth, and µ', δ cyc never exceeds 0.2 % for the characteristic earthquake of M w 6.0.Shortening could reach 5 % regarding the M w 5.0 and 4.5 events.Such small shortening results from minor CS when compared with the associated stress drops, as obtained from the application of Eq. ( 7) for both dip-slip and strikeslip faulting.

Discussion of the results
Our results primarily indicate that (i) static stress transfer would have acted as a destabilizing trigger in the phase of the seismic sequence under analysis (from 24 September to 4 October; see Fig. 2), (ii) most of the mapped faults around the Castor UGS did not experience significant variation in their stress levels because of the tremors, and (iii) the experienced events would not have significantly shortened the occurrence of a characteristic earthquake on the Main Fault.To address points (i) to (iii) made above, we have performed an analysis solely based on coseismic static stress transfer caused by the eight strongest events of the sequence, for which FM solutions could be computed.Static stress changes have recurrently been referred to as a significant destabilizing mechanism for earthquake occurrence (e.g., Reasenberg and Simpson, 1992;Stein et al., 1992;King et al., 1994;Zoback, 2007) and have also been reported in the literature to have the potential to trigger events in contexts of fluid injection (e.g., Baisch et al., 2009;Schoenball et al., 2012;Catalli et al., 2013).Catalli et al. (2016) suggest that the effect of earthquake interactions is especially important after the end of injections.With the present study, we intend to provide evidence as to why the main phase of the sequence could have responded to static stress triggering.

Static stress transfer as a triggering mechanism
We used a deterministic approach in which the slipping NPs were chosen based on a CPP criterion that selects the bestoriented NPs for shear slip under the assumption of a fluid pressure increase.A Mohr-Coulomb envelope was used for the determination of failure conditions.It is worth noting that static stress redistribution, the investigated triggering mechanism in this study, can also promote shear failure and thus may have had a role in the slipping NP as well.For this reason and to obtain a range of variation in our results, Table 5 shows how cumulative values of CS would have changed if the slipping NP had been selected based on an assumption of maximum (MAX) cumulative CS (i.e., for each pair, cumulative CS values are compared and the one with a higher value is picked as host of the next event).We obtain similar values, although for the MAX assumption, six instead of five NPs could have had positive values of CS immediately before failure.When selecting MAX NPs by trial and error (trying both NPs each time), it was observed that, except for one case, choosing the other NP did not change the resolved CS nature (positive or negative) but only its magnitude, which is something usually found in calculations (e.g., Sumy et al., 2014).In summary, our result for cumulative CS on NPs yields positive values of CS for most of the events (five to six out of seven), although only two of them (the sixth and eighth) reach values of 0.1 bar (0.01 MPa).Many papers have addressed the issue of static stress triggering by earlier earthquakes.In studies considering triggered natural events, a threshold of 0.1 bar has commonly been found (e.g., Reasenberg and Simpson, 1992;King et al, 1994;Harris, 2000), although it should be noted that this threshold value can span more than an order of magnitude.Lower threshold values have been reported for anthropogenic seismicity.Orlecka-Sikora (2010) found changes as small as 0.05 bar to be statistically significant in the triggering of events at the Rudna mine in Poland; in another mine-induced seismicity study, Kozlowska et al. (2015) also suggest that earthquakes were promoted in areas where CS increased not more than 0.1 bar.With an increase in detection capability and the analysis of smaller earthquakes (lower associated stress drops), it seems that the threshold values diminish.This could be an indicator that triggering is not a threshold process, an idea that is not new.Hardebeck et al. (1998) point out that no theoretical reason to explain the appearance of a threshold value exists.They specifically note that "an arbitrarily small static stress increase should be able to trigger an earthquake on a plane arbitrarily close to failure".Ziv and Rubin (2000) studied a dataset of 63 M ≥ 4.5 events in California and concluded that no lower threshold exists, an idea that is also supported by Ogata (2005).
Perhaps the best option to determine whether CS values can explain the triggering of a particular sequence is by comparing them with a randomized control set of events to see whether statistically significant differences are found (e.g., Ziv and Rubin, 2000;Orlecka-Sickora, 2010).Provided the Coulomb stress transfer hypothesis is valid, one would expect a remarkably higher number of events to be found in areas with positive CS for the real set when compared to the control group.A probabilistic approach in which NPs are varied can also help in determining the maximum range of variation and nature of CS when NP ambiguity cannot be disclosed.However, we followed a deterministic approach as the employed method of NP selection showed a clear preference for NP2 family (refer to Table 1).The CS values were seen to build up along the considered set of events (larger cumulative values for the events occurring at the end), which is something that responds to the geometry of slipping planes and their relative locations.An NP selection based on MAX CS, as discussed earlier and shown in Table 5, does not seem to alter the general picture of our results, which show a significant majority of the occurred events taking place on NPs in which failure is promoted by CS (71 %).Although for such a small dataset, this percentage value might not be very significant, it is in the range obtained by previous studies taking into account a larger number of events that have deemed CS to be relevant (e.g., Ziv and Rubin, 2000;Orlecka-Sikora et al., 2009;Orlecka-Sikora, 2010).
www.solid-earth.net/8/857/2017/Solid Earth, 8, 857-882, 2017 Because it was found that selected NPs were close to OOFPs and could also be very close to failure conditions, the hosting faults could have been fairly well oriented in relation to the regional stress field (see the distance of selected NPs to failure envelope in Fig. B1; P values as low as 8 bar were obtained for a best estimate of the parameters but could be even less if the entire range of possible variations as provided in Appendix B is used).As commented on in Sect.5, in a critically stressed crust (e.g., Stein and Lisowski, 1983;Bak and Tang, 1989) small-magnitude stress changes may be enough to initiate a rupture, and thus, following the comments provided, it seems coherent to think that static stress transfer played a destabilizing role in the considered events.
In fact, the development of fault rupture could be so subtle that even a minimum triggering stress at a local point may initiate a failure that, once started, could propagate throughout the total fault area (even to regions where the external perturbation does not surpass the fault's friction-based shear resistance) and generate felt events if dimensions allow it (Gischig, 2015;Piris et al., 2017).

The hosting faults
Due to depth differences, our FM solutions could not be related to the mapped faults, but even if they were pictured at the mapped fault depths, location uncertainty and fault complexity near the reservoir make it impractical to directly ascribe them to an existing structure.On the basis of betteroriented NPs for shear failure, we selected the NE-SW family of NPs (NP2 in Table 1), which contains solutions of high similarity.This approach reinforces a scenario in which various SE-dipping fractures with similar geometry (or a main unique structure) underlying the reservoir would have hosted the strongest phase of the sequence.For our next point, we assess how this result relates to information provided in the literature for the seismic sequence under analysis.Cesca et al. ( 2014) located a dataset of 73 events and obtained depths ranging from 0 to 4 km and a general orientation of the earthquake cloud NE-SW to NNE-SSW.Conversely, Gaite et al. ( 2016) located a dataset of 161 events and obtained most of the hypocenters at around 6 km in depth and a seismicity alignment NW-SE.Previous studies resolving FMs by Frontera et al. (2013) and IGN (2013) obtained solutions that are very similar to ours, while Cesca et al. (2014) obtained planes of similar strike but resolved at shallower depths and with gentler dips regarding the NP2.All of these FM solutions can be found in Table A1.Cesca et al. (2014) support the rupture of a NNE-SSW striking plane located within the reservoir (East 4) on the basis of better orientation to background stress and also discard the Montsia system (NW-SE striking) due to dip difference with that of the compatible NP family in terms of strike direction.Geometrically speaking, our NP selection could agree with the rupture model that their results mostly support, in which a nearly crit-ically stressed SE-dipping structure would have hosted the seismicity.
The analysis presented here also quantified cumulative CS values on the mapped faults around the Castor UGS resulting exclusively from slip on the FM faults (modeled according to selected NPs).We aimed to examine how much the eight main events affected those faults and what their final stress state was in terms of CS.Time analysis on all modeled patches of the mapped faults (Table D1 and Figure  D2) showed two main facts.First, though always small in magnitude, the fifth and eighth events have the highest influence in terms of CS; and second, at the end, a localized area of the Main Fault is loaded in terms of CS, while most of the East 4 fault is unloaded (Fig. 6).The most relevant CS discharge on the East 4 is obtained after the last event, represented by FM 8 (resolved at 3 km under the seabed and with approximately compatible geometry).We cannot formally exclude the fact that if a hidden structure exists at the exact location where that FM fault was placed, then the East 4 fault could be further from rupturing; but, given the proximity and similar characteristics, this result opens the door to the argument that, should one of the mapped faults have slipped, then it was most likely to have been the East 4 (i.e., at least the FM 8 solution could be the result of slip on the East 4).As reported in Sect.6.1, this result was robust to moderate variations in fault geometry and µ'.Conversely, an important depth error in FM solutions (for example, of more than 2 km) could substantially change the result on the East 4 (Fig. 7) and therefore be the basis for a different discussion.However, what remains unchanged is the link between the FM faults and the East 4, the latter always being the most sensitive to CS as a result of slip on the selected NPs.
Therefore, because of (i) similar geometry between the selected NPs and the East 4 and (ii) resolved CS values on that plane, we consider that structure to be the most likely of the mapped faults to have hosted the shallowest events in the sequence should one of them have slipped.This structure is well oriented when it comes to background stress and was found to be quite similar in strike and dip to multiple selected NPs (those of the NP2; see Table 1).Nevertheless, the latest hypocentral locations, which were mostly resolved at around 6 km of depth as indicated before, do not support its rupture (ICGC, 2015;Gaite et al., 2016).The locations by Gaite et al. (2016) were obtained with a 3-D shear wave velocity model specifically for the area.Moreover, as our FM solutions were essentially found between 5 and 8 km beneath the seabed and thus agree with the findings of Gaite et al. (2016), our study supports slip on one or various deeper structures when it comes to the strongest phase of the sequence.
As for the Main Fault, its activation during the felt events is not supported by our solutions nor the previously reported FM solutions (Table A1), none of which would be compatible with its geometry.Nevertheless, it should be stressed that the experienced seismicity could have promoted failure in this structure to some extent, as shown by its loaded patches at the end of the sequence (Fig. 6).Next, we discuss how much this CS loading could have brought the occurrence of a characteristic earthquake on the Main Fault closer.

Could the occurrence of future earthquakes have
been shortened as a result of static stress transfer?
After resolving CS on the mapped faults as generated by slip on the eight FM faults, we computed the shortening of the seismic cycle δ cyc on the Main Fault under the assumption of a characteristic earthquake.The influence of the eight strongest events on the occurrence of moderate events of M w 5.0 and 4.5 was also investigated.The results for a best estimate of the involved parameters, as presented in Sect.6.2, indicate a 0.05 % shortening for the characteristic earthquake of M w 6.0, while approximately 1 % shortening could have been attained for M w 5.0 and M w 4.5 shakes.Shifting the depths of FM solutions upwards introduced the highest deviation from our BE result.However, based on the considered ranges, shortening of the M w 6.0 event would never have increased 0.2 %, and the upper-bound limit for the M w 5.0 and 4.5 events would be 5 % (see Table 4).In addition to CS based on static stress redistribution, we worked with the hypothesis of a characteristic earthquake and calculated stress drops σ and return periods T r as described in Sects.3.2 and 5. Including T r in the results was finally discarded owing to the high uncertainty in the calculation.Taking into consideration our results summarized above, the seismic cycle of magnitude events ranging from 6.0 to 4.5 in the Main Fault could not have been brought significantly closer in time.With our configuration of selected NPs, it was seen that even when shifting FM faults 3 km shallower, the mean value of CS on the Main Fault would not increase, but rather decrease (last panel in Fig. 7).Such a result supports our statement that the shortening of the seismic cycle of a future characteristic earthquake should be minor.
The seismic cycle and earthquake occurrence are treated from a probabilistic standpoint (e.g., Parsons, 2005;Baisch et al., 2009), and hence this should be accounted for when handling return periods and shortening.Small values of relative shortening as obtained here indicate that the probability of occurrence of a particular seismic event was not remarkably increased, but this should never be used to rule out a later occurrence (e.g., McGarr, 2014;Mulargia and Bizzarri, 2014).Moreover, while we can quantify the effect of static stress transfer and evaluate the introduced shortening in the seismic cycle of the Main Fault, it should be noted that our computations have been solely based on static stress transfer.Bearing in mind that gas was being injected into the depleted Amposta reservoir when injections started, fluid overpressures are likely to have been the main introduced perturbation in the system.Because of that, considering hydromechanical coupling would probably have resulted in a higher magnitude of perturbation and therefore increased shortening, which cannot be accounted for in our analysis.

Comment on shortcomings
For the NP selection we did not take into account other mechanisms that can produce shear slip, such as poroelastic stress changes, differential compaction, mass changes, thermal effects, or chemical reactions, that may change friction conditions (e.g., Ellsworth, 2013).Our reasoning is that, for the time span and project type under analysis, (i) all other perturbations must have been less significant than fluid overpressures (probably an order of magnitude at least), and (ii) slip should have occurred along faults better oriented to regional stress for shear slip.However, the Castor UGS was placed in the depleted Amposta oil field, and thus the presence of destabilizing stress changes resulting from the exploitation of the reservoir cannot be excluded.And, most importantly, our analysis did not take into account fluid flow.While the present study has allowed for the quantifying of a specific earthquake trigger to assess its influence on the strongest phase of the sequence and in the future occurrence of earthquakes, it cannot address the issue of the origin of seismicity.As the main man-made perturbation was the injected fluid, the physical mechanism of pore pressure generation due to fluid flow should at least be included in the approach to assess this origin.This was beyond the scope of this study due to a lack of available public data and the different aims of the analysis.For this reason, the results and discussion presented here refer only to the influence of the studied triggering mechanism in the main phase of the seismic sequence, which took place after the injections had been halted, and 3 to 4 weeks after the onset of seismicity.

Conclusions
Our analysis points towards static stress transfer having been a destabilizing mechanism in the series, as shown by positive cumulative CS on most selected NPs. Background stress analysis yielded that NE-SW striking, SE-dipping planes are better oriented for shear failure than NW-SE striking, SWdipping planes.Owing to the high confidence in NP selection and the reasons given in Sect.5, the sensitivity analysis of the results was limited to the mapped faults; variations in fault geometry (strike and dip), rake, effective friction, and depth were introduced (see Sect. 6.1).A more comprehensive probabilistic uncertainty analysis would be useful but was deemed to be beyond the scope of this study.Thus, the performed static stress transfer analysis based on the combined information of background stress and FM solutions supports the following: (i) the strongest events could have been caused by the rupture of a single or various hidden NE-SW-oriented structures underlying the reservoir, and (ii) the mechanism of static stress redistribution could have triggered the strongest earthquakes on its own.
Due to an important depth difference, which is also supported by the latest study on hypocentral locations (Gaite et L. Saló et al.: Earthquake static stress transfer al., 2016), we cannot relate the events to the mapped faults around the reservoir; despite that, it was inferred that the East 4 (i.e., an SE-dipping structure within the reservoir) is the most likely to have been activated should one of the mapped faults have slipped in the series.Although Cesca et al. (2014) obtained a distribution of the earthquake cloud that could fit with certain regions of the Main Fault, its slip in the series is not supported by FM solutions in this or any of the previous studies.Consistent with our findings, the seismic cycle concerning the Main Fault characteristic earthquake (around M w 6.0) was not shortened by the experienced events.Occurrences of smaller events could have been hastened to some extent.
At the time of the present study, no analysis that incorporates fluid flow, which is most likely the main introduced perturbation, has been performed.Its consideration in a hydromechanical model should help shed light on the temporal history of the whole seismic sequence.Including fluid flow in the model would result in higher stresses near the injection wells, and its effect should definitely be incorporated when future studies assess the origin of seismicity in the series.Exhaustive knowledge regarding faults and the stress state, together with local seismic monitoring, is essential in order to diminish uncertainty in this type of analysis.
Data availability.The earthquake hypocentral locations used in this paper for FM computation of the studied events are provided by the ICGC (2015).Waveform data from ICGC stations (see Fig. 1a) are available free of charge at http://www.icgc.cat/en/Public-Administration-and-Enterprises/Services/Earthquakes/ About-the-seismic-and-accelerometric-network/Seismic-Data (see also ICGC, 2000).Data from other agencies, such as IGN and OEBR (see Fig. 1a), should be accessible by contacting the controller (Instituto Geográfico Nacional, IGN; Observatori de l'Ebre, OEBR).All remaining data required for this analysis (e.g., fault locations and geometry) are available in the text references.
Appendix A: FM solutions A summary table of the FM solutions corresponding to the eight felt events is provided (Table A1).The World Stress Map (WSM; Heidbach et al., 2008) reports both strike-slip and normal faulting stress regimes at some distance from the Castor UGS and a strike direction of the S H compression (main horizontal stress) close to 10 • (NNE).Slight preference for strike-slip regimes is, however, displayed.Schindler et al. (1998) provide similar results in their analysis with a variation in the S H from 8 to 36 • NNE at some 10 km from the area of interest in this paper.A strikeslip/normal stress regime is also reported.We follow the approach by Cesca et al. (2014) in their discussion and select a range of variation for the orientation of the principal background stresses, as explained below.To the authors' knowledge, no information regarding the regional stress magnitudes at depth has been reported in the area.However, for a given stress regime, the background stress magnitudes can be approximated through the use of frictional equations on a critically stressed crust.Further discussion on the validity of this assumption is provided in Sect. 5 in the main text.As reported by Jaeger and Cook (1979), the relation between the maximum (S 1 ) and minimum principal stress (S 3 ) magnitude for a fault at the frictional limit are given by Eq. (B1), where P H is the hydrostatic pressure and µ is the coefficient of friction: The overburden load (vertical stress, S V ) can be directly approximated by using the rock specific weight and depth of interest.According to Anderson's theory (Anderson, 1951), for a normal faulting regime S 1 equals S V , whereas for a strike-slip one S 2 = S V .Another constraint is that based on Eq. (B1), the ratio of the maximum principal effective stress (S 1 ') to the minimum principal effective stress (S 3 ') is constant and depends on µ (e.g., 3.1 for µ = 0.6).The obtained ratio is an upper bound for the relation between S 1 ' and S 3 ', which occurs when faults are critically stressed (on the verge of failure).Thus, more generally for any stress regime, In the area of interest, the dominant stress regime seems to be a combination of both normal and strike-slip, which for S h Table A1.
FM solutions summary table sorted by input magnitude (highest to lowest).Solutions obtained by Frontera et al. (2013), IGN (2013), andCesca et al. (2014) are also shown to ease comparison.The Frontera et al. ( 2013) solutions were obtained using the same method (FMNEAR) and model (default).

FM solution table
Earthquake location data (ICGC, 2015) (minimum horizontal stress) near its lower limit (S 3 min ) can be found when S H ≈ S V (S 1 ≈ S 2 ).That is the result of the inequality in Eq. (B2) that can be satisfied with both S V and S H as S 1 (see Chapter 4 in Zoback, 2007, for an insightful analysis on the topic).Owing to what is reported in the cited references in this section, a normal/strike-slip stress regime is assumed here for computation of the magnitudes of S 1 , S 2 , and S 3 , and we consider an S 2 / S 1 ratio between 0.8 and 1 and a range of S 3 from its lower limit to 1.2 of that value.As indicated before, hydrostatic pressure is considered, and the water table above the seafloor, which is less than 100 m in the Castor platform surroundings, is neglected.More comments follow.

B2 Determination of critical pore pressures
Thinking in terms of a Mohr-Coulomb failure envelope (Eq. 1 in the main text), it is widely accepted that fluid pressure decreases the effective normal stress acting on a fault plane, and if shear stress overcomes the frictional resistance, failure occurs.In the development of underground gas storage (UGS), pressurized fluid is injected into a rock formation, thus reducing the value of the right-hand term in Eq. ( 1) in the main text and bringing conditions closer to failure due to shear slip.This model is commonly used in induced seismicity analysis.In this section, we show how to determine the normal and shear stresses acting on a plane given a configuration of the background stress, and then compute the necessary pore pressures for shear slip.
For a given value of the orientation and magnitude of the principal stresses (S 1 ≥ S 2 ≥ S 3 ) and geometry of the objective fault plane (in this case, each of the two NPs provided by each FM solution), we obtain the normal and shear stress on a plane as shown in Eqs.(B3) and (B4) below (e.g., Zoback, 2007;Mukuhira et al., 2016): where, in both equations above, λ j is the directional cosine from each main stress direction to the normal vector of the fault plane.The CPP (total pore pressure necessary for shear slip) can be written as in Eq. (B5) by adding the hydrostatic pore pressure P H to the needed excess pressure P : CPP = u = P H + P .(B5) By substituting Eq. (B5) in Eq. ( 1) considering a Mohr-Coulomb envelope for shear failure with cohesion = 0, the needed P is obtained as in Eq. (B6).We compute P for the NP pairs of each FM solution and pick the NP with Each circle shows the normal and shear stress resolved on a particular NP according to number and color inside the circle (numbers are events 1 to 8, and color is NP family; 1 is blue and 2 is yellowish, as indicated in Table 1 in the main text).Note that NP2 is closer to failure, although the precise location of each NP is lost in this figure as circles have been substantially enlarged to ease identification.lower P for shear slip as the causative plane of a particular event.Clearly, lower P results in lower CPP.Those faults are closer to the failure condition for shear slip.

B3 Statistical approach
To determine CPP on an arbitrarily oriented plane, the needed parameters are the orientation of S 1 , S 2 , and S 3 , their magnitude (for which the rock specific weight γ R , depth Z, and faulting regime are needed), µ, and P H .In order to perform the computation, we assume the following (see Sect. 5): (i) S H orientation 25 • ± 25 (N to NNE), (ii) 0.8(S 1 ) ≤ S 2 ≤ S 1 , (iii) S 3 min ≤ S 3 ≤ 1.2 (S 3 min ), (iv) 0.6 ≤ µ ≤ 0.9, (v) normal or strike-slip stress regime, (vi) γ R between 20 and 26 kN m −3 , and (vii) P H as obtained at the depth of interest with a water specific weight of 10 kN m −3 .Note that for each NP pair (refer to Table 1 in the main text), the computation is performed at the depth of the FM solution, and thus principal stress magnitudes and P H do change.The selection of NPs is made via percentage statistics of 500 computations, each one choosing a random value for (i) to (vi) within the indicated ranges (Fig. B1a).All values within the specified ranges have equal probability of being selected.The results accord with our best estimate (BE) of slipping NPs (Fig. B1b), which supposes a strike-slip faulting regime (derived from both of our FM solutions, which are essentially strike-slip mechanisms, and references of the stress field reported previously), S H orientation of 23 • N (Schindler et al., 1998;Heidbach et al., 2008), S V = 0.95S H ≈ S H (as discussed in Sect.B1), and S 3 min as obtained from Eq. (B1) with µ of 0.75 (center of the assumed range); γ R = 23 kN m −3 , a low-end value for limestone rock (according to rock type in the Amposta field, where the Castor UGS was to be placed).In Fig. B1, it can be observed that the NP2 family is consistently better oriented for shear slip (i.e., lower P ), and thus is the picked one for every FM in our study.All parameter ranges for this computation are also presented in Appendix C, Table C1.
Table C1.Needed parameters for the calculation of CPP on each NP and their chosen ranges in addition to the FM solution and their depths (refer to Table A1).

Figure 1 .
Figure 1.Location of the study; the Castor platform is represented with an orange star.(a) General view showing existing permanent seismic stations around the UGS and main stations in the Catalonia area.The ECOL station was set up after the beginning of seismicity and was not used in the majority of the available earthquake locations.See legend for agencies.(b) Sketch of the UGS structure after the Shell seismic profile in Seeman et al. (1990); vertical scale is in two-way travel time.A rough approximation of the reservoir body is depicted in yellow, and stratigraphic markers are shown in blue and green (modified from Cesca et al., 2014).(c) Detailed view of known faults in the area modified from Cesca et al. (2014) according to the red area highlighted in panel (a).The map view is plotted at the approximate depth indicated by the discontinuous line in panel (b).

Figure 3 .
Figure 3. (a) Output map with the FM solution for the M w 3.9 that occurred on 30 September 2013, its location, and the stations used to compute it (green triangles).(b) Waveform fit at stations ERTA and CMAS for the named event.The recorded waveforms are plotted with a continuous line; the discontinuous line is adjusted.
The FMs obtained correspond to strike-slip mechanisms with some normal component.Excluding depth, they are similar to those in Cesca et al. (2014), though with steeper SE-dipping planes, and are in clear agreement with the Frontera et al. (2013) and IGN (2013) solutions.Detailed information on the modeling of earthquake sources from FM solutions is provided in Sect. 4.

Figure 4
Figure4.Representation of the shortening of the seismic cycle due to a perturbation that could be man-made.In this case the perturbation is of a positive nature but could be of a negative nature as well.The latter is clearly not concerning.Modified fromBaisch et al. (2009).

Figure 5 .
Figure 5. (a) Obtained solutions for the eight main events.Dots indicate the location of each beach ball plot.Mapped faults are indicated as well.In the legend, NS refers to mechanisms with normal component, whereas the others (SS) are essentially pure strikeslips.(b) 3-D view of the fault model used in Coulomb calculations with both the FM and mapped faults and their associated slip in our modeling (only the FM faults slip).View from the SE.The coordinate origin (0, 0, 0) is shown in the mesh, and it corresponds to [0 • 42 40 • 24 ] in the map above.

Figure 6 .
Figure 6.Final CS state after the eight studied events on all 15 modeled faults.Smaller rectangular patches, which are deeper, correspond to the eight FM faults.Two different views of the same figure are shown to obtain a complete view of the modeled fault system.

Figure 7 .
Figure 7. Sensitivity analysis result on the East 4 and Main Fault (MF).The figure shows average values (as obtained from all modeled patches).Squares refer to the best estimation of parameters, and lines depict the CS variation resulting from each parameter change as reported in TableC3.Except for panel (d), which directly reports the µ' value in the x axis, only the introduced variation (e.g., −15 to +15 • in the strike direction) is reported.Note that vertical scale changes for the depth variation (e).
Appendix B: Selection of NPs based on critical pore pressures (CPPs) B1 Background stress information

FM
solutions summary table.Sorted by input magnitude (highest to lowest).Solutions obtained by Frontera et al. (2013), IGN (2013) and Cesca et al. (2014) are also shown to ease comparison.Frontera's et al. (2013) solutions were obtained using the same method (solutions summary table.Sorted by input magnitude (highest to lowest).Solutions obtained by Frontera et al. (2013), IGN (2013) and Cesca et al. (2014) are also shown to ease comparison.Frontera's et al. (2013) solutions were obtained using the same method (solutions summary table.Sorted by input magnitude (highest to lowest).Solutions obtained by Frontera et al. (2013), IGN (2013) and Cesca et al. (2014) are also shown to ease comparison.Frontera's et al. (2013) solutions were obtained using the same method (solutions summary table.Sorted by input magnitude (highest to lowest).Solutions obtained by Frontera et al. (2013), IGN (2013) and Cesca et al. (2014) are also shown to ease comparison.Frontera's et al. (2013) solutions were obtained using the same method (solutions summary table.Sorted by input magnitude (highest to lowest).Solutions obtained by Frontera et al. (2013), IGN (2013) and Cesca et al. (2014) are also shown to ease comparison.Frontera's et al. (2013) solutions were obtained using the same method (

Figure
Figure B1.(a) Bar plot showing the selected NP based on CPP for each of the 500 computations.Note that NP2 is picked nearly every time for all FM solutions.(b) Normalized 3-D Mohr's circle and failure envelope for the best estimation as discussed in Sect.B3.Each circle shows the normal and shear stress resolved on a particular NP according to number and color inside the circle (numbers are events 1 to 8, and color is NP family; 1 is blue and 2 is yellowish, as indicated in Table1in the main text).Note that NP2 is closer to failure, although the precise location of each NP is lost in this figure as circles have been substantially enlarged to ease identification.

FM
Figure D1.Cumulative Coulomb stress changes resolved on the nodal planes as obtained just before slip of each one.Each subplot shows the FM solution with cumulative CS values on both the NP(s) of the old and the upcoming event.Past NPs are shown as red rectangles (faults with slip).Note that the color bar scale saturation changes as the sequence evolves.
Figure D2.Time series of CS computed onto each patch of the mapped faults.The discontinuous vertical lines indicate the occurrence of an earthquake (note that day 267 has one occurrence, and days 275 and 277 have two) with the focal mechanism indicated in the left panel (a).Red lines indicate that the final CS value on a patch is positive, while blue is used for those that are negative (inhibited) at the end.The gray weighted line on each panel represents the mean evolution as reported in TableD1.Vertical scale changes among panels.

Table 1 .
FM solutions summary showing both NP families (1 and 2) as referred to in the text and other figures.

Table 2 .
Information for the geometrical fault model used in Coulomb 3.3.Both FM and mapped faults are described.Regarding the last two columns, CS are tracked on all faults and the FM faults act as both sources (slip) and receivers, whereas the mapped faults act only as the latter.

Table 3 .
Cumulative CS evolution on each FM fault.After fault activation on a particular date (act.), it is indicated as unloaded (unl.), and no numeric value is provided.The arrows next to each value indicate whether the fault plane experiences an increase or decrease in Coulomb stress with respect to its previous state.At the beginning of the series, we suppose CS = 0. Values are in bar units.

Table 5 .
Nodal plane pairs for each FM solution.Selected NPs based on a CPP criterion and maximum cumulative CS (MAX) before occurrence are shown, as are the resolved cumulative CS values on each NP just before slip.

Table C2 .
Best estimate (BE) of parameters in order to compute CS.
r , and stress regime

Table D1 .
Cumulative CS evolution (mean values obtained by taking into account all patches of a particular fault) on each mapped fault.The arrows next to each value indicate whether the fault plane experiences an increase or decrease in Coulomb stress with respect to its previous state, while a dash indicates that there is no variation.Mean resolved values of CS on mapped faults[bar]

Table D1 .
Vertical scale changes among panels.