Simulation of seismic events induced by CO2 injection at In Salah,

Carbon capture and storage technology has the potential to reduce anthropogenic CO 2 emissions. However, the geomechanical response of the reservoir and sealing caprocks must be modelled and monitored to ensure that injected CO 2 is safely stored. To ensure conﬁdence in model results, there is a clear need to develop ways of comparing model predictions with observations from the ﬁeld. In this paper we develop an approach to simulate microseismic activity induced by injection, which allows us to compare geomechanical model predictions with observed microseismic activity. We apply this method to the In Salah CCS project, Algeria. A geomechanical reconstruction is used to simulate the locations, orientations and sizes of pre-existing fractures in the In Salah reservoir. The initial stress conditions, in combination with a history matched reservoir ﬂow model, are used to determine when and where these fractures exceed Mohr–Coulomb limits, triggering failure. The sizes and orientations of fractures, and the stress conditions thereon, are used to determine the resulting micro-earthquake focal mechanisms and magnitudes. We compare our simulated event population with observations made at In Salah, ﬁnding good agreement between model and observations in terms of event locations, rates of seismicity, and event magnitudes.


Introduction
In order to mitigate anthropogenic greenhouse gas emissions without reducing our use of fossil fuels, CO 2 can be captured at large, point source emitters (such as coal-fired power stations). This CO 2 must then be sequestered in suitable geological repositories, such as mature hydrocarbon reservoirs or deep saline aquifers. Pacala and Socolow (2004) estimate that over 3.5 billion tons of CO 2 per year must be sequestered to make a significant impact on global emissions.
Carbon capture and sequestration (CCS) is seen by many as a vital technology for meeting emissions objectives (International Energy Agency, 2012), and a number of commercial-scale CCS projects are planned in the coming decade. However, with the recent increase in seismicity induced by injection of large volumes of oilfield waste fluids (Ellsworth, 2013), concerns have been raised over the feasibility of injecting large volumes of CO 2 for CCS. In addition to the hazard posed by induced seismic events, geomechanical deformation at CO 2 storage sites, whether seismic or aseismic, Fig. 1. Schematic illustration of the In Salah gas field. In (a) we show the extent of the gas-field in the crest of the anticline. Gas is produced from 4 wells (magenta), and CO 2 is injected via 3 wells (green) to the sides of the gas reservoir. A microseismic monitoring well was placed close to CO 2 injection well KB-502. Inset in (a) shows the field location in the centre of Algeria (black square). In (b) we show a schematic cross section: the reservoir unit is 20 m thick at a depth of 1850-1900 m, overlain by 900 m of Carboniferous mudstone, which acts as a seal. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) a need to develop techniques capable of predicting seismic activity induced by geomechanical deformation.
In this paper we develop an approach that uses a geomechanical model, coupled with a history-matched reservoir model, in order to model microseismicity induced by fluid injection. We demonstrate this approach by matching our model predictions to observed seismicity at the In Salah CO 2 injection project. Fig. 1 illustrates operations at In Salah. Natural gas containing 5-10% CO 2 is produced from several gas fields. The CO 2 was stripped and re-injected into the water-leg of the Krechba gas reservoir at rates of up to 50 mmscf/day through 3 wells. In total, 3.86 million tons of CO 2 were re-injected from 2004 to 2011. The reservoir consists of a 20 m thick Carboniferous sandstone, with porosity of 13-20% and permeability averaging 10 mD . Injection depths range from 1850-1950 m. The reservoir is sealed by 950 m of Carboniferous mudstones, over which lies a sequence of fresh-water bearing Cretaceous sandstones.

