- Email: [email protected]
Microseismic assessment and fault characterization at the Sulcis (South-Western Sardinia) field laboratory
International Journal of Greenhouse Gas Control 95 (2020) 102974
Contents lists available at ScienceDirect
International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc
Microseismic assessment and fault characterization at the Sulcis (SouthWestern Sardinia) field laboratory
T
M. Anselmia,*, G. Saccorottib,1, D. Piccininib,1, C. Giunchib,1, M. Paratoreb,1, P. De Goria, M. Buttinellia, E. Maggioc,2, A. Plaisantc,2, C. Chiarabbaa a
Istituto Nazionale di Geofisica e Vulcanologia, Rome, Italy Istituto Nazionale di Geofisica e Vulcanologia, Pisa, Italy c Sotacarbo S.p.A., Carbonia (SU), Italy b
ARTICLE INFO
ABSTRACT
Keywords: CCUS Fault Lab Seismic baseline Ambient noise tomography Sardinia Sulcis basin
The general acceptance of the CO2 geological storage by stakeholders passes through the assessment and mitigation of risks, potentially induced or increased by the disposal activity. Injection of moderate to large quantities of CO2 in the sub-surface may unbalance local stress and trigger earthquakes if faults are critically stressed, condition that is not easily verifiable. Pilot sites are therefore the best way to proceed further in order to address such challenging issues. In such cases, the reconnaissance of faults and seismicity in the sub-surface, before the onset of activity, is mandatory. In this paper, we present studies carried out in the site where the Sotacarbo Fault Lab is going to be installed. This facility will be located in a very low seismic hazard region of central Mediterranean, where reports on historical large earthquakes are poor. We show results from a series of experiments aimed to monitor the background seismicity around the pilot site. As expected, seismicity is almost absent down to small magnitude close to the future injection-test well. Further seismic imaging of the subsurface layers obtained by ambient noise tomography offers the ability to resolve the presence of a seismicityfree fault located in the first 200 m below the surface, of which the last episode of activity is difficult to assess. Our results encourage the use of this site to follow the response of the system to injection of small quantity of CO2.
1. Introduction The crucial need to reduce the CO2 emission requires a multiplicity of efforts, among which carbon capture, utilization and storage (CCUS) is a viable opportunity (Dipietro et al., 2011; Boot-Handford et al., 2014; Liu et al., 2017). Geological storage of CO2 might produce environmental and safety hazards, (e.g., accidental leakage during the injection at critical pressure, slow leakage after injection in the reservoir following diffusion and upwelling processes in presence of cracking/faulting of the cap rocks) that need to be reliably assessed (Nicol et al., 2011; Vilarrasa et al., 2019). Geo-energy activities that involve fluid injection/extraction may change the state of stress in the subsurface causing deformation, microcracking and fault reactivation (Grigoli et al., 2017). Induced seismicity has received increasing attention by stakeholders (Ellsworth, 2013) and is potentially one of the main obstacles to the volumetric upscale of geo-energy projects.
This is particularly true for areas where the natural seismic hazard is low like in the Sulcis area (Gruppo di Lavoro, 2004; Woessner et al., 2015) and therefore the existing infrastructure may not be adequate to withstand the impact of potentially large anthropogenic seismicity (Takagishi et al., 2014; van Thienen-Visser and Breunese, 2015). The capacity of storing large quantities of CO2 scales with the volume that will be affected by changes in pressure and stress, increasing the risk of intercepting active faults and inducing large events (Keranen et al., 2014). Different views have emerged about up-scaling storage projects and their applicability or convenience. In particular, some opinions emerged proposing that large-scale CCUS is risky and possibly unsuccessful (Zoback and Gorelick, 2012), for the high probability that earthquakes could violate the seal integrity. Conversely, other opinions assessed that CO2 storage can be performed safely, after proper site characterization and with adequate pressure management (Vilarrasa and Carrera, 2015; Vilarrasa et al., 2019). At present, pilot and
Corresponding author. E-mail address: [email protected] (M. Anselmi). 1 Saccorotti, Piccinini, Giunchi, Paratore belong to the branch of Pisa Istituto Nazionale di Geofisica e Vulcanologia, Via Cesare Battisti, 43-56125 Pisa-Italy. 2 E. Maggio and A. Plaisant belong to the organization Sotacarbo S.p.A - Grande Miniera di Serbariu, 09013 Carbonia (SU) - Italy. ⁎
https://doi.org/10.1016/j.ijggc.2020.102974 Received 6 August 2019; Received in revised form 15 January 2020; Accepted 17 January 2020 Available online 25 January 2020 1750-5836/ © 2020 Elsevier Ltd. All rights reserved.
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
demonstration projects focus on verifying the storing capacity and the development of monitoring best practices. Considering the acceptance of a reasonable level of seismic risk, the general goal is to chart a course and defines clear and plain protocols for safe and efficient monitoring that afford reliable and accurate hazard forecasts in a real-time management of the risk (Bommer et al., 2015). One of the main challenges is to define the existence of faults capable of producing damaging earthquakes. Since activities are generally shallow, even small to moderate magnitude earthquakes can produce significant ground shaking at the surface. In addition, faults that could fail with small to moderate earthquakes are more difficult to be detected in the subsurface. Earthquakes induced by fluid injection often develop on faults barely visible in seismic imaging of the subsurface (Buttinelli et al., 2016). The reliable identification and mapping of aseismic, but possibly critically stressed faults is therefore a first key step in hazard assessment, as became clear in unconventional geothermal systems, like the Enhanced Geothermal Systems (EGS, see Zang et al., 2014). Moreover, faults can act as enhancer or obstacle to fluid flow and their reconnaissance and characterization are important to plan efficient CO2 sequestration. In this paper, we investigate such aspects for the Sotacarbo Fault Lab (South-Western Sardinia, Italy), a facility specifically designed as a test bed infrastructure in order to develop and test in-field conditions instruments and methods to be applied in the commercial-scale CO2 geological storage site elsewhere. The facility will be located in SouthWestern Sardinia, a region with low natural seismic hazard (Stucchi et al., 2011; Rovida et al., 2016). According to that, the national permanent seismic monitoring system has an intrinsic high magnitude threshold for detection of local earthquakes in Sardinia, where the Magnitude of completeness (Mc) is higher than 2.9 (Schorlemmer et al., 2010). To improve earthquake detection, various dense seismic arrays have been deployed over the past 4 years. The Sulcis temporary seismic network was installed for a period of about 15 months since July 2014 in order to monitor the South-Western Sardinia (Fig. 1a). After that, in July 2016 we installed the Sotacarbo Fault Lab permanent seismic network of 5 stations around the fault lab, complemented with a temporary network consisting of additional 11 stations operative from the first days of June 2017 to end of August 2017 (Fig. 1b). We present results of the multiscale survey and quality assessment of the overall monitoring efficiency. In order to enhance the imaging of the fault in the lab site, we performed a focused experiment for ambient noise tomography, a very low cost and efficient tool for subsurface imaging especially in case of absence of local seismicity (Zheng et al., 2008). During the first days of October 2017, we installed a temporary linear seismic array for about 4 days (see line in Fig. 1b), while the Sotacarbo Fault Lab permanent seismic network was operating. Data from both the passive seismic linear array and the Sotacarbo permanent and temporary seismic networks enable the tomographic reconstruction of the fault.
monitoring system that will operate during the future activities of the Sotacarbo Fault Lab. The South-Western Sardinia exhibits a very low seismicity level; the historical catalog reports only two earthquakes occurred in 1616 and 1771, with estimated moment magnitudes of 4.9 and 4.4, respectively (Rovida et al., 2016; Fig. 1). The low seismicity rate of the area is also confirmed by the instrumental catalog issued by Istituto Nazionale di Geofisica e Vulcanologia (INGV). For the time interval 1985-present, it reports only five earthquakes, among which the ML 4.1 occurred in 2006, a few kilometers offshore the Gulf of Cagliari (Fig. 1a). In order to increase the current seismic sensitivity of the national network toward the detection of low-magnitude earthquakes, we installed the first large temporary network of seismic stations in Sulcis, composed of high-dynamic, local-recording digitizers Reftek 130-1. The seismic network was composed of 8 Lennartz Le3D/5 s sensors and 2 Nanometrics Trillium Compact 120 s (see Fig. 1a). The network operated in continuous – recording mode between July 2014 and September 2015. The main requisite for seismic baseline is to capture the smallest magnitude occurring in the area. We therefore started analyzing the level of possible detection though the evaluation of the noise condition at each individual site, by using the Probability Density Function (PDF) of the power spectra evaluated for the three components of ground velocity. In addition, local amplification phenomena were investigated using the Horizontal-to-Vertical spectral ratios (Fig. 2). To estimate the detection capability of the temporary network, we implemented from scratch a procedure based on the comparison of the observed noise spectra with those predicted for a point source of known location and magnitude. The noise spectral amplitude was obtained by computing average spectral density for 1 week of seismic signal between 1 and 2 a.m. for the nighttime measures, and between 1 and 2 p.m. for daytime measures. First, we parameterize the Earth’s volume below the network using a regular grid of nodes spaced by 1 km along the NS and EW directions, and 5 km along the vertical direction. Each grid node represents a virtual seismic point source defined using a Brune’s model [Brune, 1970] at which we assign a moment magnitude Mw. Then we derive the scalar seismic moment and the expected displacement spectrum (following Kanamori, 1977 and Boatwright, 1978) and then we propagate the seismic amplitude in a 5.8 km/s half-space velocity model with homogeneous frequency independent attenuation (Q = 250). For a given source location and magnitude, the ground motion spectrum predicted at any individual station (taking into account the topographic differences between the stations) is then compared to the average noise spectrum observed at that particular site. If the ratio between the predicted and noise spectra is larger than 2 over the 3−25 Hz frequency band, than the procedure assumes that that station is able to detect an earthquake of that particular magnitude. When the number of stations detecting the earthquake is greater than 4, that magnitude value is accepted and attributed to the grid node under evaluation. The procedure is iterated over the entire grid of virtual sources, for moment magnitudes Mw spanning the 0.1–2.2 range, then converted in local magnitudes ML according to Munafò et al. (2016). This allows deriving maps of the minimum detectable magnitudes as a function of source location at the hypocentral depth of 5 km b.s.l.. We took into account two different scenarios, at day and night time. The estimated theoretical detection level within the Sotacarbo Fault Lab is around the ML -1.0 (see Fig. 3). This is a considerable improvement, considering the initial high detection threshold in the area ensured by the INGV national network.
2. Baseline monitoring of the Sulcis site The Sulcis is a mining area of South-Western Sardinia, where a former coalfield within the Matzaccara Quaternary basin has been selected as the site of the Sotacarbo Fault Lab (Fig. 1a and b). Aim of this facility is to study CO2 leakage processes through faults and develop advanced monitoring tools to be applied in commercial-scale CO2 geological storage applications worldwide. The characterization of the site in terms of adequateness with respect to natural hazards is therefore a pre-requisite. Following the general framework of the technical and geological feasibility evaluation of CCUS (Arts et al., 2004; Förster et al., 2006; Ringrose et al., 2013), the definition of a seismic baseline is an important aspect (Ugalde et al., 2013; Chen et al., 2015). Seismic baselines are crucial to define the reference level of natural seismicity before the onset of activities and provide hints about the presence of active faults. Data acquired during the baseline aid to design the seismic
3. Results of seismic monitoring and earthquakes location From inspection of one year of continuous recordings from the firstlarge temporary network (red and blues triangles in Fig. 1a), we 2
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 1. a) Map of the study area in the South-Western Sardinia. Red triangles represent the stations of the first-large temporary network. The blue triangles represent the broadband (BB) stations of the same network. Dashed black square indicates the location of the Sotacarbo Fault Lab (see Fig. 1b). Historical earthquakes for South-Western Sardinia are marked by purple stars. Instrumental earthquakes are marked by brown circles. b) Map of the Sotacarbo Fault Lab area. The cyan and red triangles represent the stations of the permanent network used to monitoring the Sotacarbo Fault Lab area. Green triangles represent the Sotacarbo temporary seismic network added for the noise tomographic experiment. The blue circles represent the passive linear seismic array deployed across the presumed strike of the fault zone.
3
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 2. PDF (Probability Density Function) analysis, calculated using 30 days of the ground motion continuous recordings. From top to bottom, the graphs represent the PDF of EW, NS and Vertical components. The red line indicates the 95 percentiles of the probability density function. The density of probability is indicated by the grey shaded palette. The solid lines represent the upper and lower models of the Earth’s noise after Peterson (1993). Right panels: HV spectral ratios and corresponding ± 3σ boundaries, for data samples collected at day- and night times (top and bottom panels, respectively).
obtained a dataset of 82 earthquakes, which were located using the 1-D velocity model of the INGV seismic monitoring system (Table 1) and Hypoellipse code (Lahr, 1999). Only 5 hypocenters have a high-quality location, r.m.s. less than 0.19 s., with horizontal and vertical errors < 1 and 2 km, respectively, number of stations ≥ 5 and azimuthal gap less
than 180°. Fig. 4 shows these earthquakes hypocenters (black circles), located both in the central part of the study area (North-East of Carbonia town) and in the Sant’Antioco island, South-Eastern off-shore; hypocentral depths range from 2.5 to 18 km (Table 2). The green circles represent the remaining 77 earthquakes’ hypocenters with location
Fig. 3. Sensitivity maps at day- (left) and night-time (right) for the Sulcis temporary network, assuming a hypothetical hypocentral depth of 5 km. The dashed black box indicates the Sotacarbo Fault Lab area. 4
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
the Matzaccara Quaternary basin (a part of the wider Sulcis basin) where the Sotacarbo Fault Lab is going to be installed. In May 2016, we installed the Sotacarbo Fault Lab permanent seismic network that consists of 5 seismic stations, 4 equipped with three components intermediate period (Lennartz 3D/5 s) seismic sensors and 1 with a broadband sensor (Nanometrics Trillium compact 120 s), operating in the restricted area around the site (cyan and red triangles in Fig. 1b). In Fig. 5 the sensitivity map of the Sotacarbo permanent seismic network is shown. During a three-month period in 2017, up to 11 additional seismic stations were operating to temporarily obtain a denser temporary sampling in the Sotacarbo Fault Lab site. During the recording period of dense monitoring (60 days), seismicity was not detected in the study area and its close surroundings.
Table 1 1-D velocity model used for earthquakes’ location shown in Fig. 4. This model is the mean velocity model of the Italian territory, used by INGV during the realtime seismic survey of the Italian territory (Bollettino Sismico Italiano (BSI), 2019). P-wave velocity (km/s)
Depth (km)
Vp/Vs
5.0 6.5 8.05
0 11.1 38
1.73 1.73 1.73
errors ranging from 0.3 to 10 and 0.5 to 14.1 km, for horizontal and vertical coordinates. Only for the earthquakes located with higher accuracy, we computed the local magnitude (ML) on horizontal components of ground displacement by convolving the instrument-corrected seismograms with the instrumental response of the Wood-Anderson seismograph (Table 2), following Hutton and Boore (1987), with errors estimated through the robust Huber mean. The Magnitude ML computed for these 5 earthquakes ranging from 0.42 and 1.29.
4. Seismic imaging method and results For imaging the subsurface structure around the Sotacarbo Fault Lab site, and define the position and geometry of the fault, we performed a passive experiment with ambient noise tomography. The dispersion curve of Rayleigh waves for different inter-station paths, from the correlation properties of ambient noise, are inverted to derive 2D and 3D phase velocity maps at distinct frequencies (Shapiro and Campillo, 2004; Saygin and Kennett, 2012). Throughoutthis section, with ‘noise’ or ‘ambient noise’ we refer to the background Earth’s
3.1. The dense seismic network around the fault field lab After the definition of the seismic baseline in South-Western Sardinia, we focused our attention on the small area (about 13 Km2) of
Fig. 4. Map and vertical sections of the seismicity recorded and located during the experiment. Red triangles represent the stations of the first-large temporary network equipped with a 5-s period sensors. Blue triangles represent stations equipped with a broadband sensor. Black circles represent the 5 earthquakes with highquality location. The green circles represent the other 77 earthquakes with low location accuracy. Black box indicates the Sotacarbo Fault Lab area. 5
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Table 2 Origin time (T0), coordinates and magnitudes (ML) of the 5 high-quality location earthquakes. Year(YYYY)Month(MM)Day(DD) Hour(HH)Minute(MM)Second(SS.SS)
LATITUDE
LONGITUDE
DEPTH (km)
ML
ML Error
2015-03-11 2015-03-26 2015-04-28 2015-05-06 2015-05-06
38N58.47 39N10.38 39N03.81 39N10.06 39N11.67
8E27.01 8E32.41 8E46.38 8E35.47 8E34.18
17.41 2.87 10.31 2.41 2.81
1.07 1.03 1.24 0.42 1.29
+/- 0.27 +/- 0.26 +/- 0.17 +/- 0.26 +/-0.32
16:44:43.18 10:49:28.76 03:34:16.79 11:24:35.81 11:37:07.84
Fig. 5. Sensitivity maps at day- (left) and night-time (right) for the Sotacarbo Fault Lab permanent seismic network. The dashed black box indicates the Sotacarbo Fault Lab area.
vibrations induced by either natural or artificial causes. With ‘recording’ we refer to the measurement of the vertical component of ground velocity. Furthermore, we assume that the noise wavefield is temporally stationary and spatially isotropic. The temporal average of the narrow-band correlation coefficient of the ambient noise recorded by two stations is directly related to the phase velocity of seismic waves at that particular frequency (Aki, 1957; Chávez-Garcı́a et al., 2005). Hence, by repeating correlation estimates over a set of frequencies, one may derive a dispersion curve associated with the propagation path in between the two sites. Before deriving correlation estimates, noise recordings were subjected to a pre-conditioning, followed by an amplitude-normalization procedure aimed at removing the effect of transient, large-amplitude signals of either natural or artificial origin (Bensen et al., 2007). We then calculated the frequency-dependent correlation estimates for all the independent station pairs over the 0.2–1.5 Hz frequency band. For each propagation path, a final correlation function is obtained from the averaging of individual correlation estimates calculated over several, 6hour-long time intervals. Under the assumption that the vertical-component noise wavefield is dominated by fundamental-mode Rayleigh waves, we finally inverted individual correlation functions for the corresponding dispersion curves. For this procedure, we first calculated a trial estimate of the dispersion function using the direct-search method described in Saccorotti et al. (2001). This function is then used as a starting model for an iterative, non-linear least square inversion (Menke and Jin, 2015, Fig. 6).
stations, we obtained the frequency-dependent travel-times for the corresponding ray path. These times were finally inverted in a 2-D tomographic approach for deriving maps of phase velocity at different frequency bands (see Cauchie and Saccorotti, 2013). The tomography problem is solved by parameterizing the propagation medium according to rectangular cells, each of constant velocity. Under the further assumption that linear ray theory holds, the tomographic problem is represented by a set of equations in which individual travel times are expressed as a linear combination of ray length and slowness (the inverse of velocity) in individual cells. This set of linear equation is sparse and, generally, rank-deficient. For its solution, we used a least-squares inversion subjected to the constraints of solution smoothness and minimum distance from a starting, a-priori model (Menke, 1989). At each reference frequency, we used an a-priori homogeneous model whose velocity was obtained by the linear fit of the measured travel times against the corresponding ray lengths. The damping parameter, which controls the weight assigned to the solution smoothness constrain, was set equal to 20 following a trialand-error procedure based on the examination of the trade-off curves (Aster et al., 2005). Examination of the travel-time vs. distance relationships shows that the propagation times depend principally on distance and, as expected, they do increase for increasing station separation. The discrepancies with respect to the average velocity are not systematically related to any particular inter-station azimuth. This implies that the assumption of a spatially-isotropic noise wave-field is satisfied; in case of directional noise sources, in fact, travel time residuals would be markedly correlated with specific inter-station azimuths (Fig. 7). Phase velocity maps were calculated at a set of discrete frequencies spanning the 0.5 Hz–1.5 Hz frequency range, with 0.1 Hz spacing. The
4.1. Phase velocity maps From the phase-velocity dispersion curves obtained at any pair of 6
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 6. Inversion of the correlation function for a dispersion curve at a sample station pair. The panel at the left shows the trial dispersion curve derived from the procedure described in Saccorotti et al. (2001) (gray dashed line), which is then used as a starting model for the iterative, non-linear inversion, whose result is represented by the black line. Black dots are phase velocity estimates obtained from zero-crossings of the correlation function. The panel at the right shows the timeaveraged, experimental correlation function (gray line) and that predicted on the base of the best-fitting dispersion relationship. Black circles mark the zero crossings of the correlation curve.
emerge at the edges of the plain. These latter anomalies appear as isolated spots over the lowermost frequency bands (i.e., the largest depths), but progressively extend throughout both edges of the basin as frequency increases (i.e., at shallower depths). This observation supports the presence of N–S trending normal fault, which displaced the deepest part of the basin. 4.2. Profile along the linear array The ambient noise correlation analysis is then extended to stations from the linear profile shown in Fig. 10a. In this case, we grouped the correlation coefficients obtained over the same frequency band but for different inter-station distances, then conducting a direct search for that phase velocity which minimizes the difference between the observed and predicted correlation curves. For the direct search, we considered 100 phase velocity values spanning the 0.2 km/s - 3 km/s range. We conducted the analysis over groups formed by 4 consecutive stations (distances spanning the 120 and 360 m), and using the correlation coefficients calculated for all the independent inter-station distances associated with that group. The procedure was then iterated considering the next group, which was derived by translating the index of the first station by one unit. We thus obtained a set of 8 dispersion curves, which were combined to obtain a map where phase velocity is expressed as a function of frequency and offset along the profile. The resulting phase velocity section (Fig. 10b) shows two main discontinuities, positioned at about 400−500 m and 600−700 m along the profile. We associate the shallowest low and high velocities (Vr = Rayleigh phase velocity) to the Quaternary alluvial deposits (dark blue, VR ∼ 0.600 km/s) and pyroclastic products and ignimbrites of the lower and middle Miocene (yellow-green, VR ∼0.9 km/s), respectively. The underlying, high-velocity level highlighted in the NW margin of the profile (1< VR <1.2 km/s) could instead correspond to the magmatic products of the Lower Miocene that constitute the basement of the basin. For the same reasons outlined above, we didn’t perform any further inversion for transforming the phase velocity section of Fig. 10b into a corresponding depth section of the shear-wave velocity along the
Fig. 7. Travel times as a function of inter-station distance at the sample frequency F = 0.5 Hz. An average velocity (dashed line) is obtained from the slope of a straight line fitted to experimental data in a least-square sense. Colors refer to the orientation of individual station pairs, according to the color bar at the right.
model space was defined by a set of 8 × 5 rectangular cells with sides 650 m and 660 m along the EW and NS directions, to optimize image fidelity and model resolution. The ray crossing is adequate for ensuring a good resolution of velocities in the central cells, and still acceptable at peripheral cells with a lower hit counts, according to synthetic tests (see Fig. 8b, c, d). Fig. 9 shows the phase-velocity maps obtained from the tomographic inversions at the different target frequencies. Conversion of the phase velocity maps to a 3D image of shear-wave velocities would require an inversion of the dispersion functions associated with individual nodes of the maps. Such non-unique solutions strongly depend on inversion regularizations and starting models. Since our poor a priori knowledge of the shear-wave velocities at depth, we prefer to limit our analysis to the phase velocity maps. The Matzaccara Quaternary basin is identified as a low-velocity area bounded to the west and to the east by high velocity anomalies imputable to the lower and middle Miocene volcanic products, which 7
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 8. (a) Geometry of station pairs and corresponding seismic rays used in the tomographic inversion for phase velocity maps, (b) Number of rays crossing individual cells, (c) Input model for the checkerboard validation test. The model is defined by the alternance of cells with exhibiting +/- 10 % velocity anomalies with respect to the average velocity retrieved at that particular frequency. The synthetic travel-times were perturbed by adding random errors extracted from a normal distribution having zero mean and standard deviation equal to the standard deviation of the residual distribution between observed data and travel-times predicted by the starting uniform velocity model. (d) Image recovered after the tomographic inversion of synthetic travel-times calculated for the model shown in (c). For the central part of the network, the pattern of alternating positive and negative anomalies is well recovered. The color bars indicates the phase velocity in km/s.
profile. However, following the method proposed by Socco et al. (2017), one may obtain an approximate relationship between the observed wavelengths and the resolvable depths, thus getting a qualitative indication about the depth of the main velocity discontinuities. In this manner, we estimate that the fault (or fault segments) dissecting the basin would affect a crustal portion spanning the 50−450 m depth interval. The frequency band characterized by the lowermost velocities (dark blue in Fig. 10b) becomes progressively wider as one moves toward the SE margin of the profile, indicating a progressive thickening of the sedimentary / volcanoclastic deposits in-filling the basin. The maximum thickness of these deposits is inferred to be at least on the order of 150−200 m, and this interval may be taken as a rough estimate of the minimum displacement occurred along the faults bordering the basin. However, it cannot be excluded that the sedimentary cover gets even thicker SE of the imaged section, i.e. toward the center of the basin. In that case, the above estimate on the displacement along the fault would represent a lower bound on the actual cumulated slip.
4.3. The Matzaccara buried fault The good correspondence between the phase velocity maps of Fig. 9 and the phase velocity section along the linear profile (Fig. 10b) allows identification of an east-dipping fault that dislocates the shallow layers, and contributes to the formation of the Matzaccara Quaternary basin. The low-velocity anomaly coincident with the basin elongation strikes about N-NE, indicating that even the east-dipping fault imaged in the 2D section, is elongated along the same azimuth. This latter fault is composed by a superficial splay positioned between the surface and 150−200 m deep, located around 400−500 m away from station AT11 (Fig. 10b), which connects to a deeper segment spanning the 200−400 m depth interval, positioned at a distance of 600−700 m from station AT11. The fault displaces the volcanic formations for the depth range in which it was detected, but it is not possible to verify with certainty any more recent activity related to Quaternary deposits. 5. Discussion and conclusions The development and acceptance of CCUS requires careful analyses 8
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 9. Sample phase velocity maps at different reference frequencies. Color bar indicates phase velocity in km/s. Not resolved areas (hit counts <5) are masked.
of the hazards associated with the injection and storage activity, passing through the implementation of appropriate mitigation strategies and best practices of the risk management (White and Foxall, 2016). The caprock integrity and the modification of the stress field due to the confinement of large gas volumes are primary aspects for hazard assessment (Verdon et al., 2011; Verdon, 2014; Cantucci et al., 2016). Experiments at on-shore fields laboratories (such as the Sotacarbo Fault Lab, for instance) and pilot sites are therefore the most valuable strategy to understand how the deep systems evolve in response to fluid injection, learning by doing.
Under this perspective, the level of natural seismicity before the beginning of the activities is a pre-requisite for monitoring and identifying all the possible alterations in the subsurface during the fluid injection (Cappa and Rutqvist, 2011; Stork Anna et al., 2015, 2018). Seismological applications in CCUS pilot sites focus on the evaluation of the seismicity background before injection activities (Chen et al., 2015), the definition of the target structure and the identification of the faults interest it (Alcalde et al., 2013) and could be a useful tool in the study of geomechanical response of a site to the CO2 injection and to detect the seismic footprint of the CO2 leakage (Stork Anna et al., 2018). Past 9
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 10. a) The seismic transect across the Sotacarbo Fault Lab used for deriving the 2-D phase velocity image shown in 10 b. Blue circles mark the stations, while the red boxes 1 and 2 represent the surface projection of the two branches (1 and 2) of the fault highlighted in the phase-velocity vertical section. The figure shows just a small portion of the Matzaccara basin. b) Rayleigh-wave phase velocities derived from correlation analysis. Color bar indicates phase velocity in km/s. Velocities are mapped in their dependence upon frequency and distance along the seismic profile. Note that frequency ranges of 2.4-1.6 Hz and 1.6-0.8 Hz correspond to depth intervals of 0–200 and 200−400 m, respectively. Black lines mark two evident lateral discontinuities likely attributable to the Matzaccara normal fault, with a geometry similar to the sketch in the inset. The number 1 and 2 are referred to the dashed boxes in Fig. 10a.
experiences evidenced the strict requirement of having correct baselines to interpret the microseismicity revealed during the exploitation activities (Stork Anna et al., 2015, 2018). The final aim is to use monitoring data in a proactive real time forecast system, adapted to a defined site-specific sensibility. One second key aspect for the safe storage of large CO2 volumes is the characterization of fault and fracture systems in the subsurface. Especially important is the individuation of faults quiescent since long time and capable of producing damaging seismic events. Examples of faults activation during fluid injection became evident in different energy methods and context (waste water disposal, Keranen et al., 2014; E.G.S., Kraft and Deichmann, 2014; E.G.S., Diehl et al., 2017; waste water disposal, E.O.R., hydrofracturing, C.C.S. E.G.S., McGarr et al.,
2015; C.C.S., Stork Anna et al., 2018). In several cases, the presence of these faults was unknown before the starting of the industrial injection activities, although they dissected the basement under the sedimentary cover. Seismological applications are designed to define the presence of faults at depth, but understanding if those faults are critically stressed and prone to fail is very challenging. Today, we still lack methods that clearly address the identification of critically-stressed faults among those identified by seismic imaging of the upper crust (Ellsworth, 2013; Buttinelli et al., 2016). Faults and fracture systems are also preferential paths for the CO2 plume (Chang and Bryant, 2009) and the definition of the correct asset of fractures in the subsurface is strategic to tune the location of injection and monitoring wells, and to identify possible gas escape pathways. 10
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
We have investigated these issues in the Sulcis basin (South-Western Sardinia), where the Sotacarbo Fault Lab will be installed. Several years of seismological surveys helped to define the seismicity background of the area. This task has been particularly difficult because the area is slowly deforming (Palano, 2015) and the seismic rate is low, requiring longer observation times for the definition of seismic baselines and identification of faults at depth. We performed a three-step baseline survey with a progressive focusing and densification around the target, starting from a regional scale. In fact, we ignored at the beginning the seismic level for a broad region of Sardinia, that was never the subject of a focused monitoring before. From the first, large-aperture network, we were able to identify a few earthquakes in the area with magnitude around ML -1.0 (Fig. 3), that constitutes an unprecedented result for this low seismic area. Earthquakes hypocenters are sparse and do not align on specific faults or structures. Hypocentral accuracy do not permit detailed investigations. Then, we focused on the target site with a composite network formed by a permanent backbone and extremely dense but shortduration temporary arrays. In this latter case, we did not detect any seismicity down to magnitude as low as the that indicated by the sensitivity study (ML -1 / -1.5, within the top 5 km depth). The other strategic part of the experiment was dedicated to the reconnaissance and characterization of faults in the subsurface. Given the absence of seismicity, passive seismic methods like travel-time tomography were not applicable. We therefore produced a detailed imaging of the subsurface using ambient noise tomography (Lin et al., 2013). We have thus been able to identify an east-dipping normal fault in the 400 m below the surface (Fig. 10b). Although the lack of resolution over the shallowest tens of meters prevents from estimating the last episode of slip on the fault, its last activity likely occurred during Quaternary, according with the displaced shallow sediments. The fault trends N-NE, and it possibly extends for at least 2−3 km or more, outside of the resolution provided by the tomographic model (Fig. 9). Throughout the two-year-long seismic monitoring of the area, no earthquakes have been detected in association to this fault, down to the detection threshold of ML -1.0.
References Aki, K., 1957. Space and time spectra of stationary waves with special reference to microtremors. Bull. Earthquake Res. Inst. Univ. Tokyo 35, 415–456. Alcalde, J., Martí, D., Calahorrano, A., Marzan, I., Ayarza, P., Carbonell, R., Juhlin, C., Pérez-Estaún, A., 2013. Active seismic characterization experiments of the Hontomín research facility for geological storage of CO2, Spain. Int. J. Greenh. Gas Control. 19, 785–795. https://doi.org/10.1016/j.ijggc.2013.01.039. ISSN 1750-5836. Arts, R., Eiken, O., Chadwick, A., Zweigel, P., van der Meer, L., Zinszner, B., 2004. Monitoring of CO2 injected at Sleipner using time-lapse seismic data. Energy 29 (9–10), 1383–1392. https://doi.org/10.1016/j.energy.2004.03.072. ISSN 03605442. Aster, R., Beaudoin, B., Hole, J., Fouch, M., Fowler, J., James, D., 2005. IRIS Seismology Program Marks 20 Years of Discovery. Eos, Transactions American Geophysical Union. Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M., Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad‐band surface wave dispersion measurements. Geophys. J. Int. 169, 1239–1260. https://doi.org/10.1111/j.1365-246X.2007.03374.x. Boatwright, John, 1978. Detailed spectral analysis of two small New York state earthquakes. Bull. Seismol. Soc. Am. 68 (4), 1117–1131. Bollettino Sismico Italiano (BSI) - Gruppo di Lavoro Sequenza Centro Italia (2019). Rapporto Bollettino Sismico Italiano sulla revisione della sequenza sismica del centro Italia 24 agosto 2016 - 31 agosto 2018. Bommer, J.J., Crowley, H., Pinho, R., 2015. A risk-mitigation approach to the management of induced seismicity. J. Seismol. 19 (2), 623–646. Boot-Handford, M.E., Abanades, J.C., Anthony, E.J., Blunt, M.J., Brandani, S., Mac Dowell, N., Fernandez, J.R., Ferrari, M.C., Gross, R., Hallett, J.P., Haszeldine, R.S., Heptonstall, P., Lyngfelt, A., Makuch, Z., Mangano, E., Porter, R.T.J., Pourkashanian, M., Rochelle, G.T., Shah, N., Yao, J.G., Fennell, P.S., 2014. Carbon capture and storage update. Energy Environ. Sci. 7, 130–189. Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes. J. Geophys. Res. 75 (26), 4997–5009. https://doi.org/10.1029/ JB075i026p04997. Buttinelli, M., Improta, L., Bagh, S., Chiarabba, C., 2016. Inversion of inherited thrusts by wastewater injection induced seismicity at the Val d’Agri oilfield (Italy). Sci. Rep. https://doi.org/10.1038/srep37165. Cantucci B., Buttinelli M, Procesi M., Sciarra A., Anselmi M. (2016) Algorithms for CO2 storage capacity estimation: review and case study. Book Chapter into "Geologic Carbon Sequestration: Understanding Reservoir Concepts", Section 2 “Geological Characterisation of storage sites”. published by Springer under the joint editorial handling of Prof. T. N. Singh (Group Lead, Rock Mechanics, IIT Bombay) DOI https:// doi.org/10.1007/978-3-319-27019-7. Cappa, F., Rutqvist, J., 2011. Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO2. Int. J. Greenh. Gas Control. 5 (2), 336–346. Cauchie, L., Saccorotti, G., 2013. J. Seismol. 17, 335. https://doi.org/10.1007/s10950012-9323-6. Chang, Kyung Won, Bryant, Steven, 2009. The effect of faults on dynamics of CO2 plumes. Energy Procedia 1, 1839–1846. https://doi.org/10.1016/j.egypro.2009.01. 240. Chávez-Garcı́a, F.J., Rodrı́guez, M., Stephenson, W.R., 2005. An alternative approach to the SPAC analysis of microtremors: exploiting stationarity of noise. Bull. Seism. Soc. Am. 95, 277–293. https://doi.org/10.1785/0120030179. Chen, T., Lianjie, Huang, 2015. Seismicity characterization around the Farnsworth field site for combined large-scale CO2 storage and EOR. Int. J. Greenh. Gas Control. 33, 63–68. https://doi.org/10.1016/j.ijggc.2014.11.024. ISSN 1750-5836. Diehl, T., Kraft, T., Kissling, E., Wiemer, S., 2017. The induced earthquake sequence related to the St. Gallen deep geothermal project (Switzerland): fault reactivation and fluid interactions imaged by microseismicity. J. Geophys. Res. Solid Earth 122, 7272–7290. https://doi.org/10.1002/2017JB014473. Dipietro, J., et al. “The Role of Naturally-Occurring CO2 Deposits in the Emergence of CO2 Enhanced Oil Recovery” http://co2conference.net/pdf/1.2-Slides_ DiPietroCO2Sources2011CO2FloodingConf.pdf. Ellsworth, W.L., 2013. Injection-induced earthquakes. Science 341 (6142), 1225942. https://doi.org/10.1126/science.1225942. Förster, Andrea, Norden, Ben, Zinck-Jørgensen, Kim, Frykman, Peter, Kulenkampff, Johannes, Spangenberg, Erik, Erzinger, J.örg, Zimmer, Martin, Kopp, J.ürgen, Borm, G.ünter, Juhlin, Chris, Cosma, Calin-Gabriel, 2006. Suzanne Hurter; Baseline characterization of the CO2SINK geological storage site at Ketzin. Germany. Environmental Geosciences 13 (3), 145–161. https://doi.org/10.1306/eg. 02080605016. Grigoli, F., Cesca, S., Priolo, E., Rinaldi, A.P., Clinton, J.F., Stabile, T.A., Dost, B., Fernandez, M.G., Wiemer, S., Dahm, T., 2017. Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: a European perspective. Rev. Geophys. 55, 310340. https://doi. org/10.1002/2016RG000542. Gruppo di Lavoro, 2004. Redazione della mappa di pericolosità sismica prevista dall’Ordinanza PCM 3274 del 20 marzo 2003. Rapporto conclusivo per il Dipartimento della Protezione Civile. INGV, Milano-Roma aprile 2004, 65 pp + 5 enclosures. Hutton, Boore, 1987. The ML scale in southern California. Bull. Seismol. Soc. Am. 77 (6), 2074–2094. Kanamori, H., 1977. The energy release in great earthquakes. J. Geophys. Res. 82 (20), 2981–2987. https://doi.org/10.1029/JB082i020p02981.
CRediT authorship contribution statement M. Anselmi: Conceptualization, Methodology, Data curation, Writing - original draft. G. Saccorotti: Data curation, Methodology, Writing - original draft. D. Piccinini: Methodology, Data curation. C. Giunchi: Methodology, Data curation. M. Paratore: Data curation, Software. P. De Gori: Software, Methodology, Data curation. M. Buttinelli: Conceptualization, Validation. E. Maggio: Resources, Funding acquisition. A. Plaisant: Resources, Supervision. C. Chiarabba: Conceptualization, Writing - original draft, Supervision. Declaration of Competing Interest None. Acknowledgements This study has been carried out within the “Research on Electric System” project funded by the Italian Ministry of Economic Development. The Sotacarbo Fault Lab is under construction within the “Centre of Excellence on Clean Energy” project (CUP: D83C17000370002) funded by the Regional Government of Sardinia (FSC 2014–2020). The data from the station SU26, belonging to the Sotacarbo Fault Lab seismic network, is currently archived in EIDA waveform repository (eida.rm.ingv.it). We are grateful to Gianfranco Colasanti, to the “Rete Sismica Mobile” group at INGV and to the "Corpo Forestale della Sardegna" for their helpful technical and logistics support. We thank Alberto Pettinau of Sotacarbo S.p.A. for his helpful contribution in the review of the paper. 11
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al. Keranen, K.M., Weingarten, M., Abers, G.A., Bekins, B.A., Ge, S., 2014. Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection. Science 345 (6195), 448–451. Kraft, T., Deichmann, Nicholas, 2014. High-precision relocation and focal mechanism of the injection-induced seismicity at the Basel EGS. Geothermics 52, 59–73. https:// doi.org/10.1016/j.geothermics.2014.05.014. ISSN 0375-6505. Lahr, J.C., 1999, revised 2012, HYPOELLIPSE: a computer program for determining local earthquake hypocentral parameters, magnitude, and first-motion pattern: U.S. Geological Survey Open-File Report 99–23, version 1.1, 119 p. and software, available at https://pubs.usgs.gov/of/1999/ofr-99-0023/. Lin, Fan-Chi, Li, Dunzhu, Clayton, Robert W., Hollis, Dan, 2013. High-resolution 3D shallow crustal structure in Long Beach, California: application of ambient noise tomography on a dense seismic array. Geophysics 78 (4), Q45–Q56. https://doi.org/ 10.1190/geo2012-0453.1. Liu, H.J., Were, P., Li, Q., Gou, Y., Hou, Z., 2017. Worldwide Status of CCUS Technologies and Their Development and Challenges in China Geofluids. Article ID 6126505, 25 pages. https://doi.org/10.1155/2017/6126505. McGarr, A., Bekins, B., Burkardt, N., Dewey, J., Earle, P., Ellsworth, W., Ge, S., Hickman, S., Holland, A., Majer, E., Rubinstein, J., Sheehan, A., 2015. Coping with earthquakes induced by fluid injection. Science 347, 830–831. https://doi.org/10.1126/science. aaa0494. Menke, W., 1989. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press, Inc., New York, pp. 289pp. Menke, W., Jin, G., 2015. Waveform fitting of cross spectra to determine phase velocity using Aki’s formula. Bull. Seism. Soc. Am. 105, 1619–1627. https://doi.org/10.1785/ 0120140245. Munafò, I., Malagnini, L., Chiaraluce, L., 2016. On the relationship between mw and ML for small earthquakes. Bull. Seismol. Soc. Am. 106 (5), 2402–2408. https://doi.org/ 10.1785/0120160130. Nicol, A., Carne, R., Gerstenberger, M., Christophersen, A., 2011. Induced seismicity and its implications for CO2 storage risk. Energy Procedia 4, 3699–3706. Palano, M., 2015. On the present-day crustal stress, strain-rate fields and mantle anisotropy pattern of Italy. Geophys. J. Int. 200 (2), 969–985. https://doi.org/10.1093/ gji/ggu451. Ringrose, P.S., Mathieson, A.S., Wright, I.W., Selama, F., Hansen, O., Bissell, R., Saoula, N., Midgley, J., 2013. The in Salah CO2 storage project: lessons learned and knowledge transfer. Energy Procedia 37, 6226–6236. Rovida, A., Locati, M., Camassi, R., Lolli, B., Gasperini, P. (Eds.), 2016. CPTI15, the 2015 Version of the Parametric Catalogue of Italian Earthquakes. Istituto Nazionale di Geofisica e Vulcanologia. Saccorotti, G., Almendros, J., Carmona, E., Ibáñez, J.M., Del Pezzo, E., 2001. Slowness anomalies from two dense seismic arrays at Deception Island volcano, Antarctica. Bull. Seismol. Soc. Am. https://doi.org/10.1785/0120000073. Saygin, E., Kennett, B.L.N., 2012. Crustal structure of Australia from ambient seismic noise tomography. J. Geophys. Res. 117, B01304. https://doi.org/10.1029/ 2011JB008403. Schorlemmer, D., Mele, F., Marzocchi, W., 2010. A completeness analysis of the National Seismic Network of Italy. J. Geophys. Res. 115, B04308. https://doi.org/10.1029/ 2008JB006097. Shapiro, N.M., Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophys. Res. Lett. 31, L07614. https://doi. org/10.1029/2004GL019491.
Socco, L.V., Comina, C., Anjom, F.K., 2017. Time-average velocity estimation through surface-wave analysis: Part 1 — S-wave velocity. Geophysics 82, U49–U59. https:// doi.org/10.1190/GEO2016-0367.1. Stork Anna, L., Verdon, James P., Kendall, J.-Michael, 2015. The microseismic response at the in-salah Carbon Capture and Storage (CCS) site. Int. J. Greenh. Gas Control. 32, 159–171. https://doi.org/10.1016/j.ijggc.2014.11.014. ISSN 1750-5836. Stork Anna, L., Allmark, Claire, Curtis, Andrew, Kendall, J.-Michael, White, Don J., 2018. Assessing the potential to use repeated ambient noise seismic tomography to detect CO2 leaks: application to the Aquistore storage site. Int. J. Greenh. Gas Control. 71, 20–35. https://doi.org/10.1016/j.ijggc.2018.02.007. ISSN 1750-5836. Stucchi, Massimiliano, Meletti, Carlo, Montaldo, Valentina, Crowley, Helen, Calvi, Gian Michele, Boschi, Enzo, 2011. Seismic hazard assessment (2003–2009) for the italian building code. Bull. Seismol. Soc. Am. 101 (4), 1885–1911. https://doi.org/10.1785/ 0120100130. Takagishi, Makiko, Hashimoto, Tsutomu, Horikawa, Shigeo, KinichiroKusunose, Ziqiu Xue, Susan, D.Hovorka, 2014. Microseismic Monitoring at the Large-scale CO2 Injection Site. Cranfield, MS, U.S.A. https://doi.org/10.1016/j.egypro.2014.11.476. Ugalde, A., Villaseñor, A., Gaite, B., Casquero, S., Martí, D., Calahorrano, A., Marzán, I., Carbonell, R., Estaún, A.P., 2013. Passive seismic monitoring of an experimental CO2 geological storage site in Hontomín (Northern Spain). Seismol. Res. Lett. 84 (1), 75–84. Van Thienen-Visser, K., Breunese, J.N., 2015. Induced seismicity of the Groningen gas field: history and recent developments. Lead. Edge 34https://doi.org/10.1190/ tle34060664.1. 664–666, 668–668, 670–671. Verdon, J.P., 2014. Significance for secure CO2 storage of earthquakes induced by fluid injection. Environ. Res. Lett. 9 (6). Verdon, J.P., Kendall, J.-M., White, D.J., Angus, D.A., 2011. Linking Microseismic Event Observations with Geomechanical Models to Minimise the Risks of Storing CO2 in Geological Formations. https://doi.org/10.1016/j.epsl.2011.02.048. Vilarrasa, V., Carrera, J., 2015. Geologic carbon storage is unlikely to trigger large earthquakes and reactivate faults through which CO2 could leak. Proc. Natl. Acad. Sci. 112, 5938–5943. Vilarrasa, Victor, Carrera, Jesus, Olivella, Sebastià, Rutqvist, Jonny, Laloui, Lyesse, 2019. Induced seismicity in geologic carbon storage. Solid Earth 10 (871–892). https://doi. org/10.5194/se-10-871-2019. White, J.A., Foxall, W., 2016. Assessing induced seismicity risk at CO2 storage projects: recent progress and remaining challenges. Int. J. Greenh. Gas Control. 49, 413–424. Woessner, Jochen, Laurentiu, Danciu, Giardini, Domenico, Crowley, Helen, Cotton, Fabrice, Grünthal, Gottfried, Valensise, Gianluca, Arvidsson, Ronald, Basili, Roberto, Demircioglu, Mine, Hiemer, Stefan, Meletti, Carlo, Musson, RogerW., Rovida, Andrea, Sesetyan, Karin, Stucchi, Massimiliano, 2015. The 2013 European Seismic Hazard Model: key components and results. Bull. Earthq. Eng. 13, 1–44. https://doi. org/10.1007/s10518-015-9795-1. Zang, A., Oye, Volker, Jousset, Philippe, Deichmann, Nicholas, Gritto, Roland, McGarr, Art, Majer, Ernest, Bruhn, David, 2014. Analysis of induced seismicity in geothermal reservoirs -An overview. Geothermics 52, 6–21. https://doi.org/10.1016/j. geothermics.2014.06.005. ISSN 0375-6505. Zheng, S., Sun, X., Song, X., Yang, Y., Ritzwoller, M.H., 2008. Surface wave tomography of China from ambient seismic noise correlation. Geochem. Geophys. Geosyst. 9, Q05020. https://doi.org/10.1029/2008GC001981. Zoback, M.D., Gorelick, S.M., 2012. Earthquake triggering and large-scale geologic storage of carbon dioxide. Proc. Natl. Acad. Sci. 109, 10164–10168.
12
Contents lists available at ScienceDirect
International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc
Microseismic assessment and fault characterization at the Sulcis (SouthWestern Sardinia) field laboratory
T
M. Anselmia,*, G. Saccorottib,1, D. Piccininib,1, C. Giunchib,1, M. Paratoreb,1, P. De Goria, M. Buttinellia, E. Maggioc,2, A. Plaisantc,2, C. Chiarabbaa a
Istituto Nazionale di Geofisica e Vulcanologia, Rome, Italy Istituto Nazionale di Geofisica e Vulcanologia, Pisa, Italy c Sotacarbo S.p.A., Carbonia (SU), Italy b
ARTICLE INFO
ABSTRACT
Keywords: CCUS Fault Lab Seismic baseline Ambient noise tomography Sardinia Sulcis basin
The general acceptance of the CO2 geological storage by stakeholders passes through the assessment and mitigation of risks, potentially induced or increased by the disposal activity. Injection of moderate to large quantities of CO2 in the sub-surface may unbalance local stress and trigger earthquakes if faults are critically stressed, condition that is not easily verifiable. Pilot sites are therefore the best way to proceed further in order to address such challenging issues. In such cases, the reconnaissance of faults and seismicity in the sub-surface, before the onset of activity, is mandatory. In this paper, we present studies carried out in the site where the Sotacarbo Fault Lab is going to be installed. This facility will be located in a very low seismic hazard region of central Mediterranean, where reports on historical large earthquakes are poor. We show results from a series of experiments aimed to monitor the background seismicity around the pilot site. As expected, seismicity is almost absent down to small magnitude close to the future injection-test well. Further seismic imaging of the subsurface layers obtained by ambient noise tomography offers the ability to resolve the presence of a seismicityfree fault located in the first 200 m below the surface, of which the last episode of activity is difficult to assess. Our results encourage the use of this site to follow the response of the system to injection of small quantity of CO2.
1. Introduction The crucial need to reduce the CO2 emission requires a multiplicity of efforts, among which carbon capture, utilization and storage (CCUS) is a viable opportunity (Dipietro et al., 2011; Boot-Handford et al., 2014; Liu et al., 2017). Geological storage of CO2 might produce environmental and safety hazards, (e.g., accidental leakage during the injection at critical pressure, slow leakage after injection in the reservoir following diffusion and upwelling processes in presence of cracking/faulting of the cap rocks) that need to be reliably assessed (Nicol et al., 2011; Vilarrasa et al., 2019). Geo-energy activities that involve fluid injection/extraction may change the state of stress in the subsurface causing deformation, microcracking and fault reactivation (Grigoli et al., 2017). Induced seismicity has received increasing attention by stakeholders (Ellsworth, 2013) and is potentially one of the main obstacles to the volumetric upscale of geo-energy projects.
This is particularly true for areas where the natural seismic hazard is low like in the Sulcis area (Gruppo di Lavoro, 2004; Woessner et al., 2015) and therefore the existing infrastructure may not be adequate to withstand the impact of potentially large anthropogenic seismicity (Takagishi et al., 2014; van Thienen-Visser and Breunese, 2015). The capacity of storing large quantities of CO2 scales with the volume that will be affected by changes in pressure and stress, increasing the risk of intercepting active faults and inducing large events (Keranen et al., 2014). Different views have emerged about up-scaling storage projects and their applicability or convenience. In particular, some opinions emerged proposing that large-scale CCUS is risky and possibly unsuccessful (Zoback and Gorelick, 2012), for the high probability that earthquakes could violate the seal integrity. Conversely, other opinions assessed that CO2 storage can be performed safely, after proper site characterization and with adequate pressure management (Vilarrasa and Carrera, 2015; Vilarrasa et al., 2019). At present, pilot and
Corresponding author. E-mail address: [email protected] (M. Anselmi). 1 Saccorotti, Piccinini, Giunchi, Paratore belong to the branch of Pisa Istituto Nazionale di Geofisica e Vulcanologia, Via Cesare Battisti, 43-56125 Pisa-Italy. 2 E. Maggio and A. Plaisant belong to the organization Sotacarbo S.p.A - Grande Miniera di Serbariu, 09013 Carbonia (SU) - Italy. ⁎
https://doi.org/10.1016/j.ijggc.2020.102974 Received 6 August 2019; Received in revised form 15 January 2020; Accepted 17 January 2020 Available online 25 January 2020 1750-5836/ © 2020 Elsevier Ltd. All rights reserved.
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
demonstration projects focus on verifying the storing capacity and the development of monitoring best practices. Considering the acceptance of a reasonable level of seismic risk, the general goal is to chart a course and defines clear and plain protocols for safe and efficient monitoring that afford reliable and accurate hazard forecasts in a real-time management of the risk (Bommer et al., 2015). One of the main challenges is to define the existence of faults capable of producing damaging earthquakes. Since activities are generally shallow, even small to moderate magnitude earthquakes can produce significant ground shaking at the surface. In addition, faults that could fail with small to moderate earthquakes are more difficult to be detected in the subsurface. Earthquakes induced by fluid injection often develop on faults barely visible in seismic imaging of the subsurface (Buttinelli et al., 2016). The reliable identification and mapping of aseismic, but possibly critically stressed faults is therefore a first key step in hazard assessment, as became clear in unconventional geothermal systems, like the Enhanced Geothermal Systems (EGS, see Zang et al., 2014). Moreover, faults can act as enhancer or obstacle to fluid flow and their reconnaissance and characterization are important to plan efficient CO2 sequestration. In this paper, we investigate such aspects for the Sotacarbo Fault Lab (South-Western Sardinia, Italy), a facility specifically designed as a test bed infrastructure in order to develop and test in-field conditions instruments and methods to be applied in the commercial-scale CO2 geological storage site elsewhere. The facility will be located in SouthWestern Sardinia, a region with low natural seismic hazard (Stucchi et al., 2011; Rovida et al., 2016). According to that, the national permanent seismic monitoring system has an intrinsic high magnitude threshold for detection of local earthquakes in Sardinia, where the Magnitude of completeness (Mc) is higher than 2.9 (Schorlemmer et al., 2010). To improve earthquake detection, various dense seismic arrays have been deployed over the past 4 years. The Sulcis temporary seismic network was installed for a period of about 15 months since July 2014 in order to monitor the South-Western Sardinia (Fig. 1a). After that, in July 2016 we installed the Sotacarbo Fault Lab permanent seismic network of 5 stations around the fault lab, complemented with a temporary network consisting of additional 11 stations operative from the first days of June 2017 to end of August 2017 (Fig. 1b). We present results of the multiscale survey and quality assessment of the overall monitoring efficiency. In order to enhance the imaging of the fault in the lab site, we performed a focused experiment for ambient noise tomography, a very low cost and efficient tool for subsurface imaging especially in case of absence of local seismicity (Zheng et al., 2008). During the first days of October 2017, we installed a temporary linear seismic array for about 4 days (see line in Fig. 1b), while the Sotacarbo Fault Lab permanent seismic network was operating. Data from both the passive seismic linear array and the Sotacarbo permanent and temporary seismic networks enable the tomographic reconstruction of the fault.
monitoring system that will operate during the future activities of the Sotacarbo Fault Lab. The South-Western Sardinia exhibits a very low seismicity level; the historical catalog reports only two earthquakes occurred in 1616 and 1771, with estimated moment magnitudes of 4.9 and 4.4, respectively (Rovida et al., 2016; Fig. 1). The low seismicity rate of the area is also confirmed by the instrumental catalog issued by Istituto Nazionale di Geofisica e Vulcanologia (INGV). For the time interval 1985-present, it reports only five earthquakes, among which the ML 4.1 occurred in 2006, a few kilometers offshore the Gulf of Cagliari (Fig. 1a). In order to increase the current seismic sensitivity of the national network toward the detection of low-magnitude earthquakes, we installed the first large temporary network of seismic stations in Sulcis, composed of high-dynamic, local-recording digitizers Reftek 130-1. The seismic network was composed of 8 Lennartz Le3D/5 s sensors and 2 Nanometrics Trillium Compact 120 s (see Fig. 1a). The network operated in continuous – recording mode between July 2014 and September 2015. The main requisite for seismic baseline is to capture the smallest magnitude occurring in the area. We therefore started analyzing the level of possible detection though the evaluation of the noise condition at each individual site, by using the Probability Density Function (PDF) of the power spectra evaluated for the three components of ground velocity. In addition, local amplification phenomena were investigated using the Horizontal-to-Vertical spectral ratios (Fig. 2). To estimate the detection capability of the temporary network, we implemented from scratch a procedure based on the comparison of the observed noise spectra with those predicted for a point source of known location and magnitude. The noise spectral amplitude was obtained by computing average spectral density for 1 week of seismic signal between 1 and 2 a.m. for the nighttime measures, and between 1 and 2 p.m. for daytime measures. First, we parameterize the Earth’s volume below the network using a regular grid of nodes spaced by 1 km along the NS and EW directions, and 5 km along the vertical direction. Each grid node represents a virtual seismic point source defined using a Brune’s model [Brune, 1970] at which we assign a moment magnitude Mw. Then we derive the scalar seismic moment and the expected displacement spectrum (following Kanamori, 1977 and Boatwright, 1978) and then we propagate the seismic amplitude in a 5.8 km/s half-space velocity model with homogeneous frequency independent attenuation (Q = 250). For a given source location and magnitude, the ground motion spectrum predicted at any individual station (taking into account the topographic differences between the stations) is then compared to the average noise spectrum observed at that particular site. If the ratio between the predicted and noise spectra is larger than 2 over the 3−25 Hz frequency band, than the procedure assumes that that station is able to detect an earthquake of that particular magnitude. When the number of stations detecting the earthquake is greater than 4, that magnitude value is accepted and attributed to the grid node under evaluation. The procedure is iterated over the entire grid of virtual sources, for moment magnitudes Mw spanning the 0.1–2.2 range, then converted in local magnitudes ML according to Munafò et al. (2016). This allows deriving maps of the minimum detectable magnitudes as a function of source location at the hypocentral depth of 5 km b.s.l.. We took into account two different scenarios, at day and night time. The estimated theoretical detection level within the Sotacarbo Fault Lab is around the ML -1.0 (see Fig. 3). This is a considerable improvement, considering the initial high detection threshold in the area ensured by the INGV national network.
2. Baseline monitoring of the Sulcis site The Sulcis is a mining area of South-Western Sardinia, where a former coalfield within the Matzaccara Quaternary basin has been selected as the site of the Sotacarbo Fault Lab (Fig. 1a and b). Aim of this facility is to study CO2 leakage processes through faults and develop advanced monitoring tools to be applied in commercial-scale CO2 geological storage applications worldwide. The characterization of the site in terms of adequateness with respect to natural hazards is therefore a pre-requisite. Following the general framework of the technical and geological feasibility evaluation of CCUS (Arts et al., 2004; Förster et al., 2006; Ringrose et al., 2013), the definition of a seismic baseline is an important aspect (Ugalde et al., 2013; Chen et al., 2015). Seismic baselines are crucial to define the reference level of natural seismicity before the onset of activities and provide hints about the presence of active faults. Data acquired during the baseline aid to design the seismic
3. Results of seismic monitoring and earthquakes location From inspection of one year of continuous recordings from the firstlarge temporary network (red and blues triangles in Fig. 1a), we 2
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 1. a) Map of the study area in the South-Western Sardinia. Red triangles represent the stations of the first-large temporary network. The blue triangles represent the broadband (BB) stations of the same network. Dashed black square indicates the location of the Sotacarbo Fault Lab (see Fig. 1b). Historical earthquakes for South-Western Sardinia are marked by purple stars. Instrumental earthquakes are marked by brown circles. b) Map of the Sotacarbo Fault Lab area. The cyan and red triangles represent the stations of the permanent network used to monitoring the Sotacarbo Fault Lab area. Green triangles represent the Sotacarbo temporary seismic network added for the noise tomographic experiment. The blue circles represent the passive linear seismic array deployed across the presumed strike of the fault zone.
3
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 2. PDF (Probability Density Function) analysis, calculated using 30 days of the ground motion continuous recordings. From top to bottom, the graphs represent the PDF of EW, NS and Vertical components. The red line indicates the 95 percentiles of the probability density function. The density of probability is indicated by the grey shaded palette. The solid lines represent the upper and lower models of the Earth’s noise after Peterson (1993). Right panels: HV spectral ratios and corresponding ± 3σ boundaries, for data samples collected at day- and night times (top and bottom panels, respectively).
obtained a dataset of 82 earthquakes, which were located using the 1-D velocity model of the INGV seismic monitoring system (Table 1) and Hypoellipse code (Lahr, 1999). Only 5 hypocenters have a high-quality location, r.m.s. less than 0.19 s., with horizontal and vertical errors < 1 and 2 km, respectively, number of stations ≥ 5 and azimuthal gap less
than 180°. Fig. 4 shows these earthquakes hypocenters (black circles), located both in the central part of the study area (North-East of Carbonia town) and in the Sant’Antioco island, South-Eastern off-shore; hypocentral depths range from 2.5 to 18 km (Table 2). The green circles represent the remaining 77 earthquakes’ hypocenters with location
Fig. 3. Sensitivity maps at day- (left) and night-time (right) for the Sulcis temporary network, assuming a hypothetical hypocentral depth of 5 km. The dashed black box indicates the Sotacarbo Fault Lab area. 4
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
the Matzaccara Quaternary basin (a part of the wider Sulcis basin) where the Sotacarbo Fault Lab is going to be installed. In May 2016, we installed the Sotacarbo Fault Lab permanent seismic network that consists of 5 seismic stations, 4 equipped with three components intermediate period (Lennartz 3D/5 s) seismic sensors and 1 with a broadband sensor (Nanometrics Trillium compact 120 s), operating in the restricted area around the site (cyan and red triangles in Fig. 1b). In Fig. 5 the sensitivity map of the Sotacarbo permanent seismic network is shown. During a three-month period in 2017, up to 11 additional seismic stations were operating to temporarily obtain a denser temporary sampling in the Sotacarbo Fault Lab site. During the recording period of dense monitoring (60 days), seismicity was not detected in the study area and its close surroundings.
Table 1 1-D velocity model used for earthquakes’ location shown in Fig. 4. This model is the mean velocity model of the Italian territory, used by INGV during the realtime seismic survey of the Italian territory (Bollettino Sismico Italiano (BSI), 2019). P-wave velocity (km/s)
Depth (km)
Vp/Vs
5.0 6.5 8.05
0 11.1 38
1.73 1.73 1.73
errors ranging from 0.3 to 10 and 0.5 to 14.1 km, for horizontal and vertical coordinates. Only for the earthquakes located with higher accuracy, we computed the local magnitude (ML) on horizontal components of ground displacement by convolving the instrument-corrected seismograms with the instrumental response of the Wood-Anderson seismograph (Table 2), following Hutton and Boore (1987), with errors estimated through the robust Huber mean. The Magnitude ML computed for these 5 earthquakes ranging from 0.42 and 1.29.
4. Seismic imaging method and results For imaging the subsurface structure around the Sotacarbo Fault Lab site, and define the position and geometry of the fault, we performed a passive experiment with ambient noise tomography. The dispersion curve of Rayleigh waves for different inter-station paths, from the correlation properties of ambient noise, are inverted to derive 2D and 3D phase velocity maps at distinct frequencies (Shapiro and Campillo, 2004; Saygin and Kennett, 2012). Throughoutthis section, with ‘noise’ or ‘ambient noise’ we refer to the background Earth’s
3.1. The dense seismic network around the fault field lab After the definition of the seismic baseline in South-Western Sardinia, we focused our attention on the small area (about 13 Km2) of
Fig. 4. Map and vertical sections of the seismicity recorded and located during the experiment. Red triangles represent the stations of the first-large temporary network equipped with a 5-s period sensors. Blue triangles represent stations equipped with a broadband sensor. Black circles represent the 5 earthquakes with highquality location. The green circles represent the other 77 earthquakes with low location accuracy. Black box indicates the Sotacarbo Fault Lab area. 5
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Table 2 Origin time (T0), coordinates and magnitudes (ML) of the 5 high-quality location earthquakes. Year(YYYY)Month(MM)Day(DD) Hour(HH)Minute(MM)Second(SS.SS)
LATITUDE
LONGITUDE
DEPTH (km)
ML
ML Error
2015-03-11 2015-03-26 2015-04-28 2015-05-06 2015-05-06
38N58.47 39N10.38 39N03.81 39N10.06 39N11.67
8E27.01 8E32.41 8E46.38 8E35.47 8E34.18
17.41 2.87 10.31 2.41 2.81
1.07 1.03 1.24 0.42 1.29
+/- 0.27 +/- 0.26 +/- 0.17 +/- 0.26 +/-0.32
16:44:43.18 10:49:28.76 03:34:16.79 11:24:35.81 11:37:07.84
Fig. 5. Sensitivity maps at day- (left) and night-time (right) for the Sotacarbo Fault Lab permanent seismic network. The dashed black box indicates the Sotacarbo Fault Lab area.
vibrations induced by either natural or artificial causes. With ‘recording’ we refer to the measurement of the vertical component of ground velocity. Furthermore, we assume that the noise wavefield is temporally stationary and spatially isotropic. The temporal average of the narrow-band correlation coefficient of the ambient noise recorded by two stations is directly related to the phase velocity of seismic waves at that particular frequency (Aki, 1957; Chávez-Garcı́a et al., 2005). Hence, by repeating correlation estimates over a set of frequencies, one may derive a dispersion curve associated with the propagation path in between the two sites. Before deriving correlation estimates, noise recordings were subjected to a pre-conditioning, followed by an amplitude-normalization procedure aimed at removing the effect of transient, large-amplitude signals of either natural or artificial origin (Bensen et al., 2007). We then calculated the frequency-dependent correlation estimates for all the independent station pairs over the 0.2–1.5 Hz frequency band. For each propagation path, a final correlation function is obtained from the averaging of individual correlation estimates calculated over several, 6hour-long time intervals. Under the assumption that the vertical-component noise wavefield is dominated by fundamental-mode Rayleigh waves, we finally inverted individual correlation functions for the corresponding dispersion curves. For this procedure, we first calculated a trial estimate of the dispersion function using the direct-search method described in Saccorotti et al. (2001). This function is then used as a starting model for an iterative, non-linear least square inversion (Menke and Jin, 2015, Fig. 6).
stations, we obtained the frequency-dependent travel-times for the corresponding ray path. These times were finally inverted in a 2-D tomographic approach for deriving maps of phase velocity at different frequency bands (see Cauchie and Saccorotti, 2013). The tomography problem is solved by parameterizing the propagation medium according to rectangular cells, each of constant velocity. Under the further assumption that linear ray theory holds, the tomographic problem is represented by a set of equations in which individual travel times are expressed as a linear combination of ray length and slowness (the inverse of velocity) in individual cells. This set of linear equation is sparse and, generally, rank-deficient. For its solution, we used a least-squares inversion subjected to the constraints of solution smoothness and minimum distance from a starting, a-priori model (Menke, 1989). At each reference frequency, we used an a-priori homogeneous model whose velocity was obtained by the linear fit of the measured travel times against the corresponding ray lengths. The damping parameter, which controls the weight assigned to the solution smoothness constrain, was set equal to 20 following a trialand-error procedure based on the examination of the trade-off curves (Aster et al., 2005). Examination of the travel-time vs. distance relationships shows that the propagation times depend principally on distance and, as expected, they do increase for increasing station separation. The discrepancies with respect to the average velocity are not systematically related to any particular inter-station azimuth. This implies that the assumption of a spatially-isotropic noise wave-field is satisfied; in case of directional noise sources, in fact, travel time residuals would be markedly correlated with specific inter-station azimuths (Fig. 7). Phase velocity maps were calculated at a set of discrete frequencies spanning the 0.5 Hz–1.5 Hz frequency range, with 0.1 Hz spacing. The
4.1. Phase velocity maps From the phase-velocity dispersion curves obtained at any pair of 6
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 6. Inversion of the correlation function for a dispersion curve at a sample station pair. The panel at the left shows the trial dispersion curve derived from the procedure described in Saccorotti et al. (2001) (gray dashed line), which is then used as a starting model for the iterative, non-linear inversion, whose result is represented by the black line. Black dots are phase velocity estimates obtained from zero-crossings of the correlation function. The panel at the right shows the timeaveraged, experimental correlation function (gray line) and that predicted on the base of the best-fitting dispersion relationship. Black circles mark the zero crossings of the correlation curve.
emerge at the edges of the plain. These latter anomalies appear as isolated spots over the lowermost frequency bands (i.e., the largest depths), but progressively extend throughout both edges of the basin as frequency increases (i.e., at shallower depths). This observation supports the presence of N–S trending normal fault, which displaced the deepest part of the basin. 4.2. Profile along the linear array The ambient noise correlation analysis is then extended to stations from the linear profile shown in Fig. 10a. In this case, we grouped the correlation coefficients obtained over the same frequency band but for different inter-station distances, then conducting a direct search for that phase velocity which minimizes the difference between the observed and predicted correlation curves. For the direct search, we considered 100 phase velocity values spanning the 0.2 km/s - 3 km/s range. We conducted the analysis over groups formed by 4 consecutive stations (distances spanning the 120 and 360 m), and using the correlation coefficients calculated for all the independent inter-station distances associated with that group. The procedure was then iterated considering the next group, which was derived by translating the index of the first station by one unit. We thus obtained a set of 8 dispersion curves, which were combined to obtain a map where phase velocity is expressed as a function of frequency and offset along the profile. The resulting phase velocity section (Fig. 10b) shows two main discontinuities, positioned at about 400−500 m and 600−700 m along the profile. We associate the shallowest low and high velocities (Vr = Rayleigh phase velocity) to the Quaternary alluvial deposits (dark blue, VR ∼ 0.600 km/s) and pyroclastic products and ignimbrites of the lower and middle Miocene (yellow-green, VR ∼0.9 km/s), respectively. The underlying, high-velocity level highlighted in the NW margin of the profile (1< VR <1.2 km/s) could instead correspond to the magmatic products of the Lower Miocene that constitute the basement of the basin. For the same reasons outlined above, we didn’t perform any further inversion for transforming the phase velocity section of Fig. 10b into a corresponding depth section of the shear-wave velocity along the
Fig. 7. Travel times as a function of inter-station distance at the sample frequency F = 0.5 Hz. An average velocity (dashed line) is obtained from the slope of a straight line fitted to experimental data in a least-square sense. Colors refer to the orientation of individual station pairs, according to the color bar at the right.
model space was defined by a set of 8 × 5 rectangular cells with sides 650 m and 660 m along the EW and NS directions, to optimize image fidelity and model resolution. The ray crossing is adequate for ensuring a good resolution of velocities in the central cells, and still acceptable at peripheral cells with a lower hit counts, according to synthetic tests (see Fig. 8b, c, d). Fig. 9 shows the phase-velocity maps obtained from the tomographic inversions at the different target frequencies. Conversion of the phase velocity maps to a 3D image of shear-wave velocities would require an inversion of the dispersion functions associated with individual nodes of the maps. Such non-unique solutions strongly depend on inversion regularizations and starting models. Since our poor a priori knowledge of the shear-wave velocities at depth, we prefer to limit our analysis to the phase velocity maps. The Matzaccara Quaternary basin is identified as a low-velocity area bounded to the west and to the east by high velocity anomalies imputable to the lower and middle Miocene volcanic products, which 7
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 8. (a) Geometry of station pairs and corresponding seismic rays used in the tomographic inversion for phase velocity maps, (b) Number of rays crossing individual cells, (c) Input model for the checkerboard validation test. The model is defined by the alternance of cells with exhibiting +/- 10 % velocity anomalies with respect to the average velocity retrieved at that particular frequency. The synthetic travel-times were perturbed by adding random errors extracted from a normal distribution having zero mean and standard deviation equal to the standard deviation of the residual distribution between observed data and travel-times predicted by the starting uniform velocity model. (d) Image recovered after the tomographic inversion of synthetic travel-times calculated for the model shown in (c). For the central part of the network, the pattern of alternating positive and negative anomalies is well recovered. The color bars indicates the phase velocity in km/s.
profile. However, following the method proposed by Socco et al. (2017), one may obtain an approximate relationship between the observed wavelengths and the resolvable depths, thus getting a qualitative indication about the depth of the main velocity discontinuities. In this manner, we estimate that the fault (or fault segments) dissecting the basin would affect a crustal portion spanning the 50−450 m depth interval. The frequency band characterized by the lowermost velocities (dark blue in Fig. 10b) becomes progressively wider as one moves toward the SE margin of the profile, indicating a progressive thickening of the sedimentary / volcanoclastic deposits in-filling the basin. The maximum thickness of these deposits is inferred to be at least on the order of 150−200 m, and this interval may be taken as a rough estimate of the minimum displacement occurred along the faults bordering the basin. However, it cannot be excluded that the sedimentary cover gets even thicker SE of the imaged section, i.e. toward the center of the basin. In that case, the above estimate on the displacement along the fault would represent a lower bound on the actual cumulated slip.
4.3. The Matzaccara buried fault The good correspondence between the phase velocity maps of Fig. 9 and the phase velocity section along the linear profile (Fig. 10b) allows identification of an east-dipping fault that dislocates the shallow layers, and contributes to the formation of the Matzaccara Quaternary basin. The low-velocity anomaly coincident with the basin elongation strikes about N-NE, indicating that even the east-dipping fault imaged in the 2D section, is elongated along the same azimuth. This latter fault is composed by a superficial splay positioned between the surface and 150−200 m deep, located around 400−500 m away from station AT11 (Fig. 10b), which connects to a deeper segment spanning the 200−400 m depth interval, positioned at a distance of 600−700 m from station AT11. The fault displaces the volcanic formations for the depth range in which it was detected, but it is not possible to verify with certainty any more recent activity related to Quaternary deposits. 5. Discussion and conclusions The development and acceptance of CCUS requires careful analyses 8
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 9. Sample phase velocity maps at different reference frequencies. Color bar indicates phase velocity in km/s. Not resolved areas (hit counts <5) are masked.
of the hazards associated with the injection and storage activity, passing through the implementation of appropriate mitigation strategies and best practices of the risk management (White and Foxall, 2016). The caprock integrity and the modification of the stress field due to the confinement of large gas volumes are primary aspects for hazard assessment (Verdon et al., 2011; Verdon, 2014; Cantucci et al., 2016). Experiments at on-shore fields laboratories (such as the Sotacarbo Fault Lab, for instance) and pilot sites are therefore the most valuable strategy to understand how the deep systems evolve in response to fluid injection, learning by doing.
Under this perspective, the level of natural seismicity before the beginning of the activities is a pre-requisite for monitoring and identifying all the possible alterations in the subsurface during the fluid injection (Cappa and Rutqvist, 2011; Stork Anna et al., 2015, 2018). Seismological applications in CCUS pilot sites focus on the evaluation of the seismicity background before injection activities (Chen et al., 2015), the definition of the target structure and the identification of the faults interest it (Alcalde et al., 2013) and could be a useful tool in the study of geomechanical response of a site to the CO2 injection and to detect the seismic footprint of the CO2 leakage (Stork Anna et al., 2018). Past 9
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
Fig. 10. a) The seismic transect across the Sotacarbo Fault Lab used for deriving the 2-D phase velocity image shown in 10 b. Blue circles mark the stations, while the red boxes 1 and 2 represent the surface projection of the two branches (1 and 2) of the fault highlighted in the phase-velocity vertical section. The figure shows just a small portion of the Matzaccara basin. b) Rayleigh-wave phase velocities derived from correlation analysis. Color bar indicates phase velocity in km/s. Velocities are mapped in their dependence upon frequency and distance along the seismic profile. Note that frequency ranges of 2.4-1.6 Hz and 1.6-0.8 Hz correspond to depth intervals of 0–200 and 200−400 m, respectively. Black lines mark two evident lateral discontinuities likely attributable to the Matzaccara normal fault, with a geometry similar to the sketch in the inset. The number 1 and 2 are referred to the dashed boxes in Fig. 10a.
experiences evidenced the strict requirement of having correct baselines to interpret the microseismicity revealed during the exploitation activities (Stork Anna et al., 2015, 2018). The final aim is to use monitoring data in a proactive real time forecast system, adapted to a defined site-specific sensibility. One second key aspect for the safe storage of large CO2 volumes is the characterization of fault and fracture systems in the subsurface. Especially important is the individuation of faults quiescent since long time and capable of producing damaging seismic events. Examples of faults activation during fluid injection became evident in different energy methods and context (waste water disposal, Keranen et al., 2014; E.G.S., Kraft and Deichmann, 2014; E.G.S., Diehl et al., 2017; waste water disposal, E.O.R., hydrofracturing, C.C.S. E.G.S., McGarr et al.,
2015; C.C.S., Stork Anna et al., 2018). In several cases, the presence of these faults was unknown before the starting of the industrial injection activities, although they dissected the basement under the sedimentary cover. Seismological applications are designed to define the presence of faults at depth, but understanding if those faults are critically stressed and prone to fail is very challenging. Today, we still lack methods that clearly address the identification of critically-stressed faults among those identified by seismic imaging of the upper crust (Ellsworth, 2013; Buttinelli et al., 2016). Faults and fracture systems are also preferential paths for the CO2 plume (Chang and Bryant, 2009) and the definition of the correct asset of fractures in the subsurface is strategic to tune the location of injection and monitoring wells, and to identify possible gas escape pathways. 10
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al.
We have investigated these issues in the Sulcis basin (South-Western Sardinia), where the Sotacarbo Fault Lab will be installed. Several years of seismological surveys helped to define the seismicity background of the area. This task has been particularly difficult because the area is slowly deforming (Palano, 2015) and the seismic rate is low, requiring longer observation times for the definition of seismic baselines and identification of faults at depth. We performed a three-step baseline survey with a progressive focusing and densification around the target, starting from a regional scale. In fact, we ignored at the beginning the seismic level for a broad region of Sardinia, that was never the subject of a focused monitoring before. From the first, large-aperture network, we were able to identify a few earthquakes in the area with magnitude around ML -1.0 (Fig. 3), that constitutes an unprecedented result for this low seismic area. Earthquakes hypocenters are sparse and do not align on specific faults or structures. Hypocentral accuracy do not permit detailed investigations. Then, we focused on the target site with a composite network formed by a permanent backbone and extremely dense but shortduration temporary arrays. In this latter case, we did not detect any seismicity down to magnitude as low as the that indicated by the sensitivity study (ML -1 / -1.5, within the top 5 km depth). The other strategic part of the experiment was dedicated to the reconnaissance and characterization of faults in the subsurface. Given the absence of seismicity, passive seismic methods like travel-time tomography were not applicable. We therefore produced a detailed imaging of the subsurface using ambient noise tomography (Lin et al., 2013). We have thus been able to identify an east-dipping normal fault in the 400 m below the surface (Fig. 10b). Although the lack of resolution over the shallowest tens of meters prevents from estimating the last episode of slip on the fault, its last activity likely occurred during Quaternary, according with the displaced shallow sediments. The fault trends N-NE, and it possibly extends for at least 2−3 km or more, outside of the resolution provided by the tomographic model (Fig. 9). Throughout the two-year-long seismic monitoring of the area, no earthquakes have been detected in association to this fault, down to the detection threshold of ML -1.0.
References Aki, K., 1957. Space and time spectra of stationary waves with special reference to microtremors. Bull. Earthquake Res. Inst. Univ. Tokyo 35, 415–456. Alcalde, J., Martí, D., Calahorrano, A., Marzan, I., Ayarza, P., Carbonell, R., Juhlin, C., Pérez-Estaún, A., 2013. Active seismic characterization experiments of the Hontomín research facility for geological storage of CO2, Spain. Int. J. Greenh. Gas Control. 19, 785–795. https://doi.org/10.1016/j.ijggc.2013.01.039. ISSN 1750-5836. Arts, R., Eiken, O., Chadwick, A., Zweigel, P., van der Meer, L., Zinszner, B., 2004. Monitoring of CO2 injected at Sleipner using time-lapse seismic data. Energy 29 (9–10), 1383–1392. https://doi.org/10.1016/j.energy.2004.03.072. ISSN 03605442. Aster, R., Beaudoin, B., Hole, J., Fouch, M., Fowler, J., James, D., 2005. IRIS Seismology Program Marks 20 Years of Discovery. Eos, Transactions American Geophysical Union. Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M., Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad‐band surface wave dispersion measurements. Geophys. J. Int. 169, 1239–1260. https://doi.org/10.1111/j.1365-246X.2007.03374.x. Boatwright, John, 1978. Detailed spectral analysis of two small New York state earthquakes. Bull. Seismol. Soc. Am. 68 (4), 1117–1131. Bollettino Sismico Italiano (BSI) - Gruppo di Lavoro Sequenza Centro Italia (2019). Rapporto Bollettino Sismico Italiano sulla revisione della sequenza sismica del centro Italia 24 agosto 2016 - 31 agosto 2018. Bommer, J.J., Crowley, H., Pinho, R., 2015. A risk-mitigation approach to the management of induced seismicity. J. Seismol. 19 (2), 623–646. Boot-Handford, M.E., Abanades, J.C., Anthony, E.J., Blunt, M.J., Brandani, S., Mac Dowell, N., Fernandez, J.R., Ferrari, M.C., Gross, R., Hallett, J.P., Haszeldine, R.S., Heptonstall, P., Lyngfelt, A., Makuch, Z., Mangano, E., Porter, R.T.J., Pourkashanian, M., Rochelle, G.T., Shah, N., Yao, J.G., Fennell, P.S., 2014. Carbon capture and storage update. Energy Environ. Sci. 7, 130–189. Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes. J. Geophys. Res. 75 (26), 4997–5009. https://doi.org/10.1029/ JB075i026p04997. Buttinelli, M., Improta, L., Bagh, S., Chiarabba, C., 2016. Inversion of inherited thrusts by wastewater injection induced seismicity at the Val d’Agri oilfield (Italy). Sci. Rep. https://doi.org/10.1038/srep37165. Cantucci B., Buttinelli M, Procesi M., Sciarra A., Anselmi M. (2016) Algorithms for CO2 storage capacity estimation: review and case study. Book Chapter into "Geologic Carbon Sequestration: Understanding Reservoir Concepts", Section 2 “Geological Characterisation of storage sites”. published by Springer under the joint editorial handling of Prof. T. N. Singh (Group Lead, Rock Mechanics, IIT Bombay) DOI https:// doi.org/10.1007/978-3-319-27019-7. Cappa, F., Rutqvist, J., 2011. Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO2. Int. J. Greenh. Gas Control. 5 (2), 336–346. Cauchie, L., Saccorotti, G., 2013. J. Seismol. 17, 335. https://doi.org/10.1007/s10950012-9323-6. Chang, Kyung Won, Bryant, Steven, 2009. The effect of faults on dynamics of CO2 plumes. Energy Procedia 1, 1839–1846. https://doi.org/10.1016/j.egypro.2009.01. 240. Chávez-Garcı́a, F.J., Rodrı́guez, M., Stephenson, W.R., 2005. An alternative approach to the SPAC analysis of microtremors: exploiting stationarity of noise. Bull. Seism. Soc. Am. 95, 277–293. https://doi.org/10.1785/0120030179. Chen, T., Lianjie, Huang, 2015. Seismicity characterization around the Farnsworth field site for combined large-scale CO2 storage and EOR. Int. J. Greenh. Gas Control. 33, 63–68. https://doi.org/10.1016/j.ijggc.2014.11.024. ISSN 1750-5836. Diehl, T., Kraft, T., Kissling, E., Wiemer, S., 2017. The induced earthquake sequence related to the St. Gallen deep geothermal project (Switzerland): fault reactivation and fluid interactions imaged by microseismicity. J. Geophys. Res. Solid Earth 122, 7272–7290. https://doi.org/10.1002/2017JB014473. Dipietro, J., et al. “The Role of Naturally-Occurring CO2 Deposits in the Emergence of CO2 Enhanced Oil Recovery” http://co2conference.net/pdf/1.2-Slides_ DiPietroCO2Sources2011CO2FloodingConf.pdf. Ellsworth, W.L., 2013. Injection-induced earthquakes. Science 341 (6142), 1225942. https://doi.org/10.1126/science.1225942. Förster, Andrea, Norden, Ben, Zinck-Jørgensen, Kim, Frykman, Peter, Kulenkampff, Johannes, Spangenberg, Erik, Erzinger, J.örg, Zimmer, Martin, Kopp, J.ürgen, Borm, G.ünter, Juhlin, Chris, Cosma, Calin-Gabriel, 2006. Suzanne Hurter; Baseline characterization of the CO2SINK geological storage site at Ketzin. Germany. Environmental Geosciences 13 (3), 145–161. https://doi.org/10.1306/eg. 02080605016. Grigoli, F., Cesca, S., Priolo, E., Rinaldi, A.P., Clinton, J.F., Stabile, T.A., Dost, B., Fernandez, M.G., Wiemer, S., Dahm, T., 2017. Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: a European perspective. Rev. Geophys. 55, 310340. https://doi. org/10.1002/2016RG000542. Gruppo di Lavoro, 2004. Redazione della mappa di pericolosità sismica prevista dall’Ordinanza PCM 3274 del 20 marzo 2003. Rapporto conclusivo per il Dipartimento della Protezione Civile. INGV, Milano-Roma aprile 2004, 65 pp + 5 enclosures. Hutton, Boore, 1987. The ML scale in southern California. Bull. Seismol. Soc. Am. 77 (6), 2074–2094. Kanamori, H., 1977. The energy release in great earthquakes. J. Geophys. Res. 82 (20), 2981–2987. https://doi.org/10.1029/JB082i020p02981.
CRediT authorship contribution statement M. Anselmi: Conceptualization, Methodology, Data curation, Writing - original draft. G. Saccorotti: Data curation, Methodology, Writing - original draft. D. Piccinini: Methodology, Data curation. C. Giunchi: Methodology, Data curation. M. Paratore: Data curation, Software. P. De Gori: Software, Methodology, Data curation. M. Buttinelli: Conceptualization, Validation. E. Maggio: Resources, Funding acquisition. A. Plaisant: Resources, Supervision. C. Chiarabba: Conceptualization, Writing - original draft, Supervision. Declaration of Competing Interest None. Acknowledgements This study has been carried out within the “Research on Electric System” project funded by the Italian Ministry of Economic Development. The Sotacarbo Fault Lab is under construction within the “Centre of Excellence on Clean Energy” project (CUP: D83C17000370002) funded by the Regional Government of Sardinia (FSC 2014–2020). The data from the station SU26, belonging to the Sotacarbo Fault Lab seismic network, is currently archived in EIDA waveform repository (eida.rm.ingv.it). We are grateful to Gianfranco Colasanti, to the “Rete Sismica Mobile” group at INGV and to the "Corpo Forestale della Sardegna" for their helpful technical and logistics support. We thank Alberto Pettinau of Sotacarbo S.p.A. for his helpful contribution in the review of the paper. 11
International Journal of Greenhouse Gas Control 95 (2020) 102974
M. Anselmi, et al. Keranen, K.M., Weingarten, M., Abers, G.A., Bekins, B.A., Ge, S., 2014. Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection. Science 345 (6195), 448–451. Kraft, T., Deichmann, Nicholas, 2014. High-precision relocation and focal mechanism of the injection-induced seismicity at the Basel EGS. Geothermics 52, 59–73. https:// doi.org/10.1016/j.geothermics.2014.05.014. ISSN 0375-6505. Lahr, J.C., 1999, revised 2012, HYPOELLIPSE: a computer program for determining local earthquake hypocentral parameters, magnitude, and first-motion pattern: U.S. Geological Survey Open-File Report 99–23, version 1.1, 119 p. and software, available at https://pubs.usgs.gov/of/1999/ofr-99-0023/. Lin, Fan-Chi, Li, Dunzhu, Clayton, Robert W., Hollis, Dan, 2013. High-resolution 3D shallow crustal structure in Long Beach, California: application of ambient noise tomography on a dense seismic array. Geophysics 78 (4), Q45–Q56. https://doi.org/ 10.1190/geo2012-0453.1. Liu, H.J., Were, P., Li, Q., Gou, Y., Hou, Z., 2017. Worldwide Status of CCUS Technologies and Their Development and Challenges in China Geofluids. Article ID 6126505, 25 pages. https://doi.org/10.1155/2017/6126505. McGarr, A., Bekins, B., Burkardt, N., Dewey, J., Earle, P., Ellsworth, W., Ge, S., Hickman, S., Holland, A., Majer, E., Rubinstein, J., Sheehan, A., 2015. Coping with earthquakes induced by fluid injection. Science 347, 830–831. https://doi.org/10.1126/science. aaa0494. Menke, W., 1989. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press, Inc., New York, pp. 289pp. Menke, W., Jin, G., 2015. Waveform fitting of cross spectra to determine phase velocity using Aki’s formula. Bull. Seism. Soc. Am. 105, 1619–1627. https://doi.org/10.1785/ 0120140245. Munafò, I., Malagnini, L., Chiaraluce, L., 2016. On the relationship between mw and ML for small earthquakes. Bull. Seismol. Soc. Am. 106 (5), 2402–2408. https://doi.org/ 10.1785/0120160130. Nicol, A., Carne, R., Gerstenberger, M., Christophersen, A., 2011. Induced seismicity and its implications for CO2 storage risk. Energy Procedia 4, 3699–3706. Palano, M., 2015. On the present-day crustal stress, strain-rate fields and mantle anisotropy pattern of Italy. Geophys. J. Int. 200 (2), 969–985. https://doi.org/10.1093/ gji/ggu451. Ringrose, P.S., Mathieson, A.S., Wright, I.W., Selama, F., Hansen, O., Bissell, R., Saoula, N., Midgley, J., 2013. The in Salah CO2 storage project: lessons learned and knowledge transfer. Energy Procedia 37, 6226–6236. Rovida, A., Locati, M., Camassi, R., Lolli, B., Gasperini, P. (Eds.), 2016. CPTI15, the 2015 Version of the Parametric Catalogue of Italian Earthquakes. Istituto Nazionale di Geofisica e Vulcanologia. Saccorotti, G., Almendros, J., Carmona, E., Ibáñez, J.M., Del Pezzo, E., 2001. Slowness anomalies from two dense seismic arrays at Deception Island volcano, Antarctica. Bull. Seismol. Soc. Am. https://doi.org/10.1785/0120000073. Saygin, E., Kennett, B.L.N., 2012. Crustal structure of Australia from ambient seismic noise tomography. J. Geophys. Res. 117, B01304. https://doi.org/10.1029/ 2011JB008403. Schorlemmer, D., Mele, F., Marzocchi, W., 2010. A completeness analysis of the National Seismic Network of Italy. J. Geophys. Res. 115, B04308. https://doi.org/10.1029/ 2008JB006097. Shapiro, N.M., Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophys. Res. Lett. 31, L07614. https://doi. org/10.1029/2004GL019491.
Socco, L.V., Comina, C., Anjom, F.K., 2017. Time-average velocity estimation through surface-wave analysis: Part 1 — S-wave velocity. Geophysics 82, U49–U59. https:// doi.org/10.1190/GEO2016-0367.1. Stork Anna, L., Verdon, James P., Kendall, J.-Michael, 2015. The microseismic response at the in-salah Carbon Capture and Storage (CCS) site. Int. J. Greenh. Gas Control. 32, 159–171. https://doi.org/10.1016/j.ijggc.2014.11.014. ISSN 1750-5836. Stork Anna, L., Allmark, Claire, Curtis, Andrew, Kendall, J.-Michael, White, Don J., 2018. Assessing the potential to use repeated ambient noise seismic tomography to detect CO2 leaks: application to the Aquistore storage site. Int. J. Greenh. Gas Control. 71, 20–35. https://doi.org/10.1016/j.ijggc.2018.02.007. ISSN 1750-5836. Stucchi, Massimiliano, Meletti, Carlo, Montaldo, Valentina, Crowley, Helen, Calvi, Gian Michele, Boschi, Enzo, 2011. Seismic hazard assessment (2003–2009) for the italian building code. Bull. Seismol. Soc. Am. 101 (4), 1885–1911. https://doi.org/10.1785/ 0120100130. Takagishi, Makiko, Hashimoto, Tsutomu, Horikawa, Shigeo, KinichiroKusunose, Ziqiu Xue, Susan, D.Hovorka, 2014. Microseismic Monitoring at the Large-scale CO2 Injection Site. Cranfield, MS, U.S.A. https://doi.org/10.1016/j.egypro.2014.11.476. Ugalde, A., Villaseñor, A., Gaite, B., Casquero, S., Martí, D., Calahorrano, A., Marzán, I., Carbonell, R., Estaún, A.P., 2013. Passive seismic monitoring of an experimental CO2 geological storage site in Hontomín (Northern Spain). Seismol. Res. Lett. 84 (1), 75–84. Van Thienen-Visser, K., Breunese, J.N., 2015. Induced seismicity of the Groningen gas field: history and recent developments. Lead. Edge 34https://doi.org/10.1190/ tle34060664.1. 664–666, 668–668, 670–671. Verdon, J.P., 2014. Significance for secure CO2 storage of earthquakes induced by fluid injection. Environ. Res. Lett. 9 (6). Verdon, J.P., Kendall, J.-M., White, D.J., Angus, D.A., 2011. Linking Microseismic Event Observations with Geomechanical Models to Minimise the Risks of Storing CO2 in Geological Formations. https://doi.org/10.1016/j.epsl.2011.02.048. Vilarrasa, V., Carrera, J., 2015. Geologic carbon storage is unlikely to trigger large earthquakes and reactivate faults through which CO2 could leak. Proc. Natl. Acad. Sci. 112, 5938–5943. Vilarrasa, Victor, Carrera, Jesus, Olivella, Sebastià, Rutqvist, Jonny, Laloui, Lyesse, 2019. Induced seismicity in geologic carbon storage. Solid Earth 10 (871–892). https://doi. org/10.5194/se-10-871-2019. White, J.A., Foxall, W., 2016. Assessing induced seismicity risk at CO2 storage projects: recent progress and remaining challenges. Int. J. Greenh. Gas Control. 49, 413–424. Woessner, Jochen, Laurentiu, Danciu, Giardini, Domenico, Crowley, Helen, Cotton, Fabrice, Grünthal, Gottfried, Valensise, Gianluca, Arvidsson, Ronald, Basili, Roberto, Demircioglu, Mine, Hiemer, Stefan, Meletti, Carlo, Musson, RogerW., Rovida, Andrea, Sesetyan, Karin, Stucchi, Massimiliano, 2015. The 2013 European Seismic Hazard Model: key components and results. Bull. Earthq. Eng. 13, 1–44. https://doi. org/10.1007/s10518-015-9795-1. Zang, A., Oye, Volker, Jousset, Philippe, Deichmann, Nicholas, Gritto, Roland, McGarr, Art, Majer, Ernest, Bruhn, David, 2014. Analysis of induced seismicity in geothermal reservoirs -An overview. Geothermics 52, 6–21. https://doi.org/10.1016/j. geothermics.2014.06.005. ISSN 0375-6505. Zheng, S., Sun, X., Song, X., Yang, Y., Ritzwoller, M.H., 2008. Surface wave tomography of China from ambient seismic noise correlation. Geochem. Geophys. Geosyst. 9, Q05020. https://doi.org/10.1029/2008GC001981. Zoback, M.D., Gorelick, S.M., 2012. Earthquake triggering and large-scale geologic storage of carbon dioxide. Proc. Natl. Acad. Sci. 109, 10164–10168.
12