CO 2 injection at In Salah
CO 2 injection has lead to substantial increases in pore pressure around the injection sites (Bissell et al., 2011), which in turn has produced geomechanical deformation. Satellite InSAR measurements of surface uplift at the injectors first indicated scale of geomechanical deformation (e.g. Onuma and Ohkawa, 2009). Subsequent inversions of surface uplift data indicated the possibility of active fracture networks in the reservoir (Vasco et al., 2010).
In 2009, 5 yr after the start of injection, a microseismic monitoring array was installed near to injection well KB-502 (Stork et al., 2015), where the most intense fracturing had been inferred from InSAR observations (Vasco et al., 2010). A string of geophones was deployed in a vertical borehole extending to a depth of 500 m. However, logistical issues with the deployment (described in Stork et al., 2015) meant that only the uppermost station, at a depth of 80 m, provided fully useable data, while a second geophone provided useable data from the vertical component (but not the horizontal components). Nevertheless, Stork et al. (2015) generated P -and S-wave picks for 6280 events between August 2009 to June 2011, using P -wave hodogram analysis in combination with the differential arrival time between P -and S-waves (t S-P ) to constrain event locations. Stork et al. (2015) also computed event magnitudes, finding that the largest event had a magnitude of M W = 1.7, and that the events followed a Gutenberg-Richter distribution with b = 2.17 ± 0.09. Elevated b values are commonly observed for seismicity induced by injection (Verdon, 2013). Goertz-Allmann et al. (2014) also analysed microseismic data from In Salah, using a master-event cross-correlation technique to detect events, finding a total of 3651 events with P -and S-wave picks. We attribute the substantial difference in the number of events detected to differences in the event detection methods used by Stork et al. (2015) and Goertz-Allmann et al. (2014) respectively. Nevertheless, both studies find similar changes to the rate of seismicity through time, a similar range of event magnitudes, and similar Gutenberg-Richter b-values.

Method -modelling induced seismic events using geomechanical models
It is commonly assumed that seismic activity in hydrocarbon reservoirs occurs on pre-existing faults and/or fractures (e.g. Segall, 1989). Pore pressure changes directly change the effective stress in the reservoir, and compaction or inflation of the reservoir leads to stress changes in the rocks surrounding the reservoir. Changes in the effective stress field can be resolved into changes in effective shear and normal stresses acting on potential slip planes. Where these stresses exceed the Mohr-Coulomb failure criteria, slip can occur and a seismic event is triggered.
On this understanding, if the location, size, orientation and frictional properties of faults and fractures in a reservoir are known, or can be approximated; and the initial stress conditions, and the changes in stress caused by injection or production can be modelled, then it should be possible to simulate when and where fractures exceed the failure envelope and trigger a seismic event. This understanding forms the basis of the workflow developed here.

Modelling fractures at In Salah
Our approach requires that the locations, sizes and orientations of fractures in the reservoir are known or can be modelled. For our application to In Salah, we use the results described by Bond et al. (2013). Where such information is not available, statistical fracture distributions can be generated (e.g. Goertz-Allmann and Wiemer, 2013). Bond et al. (2013) use a geomechanical reconstruction method to model fractures at In Salah. A structural model of the faulted, folded reservoir was generated from seismic reflection data. This deformed model was restored to an undeformed template using a geomechanical algorithm (Shackleton et al., 2009). This algorithm is based on a mass-spring solver (Provot, 1995), where each mesh element acts as a spring with specified stiffness, and the joining points as masses. A final, 'restored' state is found by iteratively minimising the energy in the springs.  Bond et al., 2013). In (a) we map the spatial variation in fracture intensity (contours show log 10 (N) per 50 × 50 m bin). Injection wells are marked in green, production wells in magenta. In (b) we show histograms of the strikes of fault (blue) and fold (red) associated fractures, and in (c) we show dip histograms. In (d) we show the cumulative length distribution for fault (blue) and fold (red)-associated fractures. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Bond et al. (2013) then forward modelled the restored surface, using the same algorithm, to predict the spatial distribution of strain intensity and orientation through time. The strains generated by the forward model were used to produce the fracture model. Two fracture sets were generated, one associated with folding and the other with fracturing. In total, over 300 000 fractures were modelled. Fig. 2 shows the spatial variation in fracture intensity, and fracture strikes, dips and lengths. Fault-associated fractures were modelled with a maximum length of 150 m. This corresponds to the maximum fault offsets. Fold-associated fractures were modelled with a maximum length of 200 m. Given a reservoir thickness of 20 m overlain by mechanically heterogeneous strata (inter-bedded sands and shales), and that mechanical interfaces inhibit frac-ture propagation longer fractures are thought unlikely. For both sets, the minimum modelled length was 50 m. These bandwidths (50-200 m) also correspond to prior work by Iding and Ringrose (2010), who used discrete fracture networks (DFN) to model injectivity at In Salah. Their DFN was based on seismic data, core samples and image logs, and was able to reproduce observed well tests results from In Salah.

Modelling of reservoir fluid flow
We use a history matched reservoir model to simulate pore pressure changes caused by gas production and CO 2 injection (Bissell et al., 2011). The simulation covered an area of 1600 km 2 (40 km × 40 km), to a depth of 3.5 km, containing 1.7 million grid cells (133 × 188 × 69). Initial estimates for porosity and permeability for each cell were determined from the operator's in-house geological model. An anisotropic permeability field is used to represent the effects of aligned sets of fractures, with maximum horizontal permeability aligned with the maximum stress direction.
Four components were modelled: H 2 O, NaCl, CO 2 and CH 4 . CO 2 is able to dissolve into the water. Thermal effects were modelled using an energy transport equation. The fluid flow boundary conditions were no-flow boundaries at the edges of the model. The reservoir was simulated from the onset of injection in 2004 to the present, and forward modelled out to the year 3000, outputting pore pressure maps at monthly intervals.
The model was history matched by comparing observed and calculated bottom hole pressures (BHP) for each well. BHPs for the field were estimated using commercial well performance code (PROSPER, Petex Edinburgh). BHP is determined through integration of fluid density down the borehole with an additional correction to account for friction in the borehole. The density of the CO 2 in the borehole is determined as a function of temperature and pressure, which are both calculated down the borehole. The flow simulation model calculates BHP at each well. The model was adjusted by altering the porosities and permeabilities of sets of grid cells until there was an acceptable agreement between calculated BHP and those determined by the well performance code. A good match was obtained. The above workflow represents standard industry practice in history matching reservoir flow models with observed pressures and flow rates from the field.

Modelling initial stress conditions
Initial work on this project used a coupled fluid-flow/geomechanical algorithm provided by the operator (BP) to fully simulate changes in effective stress caused not just by pore pressure changes but also by the poroelastic expansion of the reservoir during injection. However, this model only covered the earlier injection stages (2004)(2005)(2006)(2007)(2008), and provided output only at quarterly intervals.
Our analysis of this coupled model indicated that in the reservoir the poroelastic coupling effects were not significant, and the changes in effective stress could be reasonably approximated by subtracting the changing pore pressure from an initial stress condition. We take this approach in our work, modelling initial stress conditions using a geomechanical solver, and using pore pressures computed by the reservoir flow model to compute the effective stress at monthly intervals.
The geomechanical solver to determine the initial stress conditions uses the same mesh as the fluid flow solver. The model boundary conditions were no-strain boundaries at the edges and base of the model, while the ground surface is modelled as a free surface. The initial vertical stress gradient was based on integration of density logs from the field. Horizontal stress gradients were determined by fitting modelled initial stresses to a variety of geomechanical observations from the field, including leak-off tests and borehole breakouts. The vertical stress gradient was taken as 22.5 kPa/m, the maximum horizontal stress gradient was 25.1 kPa/m, and the minimum horizontal stress gradient was 15.8 kPa/m. The orientation of the maximum horizontal stress direction was 135 • (SE).

Stress changes on fractures and modelling induced seismicity
Each modelled fracture takes its initial stress conditions, σ f ij , from the applied stress at the nearest node of the reservoir model, σ nn ij . Langenbruch and Shapiro (2013) show how natural material heterogeneity leads to variability in both principal stress magnitude and orientation when uniform external stresses are applied. Therefore, the stresses from the reservoir model are modulated for each fracture. Stress modulation factors, S mod i , and angles, θ mod i (i = 1, 2, 3) are defined for every fracture, with S mod i chosen at random from a normal distribution with mean of 0 and standard deviation of 10%, and θ mod i chosen at random from a normal distribution with mean of 0 and standard deviation of 5 • . The magnitude of each principal stress is modulated by S mod after which σ fm ii is rotated by Euler angles θ mod cos θ 1 cos θ 2 − cos θ 3 sin θ 2 sin θ 1 cos θ 1 sin θ 2 + cos θ 3 cos θ 2 sin θ 1 cos θ 1 sin θ 3 − sin θ 1 cos θ 2 − cos θ 3 sin θ 2 cos θ 1 − sin θ 1 sin θ 2 − cos θ 3 cos θ 2 cos θ 1 cos θ 1 sin θ 3 to give σ fr ii , the rotated, modulated stress tensor acting on a fracture.
The reservoir model provides a map of pore pressure at monthly time intervals. We compute the effective stress acting on each fracture, σ e ij from σ fr ij and the pore pressure at the nearest node of the fluid flow model, P nn : where δ ij is the Kronecker delta. Having computed the effective stress for each individual fracture, we resolve this stress into shear and normal stresses, τ f and σ f n , acting on the fracture face: where n is a unit vector normal to the fracture. The Mohr-Coulomb properties of each fracture must be defined to determine whether τ f and σ f n are sufficient to induce failure. In the absence of data from the In Salah field, we use generic values in our model. For each fracture the cohesion and friction coefficient, C f and μ f , are chosen at random from a normal distribution, with mean of 2.2 MPa and standard deviation of 0.5 MPa for C f , and mean of 0.6 and standard deviation of 0.1 for μ f . If stress conditions are such that the Mohr-Coulomb envelope, τ 0 > μ f σ f n + C f , is exceeded, then an event is deemed to have occurred.

Event source characteristics
When events have been triggered in our model, we compute event source characteristics for comparison with observation. Event magnitudes can be determined from the area of rupture, A, and the stress drop, σ (Kanamori and Brodsky, 2004): To model event magnitudes we determine A from the length of the triggered fracture, multiplied by 20 m, the thickness of the In Salah reservoir. We limit fracture heights to the height of the reservoir to account for the fact that (1) the differences in mechanical properties between layers are likely to act as a barrier to fracture propagation, and (2) the stress conditions in the overburden are likely to be different, so although a seed may be triggered in the reservoir, Mohr-Coulomb criteria may not be exceeded in the overburden.
We also account for the fact that the full area of the fracture may not rupture in a given event, multiplying A by a modifying factor chosen at random from a uniform distribution ranging from 1-100%. The stress drop is taken from τ f , the shear stress acting on the fracture. An event may only release a portion of the accumulated shear stress, so σ is taken as τ f multiplied by a modifying factor chosen at random from a uniform distribution ranging from 1-100%.
We assume double-coupled source mechanisms for our events, defined by the strike and dip of the nodal plane, and the rake, which describes the slip direction along this plane. For each event the strike and dip are taken from the fracture plane that has triggered. The rake is determined from the orientation of the shear stress vector on the fracture plane.

Stress drops and repeated events
The occurrence of an event will serve to release stress on a fracture, reducing the likelihood of further seismicity on this plane. This stress drop must be accounted for at subsequent time-steps after an event has occurred (e.g. Goertz-Allmann and Wiemer, 2013). For fractures that have experienced failure in a previous time-step, we rotate σ e ij into a coordinate system defined by the shear and normal stress vectors associated with that event, before subtracting σ from the off-diagonal component of the rotated stress tensor. Once the stress released by the previous event has been removed, the updated stress tensor is returned to the global coordinate system, and equations (6)-(8) are used as above to determine whether the fracture has again exceeded the Mohr-Coulomb threshold, triggering a repeat event. This means that a single point has the potential to slip multiple times during the simulation if pore pressure continues to increase.
The Mohr-Coulomb parameters μ and C are chosen from random distributions. This means that there is a chance that some fractures already exceed the Mohr-Coulomb limits at the start of the simulation. We interpret these as fractures that would have slipped and released stress naturally during the geologic history of the field. We therefore perform an initiation step, identifying any fractures that exceed Mohr-Coulomb criteria at the initial stress conditions, and removing all shear stress on these fractures following the method outlined above.
The process is probabilistic in nature, as the various parameters described above are selected at random from statistical distributions, meaning that multiple realisations must be considered to ensure statistical robustness. After 10 simulations we noted that the results were relatively consistent between simulations, removing the need for such extended modelling. In the following section we compare our simulated event populations with observations made at In Salah by Stork et al. (2015).

Results and comparison with observations
Microseismic monitoring at In Salah was conducted from a single shallow borehole near to injector KB-502. This monitoring array was switched on in August 2009, monitoring until the end of injection in 2011. Our model covers the whole of the field, including all three injection wells and the 4 gas producers. Many of the modelled events occur in proximity to injection wells KB-501 and KB-502. However we focus our attention on modelled events located within 5 km of the microseismic monitoring well (i.e. in proximity to well KB-502), occurring after August 2009, as these can be directly compared with observations from the field. Fig. 3a shows the observed rate of seismicity at In Salah. After a pause of approximately 2 yr, injection in well KB-502 recommenced in November 2009. Injection rates and pressures increased during the first few months of 2010, and small amounts of seismicity were detected. A significant increase in the rate of seismicity was observed in April 2010, which peaked in June, returning to a lower rate by August 2010. Fig. 3b shows the modelled rate of seismicity. We note good agreement between changes in observed and modelled seismicity rates through time. Our model predicts a limited degree of seismicity as injection recommences in November 2009. This rate increases significantly in April 2010, peaking in June and returning to the lower rate after August 2010. The increase in seismicity between April and August 2010 is associated with an increase in injection pressure. However, it is notable that injection pressures reach these levels in February and March 2010, without a notable increase in seismicity. The observed behaviour -the delay between injection pressure increase and the increase in seismicity -is captured by our model.

Event occurrence rate
The major discrepancy between model and observation lies in the total number of events. While over 6000 events were robustly identified at In Salah, our 10 models simulated a mean of 517 events (with a minimum of 491 and a maximum of 555). We examine potential reasons for this discrepancy in the Discussion section. Fig. 4 shows the locations of modelled events from a single model realisation. Technical issues at In Salah left only a single geophone providing useful data. As such, it is not possible to determine accurate locations for observed events to be compared with our modelled locations. Instead, Stork et al. (2015) constrain approximate locations for In Salah events by performing a hodogram analysis to determine P -wave arrival azimuths, θ P , and inclinations, ϕ P , as well as measuring the differential P -to S-wave arrival times, t S-P (Fig. 5a and c). Stork et al. (2015) observe two main event clusters. The largest, by number of events, had t S−P ≈ 0.65 s, arrival azimuths trending NW-SE, and inclinations from 0 • (i.e. directly underneath the geophone array) to 20 • . A second cluster was observed at greater distance from the monitoring array with t S−P ≈ 0.95 s, located to the NW with θ P ≈ 300 • , and ϕ P ≈ 25 • .

Event locations
In order to make comparisons between these observations and our model results, we must simulate arrival angles and t S-P for our modelled event locations. We use a finite difference waveform simulation tool E3D (Larsen and Schultz, 1995), performing 2D simulations to model waveforms from each event to the single recording geophone. We then determine t S-P and ϕ P from these waveforms. Because the simulations are performed in 2D, using a 1D layer-cake velocity model, values for θ P are taken from the geographic azimuth between event location and geophone. The modelled arrival parameters are shown in Fig. 5b and d.
We note that modelled t S-P and ϕ P are very sensitive to the choice of velocity model, which at In Salah is based solely upon sonic log data. Standard practice when using microseismic monitoring is to calibrate the velocity model with perforation shots, whose locations are known. Because this data was not collected at In Salah, the velocity model is not calibrated.
Nevertheless, we observe broad agreement with observed modelled parameters. In both cases the trend of events azimuths is NW-SE: the mean observed θ P was 110 • , while the mean modelled θ P was 111 • . 0.65 < t S-P < 1.0 s for both observed and modelled populations, with t S-P increasing with incidence angle.
The majority of events with t S-P ≈ 1.0 s are to the NW of the monitoring well, but there are also a handful of such events to the SSE as well.
The most notable discrepancy between observations and our model results lies in the clustering of observed travel times (5c and d). Observed t S-P times have a distinctly bimodal population, mostly t S-P ≈ 0.65 s, with a smaller grouping at 1.0 s. In contrast, our modelled t S-P times show a more continuous distribution between 0.6-1.0 s. We are not able to determine from where this discrepancy has been caused. Regardless, the overall agreement that we observe between event arrival azimuths and travel times, is notable.

Stork et al. (2015) compute magnitudes for events observed at
In Salah (Fig. 6a). The largest event had a magnitude of M W = 1.7. The Gutenberg-Richter b-value for the whole population was 2.17 ± 0.09, and the minimum magnitude of completeness was M W = 0.05. Fig. 6b shows the frequency-magnitude distribution for our modelled events. The mean largest event for our 10 model realisations was M W = 1.73, with a minimum of 1.61 and maximum of 1.82. The mean b value was 2.3, with a minimum of 2.1 and a maximum of 2.7. Fig. 7 shows the strike, dip and rake of our modelled focal mechanisms. The population is dominated by left-and right-lateral strike slip behaviour (rake ≈0 • or ≈180 • ) on subvertical fractures with strikes ranging from 90 • -180 • , with a mean of 135 • (NW-SE). The limitations of the recording geometry at In Salah means that focal mechanisms for the observed events cannot be determined, so we are not able to compare our modelled source mechanisms with observation.
In Fig. 8 we examine stress drops and scaling relationships for our modelled event population. Stress drops range from 0.01-4 MPa, with a mean of 1.5 MPa. We note that our modelled events follow scaling relationships expected of small-magnitude earthquake populations (e.g. Abercrombie, 1995). Again, due to the limitations of the recording array, stress drop calculations were not possible for the In Salah events (e.g. Stork et al., 2014), so we are not in this case able to directly compare our model results with observations.

Discrepancy between number of modelled and observed events
The largest discrepancy between our model results and observed events at In Salah is the number of events. While the relative rates of seismicity are in good agreement, our models predict only ∼500 events during the period of interest, while over 6000 events were robustly identified at In Salah.

Number of seed points
Perhaps the simplest way to increase the number of modelled events is to increase the number of seed points. The number of seed fractures is ultimately constrained by the strain modelled by our geomechanical reconstruction. However, the smallest fracture length was arbitrarily constrained at 50 m. The modelled strain could instead be modelled by a larger number of smaller fractures.
This would increase the number of modelled events, and especially events of smaller magnitude. From the frequency-magnitude plots shown in Fig. 6, the number and b value for larger events is well-matched with observations, but the modelled distribution suffers from a lack of smaller events. Accommodating the modelled strain with a larger number of smaller fractures would address this discrepancy, and is therefore a plausible explanation for our under-prediction of the overall number of events. Nevertheless, in the following sections we consider other mechanisms that may result in a smaller number of modelled events than observed.

Mohr-Coulomb property distribution
The number of modelled events is sensitive to the Mohr-Coulomb failure properties of the fractures. We investigated the sensitivity of the number of events to these properties, performing a simulation with the cohesion set to 0 MPa for all seed points, and a reduced mean μ f of 0.4. This increased the number of events, but only to a total of approximately 1000, which still falls far short of the number of observed events.   9. In (a) we plot CFS contours (in MPa) generated by the largest modelled event. In (b) we project the modelled stress changes onto surrounding seed points (the event is denoted by the purple star). In (c) we plot the changes in pore pressure at each seed point from initial conditions to May 2010. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Temporal resolution of input data
We note that our model consists of pore pressures updated on a monthly basis. In fact, injection rates and pore pressures are changing continually at In Salah (e.g. Fig. 3). These variations are not included in our model because modelled pore pressures were only available on a monthly basis.
In our modelling approach, events are triggered by changes in pore pressure, which cause changes in the effective stress. Therefore the total number of events may be modulated by the temporal resolution of our input pore pressure data. Finer temporal resolution may capture shorter-term pore pressure changes capable of triggering seismicity that are not captured by monthly average data, leading to an underestimation of the total number of events.
We investigated the potential impact of the temporal resolution issue by re-running our simulations with lower temporal resolution. We decimated our input data to a quarterly (3 month) input. In doing so, the number of events was reduced to approximately 450, a small reduction of 14% fewer than the mean value for our original models. However, the difference between modelled and observed numbers of events is of the order of 1000%, suggesting that improvement of temporal resolution alone might not be sufficient to account for the differences between observed and mod-elled numbers of events, although without finer resolution data we cannot confirm this directly.

Fracture-fracture interaction
An assumption implicit in our approach is that all events are triggered by pore pressure changes, and none by interactions between fractures. In reality, a mix of triggering mechanisms would be expected, and Sumy et al. (2014) have shown evidence of events induced by waste water disposal being initially triggered by pore pressure changes, but subsequently promoted by the static stress changes generated by the initial events.
In order to test the relative impact of pore pressure changes and static stress transfer, we compute the static stress change produced by the largest modelled event, and compare it to pore pressure changes. We model the Coulomb failure stress (CFS) generated by our largest event using a standard approach (Wang et al., 2006). This event has a length of 136 m, a stress drop of 3.7 MPa, and strike/dip/rake of 101 • /85 • /175 • . We resolve the CFS on to a generic seed point with a strike/dip of 140 • /85 • from horizontal. Fig. 9a shows the resulting changes in CFS produced by this event. In Fig. 9b we show how this CFS would impact nearby seed points. For comparison, in Fig. 9c we show the pore pres-sure change from initial conditions to May 2010 (a period with elevated injection pressure). We note that the change in pore pressure reaches 10 MPa at the injection well, and is approximately 4 MPa at the site of the largest event (which is approximately 1 km from the injection well). The change in CFS reaches ±5 MPa within half a source length of the event, but soon decays to <0.5 MPa at distances greater than a single source length (136 m).
From this we conclude that static stress transfer between events may rival the impact of pore pressure changes only within a very limited distance from the source point. Across the majority of modelled seed points, pore pressure effects appear to dominate in our analysis, even for the largest of the modelled events. This provides a justification for neglecting fracture-fracture interactions in our analysis.
In comparison with other injection-induced events (e.g. Verdon, 2014), seismicity at In Salah is not of particularly large magnitude, while the pore pressure changes induced by injection are relatively large in comparison to other injection projects (e.g. Verdon et al., 2013). In situations where the pore pressure change induced by injection is smaller, but the induced events are larger, then it may be that static stress transfer becomes more important, as found by Sumy et al. (2014).
In the light of our observations here in combination with those made by Sumy et al. (2014), we suggest that future research should consider combining both triggering by pore pressure and through static stress changes to fully understand seismicity triggered by fluid injection.

Propagation of new fractures by hydrofracture
Our modelling assumes that all seismicity occurs on the preexisting fractures modelled by Bond et al. (2013). However, it is possible, even likely, that new fractures have been created during the injection process. Goertz-Allmann et al. (2014) examine injection rates and pressures, noting two different behaviours, referred to as "matrix injection" at lower injection pressures, where fluid enters the rock matrix without altering the permeability, and "fracture injection" at higher injection pressures, where permeability increases with injection, implying that fractures are being opened or created. This in itself does not necessarily imply that new fractures are being created, as it is possible that pre-existing fractures are being opened without the creation of new hydrofractures. However, White et al. (2014) examine bottom-hole pressures (BHP) in detail, concluding that at times BHP in well KB-502 exceeded the fracture initiation pressure, and therefore that the creation of new fractures by hydrofracturing is likely.
The tensile failure that accompanies fracture propagation during hydraulic stimulation typically produces smaller magnitude microseismic events than the accompanying shear failure. Indeed, tensile events often go unobserved during hydraulic stimulation (e.g. Eisner et al., 2010;Maxwell, 2010). Nevertheless, the propagation of new fractures will produce new seed points on which detectable events can occur, increasing the number of events. White et al. (2014) conclude that the relative importance of pre-existing fractures versus newly created fractures at In Salah is a large uncertainty: without better constraint in this regard, it is impossible to determine the extent to which this effect contributes to an underestimation in the number of modelled events. However, using shear-wave splitting analysis, Stork et al. (2015) found that increases in anisotropy (and by inference fracture intensity) during elevated injection periods reduced back to original values when the pressure dropped, implying that the dominant mode of deformation during 2010 was the opening of pre-existing fractures, rather than the creation of new fractures.

Out of zone fractures
The limitations of the input data used in this study means that our modelling approach can only simulate stress changes on fractures within the reservoir unit. The stability of observed t S-P times implies that most of the seismicity occurs in and around the reservoir. However, it is clear both from the increased variation in observed t S-P during high injection periods, and from inversions of InSAR surface deformation measurements, that some injectioninduced deformation extends into the caprock. By extending our modelling to include fractures and stress changes in the overburden, we would undoubtedly increase the number of modelled events.

Events at injection wells KB-501
In order to compare our model with observations, we have focused our analysis on modelled events located close to where the microseismic monitoring array was installed. However, our model also predicts substantial numbers of events around injection well KB-501, which is approximately 10 km to the southeast of the microseismic monitoring well. Stork et al. (2015) detected 11 events that, based on them having larger t S-P times and azimuths of ∼135 • might be linked to well KB-501.
The fact that our model predicts a substantial number of events around these wells, but only a few are observed, implies that they may be too far from the monitoring well to be detected. We note that high velocity layers are present in the overburden at In Salah. Stork et al. (2015) traced rays through this velocity model, finding that for oblique take-off angles the waves became trapped by these layers and did not reach the surface at larger distances. This implies that in addition to geometric spreading and intrinsic attenuation, the amplitude of events near to KB-501 may be further reduced on detection at the microseismic monitoring well.
We use reflectivity code based on Kennett (1993) to model the change in amplitude as a function of distance for an event at 1.85 km depth for the layered In Salah velocity model. We find that the amplitude of event recorded at 10 km distance from the monitoring array is approximately 1.8% of an event located underneath the array. The largest event recorded at In Salah had a recorded amplitude of 1150 nm. The same size event occurring at 10 km distance would therefore be expected to have an amplitude of 20 nm at the array. We estimate the magnitude of completeness at In Salah to be 0.05 (Fig. 6). To establish the amplitude that results in event detection, we therefore examine the recorded amplitudes of observed events with magnitudes ranging from 0.0 to 0.1. We find that these events range in amplitude from 3 nm to 15 nm.
The modelled reduction in amplitude as a function of distance is strongly dependent on the velocity model used, as well as other factors such as attenuation. Therefore the above discussion is only intended to provide an approximate estimation of detection thresholds. Nevertheless, we find that the expected amplitude for the largest event at 10 km distance produces a similar amplitude (20 nm cf. 15 nm) to an event whose magnitude is at the detection threshold located directly below the monitoring array. The observations are therefore not inconsistent with model results: it is plausible to suggest that events have occurred at KB-501 but gone undetected due to the distance from KB-501 to the monitoring well. For obvious reasons, this assertion cannot be confirmed with existing data: a more comprehensive monitoring program would have been required to robustly assess seismicity, if any, at injection well KB-501.

Implications for geomechanical behaviour at In Salah
White et al. (2014) consider the observed geomechanical behaviour at In Salah, focusing principally on the observed surface uplift, pressure changes, and seismic reflection imaging, and posit a number of hypotheses: (1) that all behaviour can be explained by deformation in the reservoir; (2) that a fault intersects well KB-502 and is leaking CO 2 into the overburden; (3) that injection caused hydraulic fracturing in the reservoir and lower caprock units; or (4) that pre-existing fractures have been re-activated by injection. White et al. (2014) conclude that hypothesis (3), that injection has caused hydraulic fracturing in the reservoir and lower caprock units, is the most probable, but remain uncertain as to the contribution from pre-existing fracture networks.
In this paper we model the position and orientation of these pre-existing fracture networks to simulate microseismic activity, finding that our simulations provide a good match to field observations. The fact that we are able to do this indicates that the pre-existing fractures are likely to have played a significant role in the deformation induced by CO 2 injection at In Salah.
The caprock at In Salah is 950 m thick, comprising a number of thick, resilient seals. White et al. (2014) conclude that, even though fracturing may penetrate the lower caprock units, the overall integrity of the storage site is maintained. Both reflection seismic imaging and inversion of InSAR surface deformation measurements indicate that deformation extends no more than 100-200 m into the overburden. Similarly, while event locations are uncertain, there is no evidence from the microseismic data that fractures have propagated a significant distance into the caprock. Therefore our conclusions match those made by White et al. (2014): the thickness of the seal means that it is unlikely that injection-induced deformation has compromised storage integrity.

Conclusions
Workflows that model induced seismicity are required to link geomechanical simulations with field observations. The benefits from doing so are twofold -as a calibration, demonstrating that models are matching observations from the field; and as an interpretive step, to help understand why microseismic events occurred where and when they did from a geomechanical perspective.
We have developed a workflow using geomechanical models to simulate injection-induced seismicity, applying this workflow to the In Salah CCS project, Algeria. Our method is based on the assumption that seismicity occurs on pre-existing fractures. Thus, if the location, orientation and size of fractures is known, or can be modelled, and if the stress changes caused by injection can be modelled, then by resolving stress changes into changes in normal and shear stress on these fractures we can determine where and when seismicity is likely to occur, as well as event magnitudes and source mechanisms.
We compare our modelled events with observations from In Salah, noting good agreement on a number of counts. Our model does a good job of capturing changes in the rate of seismicity through time, as well as the locations of events. Our model also performs well with respect to the magnitudes of induced events, with the largest modelled event closely matching the largest event observed at In Salah.
It is clear that geomechanical deformation can pose a risk to CCS projects. This risk is two-fold: geomechanical deformation can compromise the integrity of the caprock; and it can trigger seismic events. Potential CCS reservoirs should be carefully appraised from a geomechanical perspective, and models created to simulate the impact of CO 2 injection.
Furthermore, geomechanical deformation should be monitored once injection begins. Geomechanical monitoring at CCS sites has primarily consisted of geodetic methods (InSAR, GPS and tiltmeters) and microseismic monitoring. Initial models must be updated so that they provide a good agreement with field observations as injection continues. Our workflow, which requires the kind of careful appraisal described above as input parameters, is capable of simulating induced seismicity, allowing comparison between geomechanical model and field appraisal.