4D Travel-Time Tomography as a Tool for Tracking Fluid-Driven Medium Changes in Offshore Oil–Gas Exploitation Areas

Abstract

The monitoring of rock volume where offshore exploitation activities take place is crucial to assess the corresponding seismic hazard. Fluid injection/extraction operations generate a pore fluid pressure perturbation into the volume hosting the reservoir which, in turn, may trigger new failures and induce changes in the elastic properties of rocks. Our purpose is to evaluate the feasibility of reconstructing pore pressure perturbation diffusion in the host medium by imaging the 4D velocity changes using active seismic. We simulated repeated active offshore surveys and imaged the target volume. We constructed the velocity model perturbed by the fluid injection using physical modeling and evaluated under which conditions the repeated surveys could image the velocity changes. We found that the induced pressure perturbation causes seismic velocity variations ranging between 2–5% and 15–20%, depending on the different injection conditions and medium properties. So, in most cases, time-lapse tomography is very efficient in tracking the perturbation. The noise level characterizing the recording station sites is a crucial parameter. Since we evaluated the feasibility of the proposed 4D imaging strategy under different realistic environmental and operational conditions, our results can be directly applied to set up and configure the acquisition layout of surveys aimed at retrieving fluid-induced medium changes in the hosting medium. Moreover, our results can be considered as a useful starting point to design the guidelines to monitor exploitation areas.

1. Introduction

The exploration and exploitation of geo-resources play an increasingly significant role in energy production. However, during industrial operations, the extraction or re-injection of fluids may create new fractures and/or alter the frictional condition of existing faults, triggering low to moderate earthquakes [1]. Therefore, monitoring areas of geo-resource exploitation is fundamental since it allows us to follow the evolution in space, time, and magnitude of seismicity, and to reschedule or suspend industrial activities [2].

Fluid injection/extraction operations generate a pore fluid pressure perturbation into the volume hosting the reservoir. Even if the production reservoir is only a few hundred meters thick, the pore pressure perturbation diffusion can extend several kilometers [3]. In fact, induced seismicity, with a magnitude that can be larger than M 3, has been reported to occur kilometers to tens of kilometers away from injection locations via pore pressure increase (e.g., [4,5]).

The pore fluid perturbation, whose features depend on the injection source time history, expands according to the characteristics of the medium and of the fluid itself, and modifies the effective pressure state in the rock volume. In turn, to the first order, the effective pressure affects the elastic properties of rocks, i.e., the seismic P- and S-wave velocities. Using rock physics approaches, the diffusive propagation of the pore pressure perturbation can be reconstructed by analyzing the tomographic images of seismic velocity in consecutive time epochs (time-lapse travel-time tomography) [6]. This application allows identifying and tracking any changes or anomalies possibly induced by pore fluid pressure diffusion that, in turn, can modify the stress state and trigger a new failure.

Time-lapse tomography consists of processing the data acquired during repeated active or passive seismic surveys and comparing the resulting tomographic images. In this case, each image can be considered as a snapshot of the investigated medium at the time (epoch) data were collected. The acquisition layout is kept fixed to not introduce a condition that can modify the resolved volumes between surveys. This condition is obtained in the case of active seismic data, but not with the use of passive seismicity, as the position of sources may change. In order to mitigate the seismic hazard in highly vulnerable areas, the 4D travel-time tomographic technique is now largely applied in volcanic, tectonic, and exploitation environments [7,8,9] by using passive seismic sources.

Gunasekera et al. [7] applied 4D travel-time tomographic analysis to the micro-seismicity occurring at The Geysers geothermal field (California, USA) recorded by a dense, permanent network. The authors performed tomographic inversions for the 3D velocity structure of the reservoir for five time epochs between 1991 and 1998. They inferred that the retrieved anomalies in the velocity fields are caused by exploitation. In particular, water depletion and pressure reduction are the dominant processes in the volume embedding the reservoir.

Valoroso et al. [8] used local earthquake tomography to investigate spatio-temporal anomalies of the crustal velocity structure in Val d’Agri (southern Apennines, Italy). They observed changes of velocity images between 3 and 6 km depth, where fluids are contained in a highly fractured medium. These changes identified that transient variations of pore fluid pressure related, in turn, to a rise in the seismicity rate. These results allow the authors to infer a mechanism of reservoir seismicity induced by fluid pressure increase and pore pressure diffusion.

Calò and Dorbath [9] presented the results of a 4D seismic tomography obtained by inverting the P-wave arrival times of seismic events recorded in 2003 during the stimulation (hydraulic injections) of a well (depth of about 5 km) in the Soultz-sous-Forêts (northeastern France) geothermal fields. The events were divided into 13 subsets identified with respect to the injection scheme. The authors observed a lack of large low-P-wave velocity anomalies during the stimulation. They proposed that a large network of faults represents the main path of the injected water, avoiding the accumulation of effective stresses in the reservoir close to the well. Moreover, the water injection affects a large region and activates structures far from the wells.

After the 2012 Emilia earthquake, the Italian Ministry of Economic Development (MISE) introduced the guidelines (https://unmig.mise.gov.it/index.php/it/sicurezza/geomonitoraggi/linee-guida) for monitoring induced seismicity possibly generated by the anthropic activities of onshore cultivation, re-injection, and storage of hydrocarbons. The guidelines prescribe that monitoring must be performed before and during the period of exploitation activity. Since there are many concessions offshore, no guidelines for reservoir monitoring have been established yet and few are the case studies in these kinds of sites. De Landro et al. [10] proposed an optimization study for seismic monitoring networks based on the MISE guideline requirements in terms of detection threshold (i.e., M = 0) and location accuracy of seismic events (i.e., location errors of few hundred meters). The results show that the station density and the noise level represent the crucial parameters for the seismic network performances in locating micro-earthquakes and determining their magnitude. In this work, we considered the network configuration of the work of De Landro et al. [10], optimized to meet the MISE guidelines requirements.

The offshore areas in Italy where hydrocarbon research and exploitation are currently performed are mainly located in the Adriatic, Ionian, and Tyrrhenian Seas (https://unmig.mise.gov.it/images/cartografia/carta-titoli-31-luglio-2020.pdf). As a case study, we focused our analysis on the Adriatic Sea areas since they host the majority of oil–gas offshore exploitation fields.

Although we chose a specific Italian area, we believe that the proposed methodology, which allows for the tracking of medium velocity changes even before the occurrence of induced earthquakes, can be of general interest for the mitigation of seismic hazards in exploited areas.

In this work, we evaluate the capability of 4D P-wave seismic tomography to provide images of the volume containing the reservoir in order to detect and track the pore pressure perturbation induced in the host medium by the fluid injection. To avoid being limited to the natural (or induced) seismicity rate, which might not be sufficient for a seismic tomography, we simulated several active offshore seismic surveys (see the scheme in Figure 1), in different time epochs, after the start of the continuous injection. The design of the source–station acquisition layout was optimized by considering the spatial extension of the target volume and by building realistic phase catalogues based on reliable attenuation laws. Then, we modeled the pore pressure perturbation diffusion in the reservoir volume after a continuous fluid injection in order to retrieve the perturbed velocity model. Finally, we evaluated in which of the simulated operational and environmental conditions the active offshore seismic survey is capable of imaging the velocity changes due to the fluid injection and, therefore, able to properly monitor the conditions of the volume hosting the reservoir.

In the following, we first describe the structure of the feasibility study in terms of setting up the active survey parameters, construction of the perturbed velocity model, and the 3D tomographic technique. Then, we analyze the obtained results, evaluating the method performance in reconstructing the shape, amplitude, and location of velocity anomalies in the different time epochs. Finally, we discuss the main results and their implication, and highlight the main findings in the conclusion.

2. Materials and Methods

We investigated a volume of 30 × 30 × 10 km³ embedding the injection point, located at 3 km depth.

We performed the feasibility study by simulating repeated active seismic surveys after 1 day, 1 week (i.e., 7 days), and 1 month (i.e., 30 days) from the beginning of the continuous water injection.

We followed the steps shown in Figure 2:

• We defined the source–station layout. We considered a seismic network with 24 ocean-bottom seismometers (OBS). The number of stations and their spacing were designed to ensure a high resolution in the volume embedding the reservoir (for further details, see [10]). We fixed the number of active sources and their location according to the dimension and depth of the target volume.
• We introduced an attenuation law. The pressure waves emitted by the air gun propagate through the water and reach the sea bottom interface, where they generate elastic waves. The attenuation law allows us to properly evaluate the amplitude of the transmitted P-wave and build reliable phase catalogues in different conditions of sea bottom interface characteristics and noise level at the stations.
• We constructed the velocity model perturbed by fluid pressure diffusion. The pore pressure perturbation induced by the fluid injection is modeled using the Biot [11] equations of poroelasticity. Then, we applied the Eberhart-Phillips et al. [12] empirical relations to convert the pore pressure perturbation in seismic velocity perturbations. We constructed the models in different conditions of medium diffusivity and injection rates.
• In this perturbed velocity model, we computed the P-wave travel times. In order to have a more reliable catalogue, we added noise and assigned a weight to arrival time estimates based on the signal to noise ratio (hereinafter S/N ratio). The S/N ratio was used to perform a selection of the total amount of phases (S/N ratio greater than 2). We repeated this step for each simulated scenario with two different noise levels and sea bottom characteristics (i.e., we built four different catalogues).
• We inverted the simulated datasets with a linearized iterative scheme, obtaining 3D tomographic velocity images. Finally, we evaluated the feasibility in terms of differences between the retrieved and the “true” velocity anomalies in the target volume for each simulated scenario.

Details about the previous steps follow in the next sections.

2.1. Setting up the Parameters of an Ideal Active Offshore Seismic Survey

Compressed air sources, referred to as air guns, have been the dominant marine sound source since the 1960s for offshore seismic exploration (for further details on air guns, see [13]). In our study, we considered batteries of air guns dragged by a vessel as the active sources in our repeated seismic experiments.

In order to optimize the layout of the active seismic experiment, we preliminarily defined the volume that has to be investigated by the tomography. The injection point is assumed to be at 3 km of depth. According to our numerical simulations (see the following section), in the analyzed scenarios, the pore pressure perturbation can spread up to a maximum depth of 6 km, so we are interested in reaching a good resolution down to this depth.

Then, we traced the seismic ray paths in a 2D layered medium to define the minimum source–receiver distance, in order to properly image the depth of interest. We used the 1D velocity model of Venisti et al. [14], optimized for the Adriatic area. The resulting maximum distance is about 45 km (see Figure S1 in Supplementary Material, hereinafter SM).

We introduced an attenuation law to assign an error and a weight, based on the S/N ratio, to each recorded first P arrival time, and therefore to obtain a phase catalogue. Moreover, the attenuation law allows us to evaluate the signal amplitude expected at the maximum distance (i.e., 45 km) and to verify that this amplitude is sufficiently larger than the considered noise levels.

Considering the expected attenuation in a homogeneous and elastic medium, the law describing the amplitude decay for the P-wave ($A_p$) with distance can be expressed as

$$\log A_p\left(R,\omega \right)=\log A_0-k\log R-0.43\frac{{\omega }_{0}R}{2\alpha Q_p}$$

(1)

where $A_0$ is the amplitude at the source, $k$ is the distance attenuation parameter, $R$ is the source–receiver distance, ${\omega }_0$ is the wavelet dominant frequency, $\alpha$ is the P-wave velocity, and $Q_p$ is the P-wave quality factor. This equation expresses in logarithmic form the far-field P-wave amplitude radiated by a point source in a homogeneous elastic medium that underwent weak anelastic attenuation ([15], ch.5). It has the classical form of a ground motion prediction equation, where the terms depending on distance and frequency are given by the theory of a spherical radiating wave in a weakly attenuating medium characterized by a uniform seismic velocity $V_p$ and quality factor $Q_p$.

When the sound pressure wave emitted by the air gun reaches the sea bottom, it generates an elastic wave in the sub-marine sediments. In order to evaluate the amplitude of the transmitted P-wave, we calculated the water particle velocity corresponding to the pressure at the air gun source [16]. Then, the amplitude of the P-wave transmitted to the sediment layer was determined by multiplying the particle velocity by the transmission coefficient [17]. The obtained value constitutes the term $A_0$ in the attenuation law.

We considered air guns developing a pressure of 260 dB, referred to as a pressure of 1 μPa [18]. We used two transmission coefficients associated with different rheological behaviors of the seafloor (soil characteristics in

Table 1. Details regarding the input parameters in the rock modeling.

Z (km)	Porosity	Clay Content	P Wave Velocity (Km/s)	Shear Modulus (GPa)
1	0.4	0.3	4.1	35
5	0.3	0.3	5.9	35
7	0.2	0.3	6.1	35

): 0.54 named soft (SF), and 0.26 named hard (HD) [17]. We fixed the P-wave velocity at 5500 m/s [14] and the $Q_p$ at 200. This is a typical value of the attenuation factor in the first kilometer of the seafloor [19]. We expect the $Q_p$ to be greater than 200 below the first kilometers, so our choice is conservative since the wave attenuation effects reduce with increasing $Q$.

For a reliable definition of the S/N ratio, we considered two different noise levels: 10⁻⁷ m/s, named the maximum noise level (MAX), and 10⁻⁸ m/s, named the minimum noise level (MIN). These values were extracted from the power spectral density function obtained in different marine environments (i.e., Campi Flegrei and Pozzuoli Bay, Italy, [20]; Norwegian Sea, [21]; North Atlantic and Tyrrhenian Sea, [22]). It is worth noting that noise levels ranging from 10⁻⁷ to 10⁻⁶ m/s were observed at shallow-depth water in Pozzuoli Bay, an area that is characterized by high anthropic seismic noise, which for our purposes represents the worst-case scenario.

Considering the parameters chosen in our simulations (the red curves in Figure 3a,b), in the worst-case scenario (i.e., HD MAX), the signal amplitude ratio is higher than that of the noise up to a 40-km distance, just below the required one.

In order to account for a larger medium parameter variability, we also evaluated the attenuation law for a wider range of the $V_p$ and the $Q_p$ values (Figure 3a,b). We can see that the distance up to which the S/N ratio is higher than one varies over tens of km, in correspondence with the variation in the noise level and in the velocity and attenuation parameters.

It is clear that the parameters in the attenuation law have to be fixed reasonably to properly design the active survey. The MISE guidelines, for example, require a 3D seismic survey before the start of the industrial operations. This survey is likely also necessary in the offshore monitoring cases, and it can be used to determine suitable values for these parameters for the area of interest.

Finally, we fixed the source–receiver layout and reconstructed the 3D ray path distribution in order to evaluate the ray coverage and resolution. The station geometry is designed according to the MISE guidelines for proper reservoir monitoring in an offshore area [10]. We made this choice to reduce the parameters of the simulation. As for the source geometry, we used a spiral pattern, considered as one of the most efficient shot geometries in sub-surface seismic imaging applications [23]. This allowed us to image the small volume at the center of the spiral with a high resolution. The path starts above the center of the area covered by the stations and proceeds with a rate of one shot every 125 m for a period of three days. We therefore plan about 5000 shots located at a maximum distance of about 45 km from the well injection location and a total ship travel path of about 600 km. These characteristics reproduce the real offshore active survey SERAPIS [24]. The retrieved ray paths confirm that the target volume, up to 6 km of depth, is adequately sampled by the seismic rays (see Figure 3c,d).

2.2. Construction of the Perturbed Velocity Model

The differential pressure, also called effective pressure, is the difference between the confining and the pore pressure (see, for instance, [25]).

$$P=P_{conf}-P_{pore}$$

(2)

To the first order, this difference affects the elastic properties of rocks. The P-wave (compressional) and S-wave (shear) velocities increase with the increase in the effective pressure because of the closing of cracks, flaws, and grain boundaries, which elastically stiffen the mineral frame of the rocks. Conversely, an increase in pore pressure opens crack-like pores, flaws, and grain boundaries, by softening the elastic mineral frame, and consequently decreases the P- and S-wave velocities. Since we simulated an increase in pore pressure, we expect a decrease in velocity.

When a volume of fluid is injected into the reservoir, a pore pressure pulse is generated, whose time shape closely relates to the injection source time history. Its diffusive behavior depends on the characteristics of the medium and of the fluid itself. By using a simple rock physics approach, we modeled the diffusive propagation of the pore pressure perturbation from the injection point through the subsoil.

The two main equations used for the modeling are the Biot [11] equations of poroelasticity and the Eberhart-Phillips et al. [12] (hereafter referred to as EP) empirical relations between seismic velocity and pore pressure in sandstone rocks.

The linear dynamics of poroelastic deformation is described by the Biot [11] equations. In the general case, these equations imply the existence of two compressional waves (the P- and Biot wave) and one shear wave (the S-wave) in the system. The second compressional wave (Biot wave) is a diffusional one that corresponds to the process of pore pressure diffusion:

$$\frac{\partial p}{\partial t}=D{\nabla }^{2}p$$

(3)

where $p$ is the pore pressure and $D$ is the hydraulic diffusivity, which increases with the medium permeability and the fluid compressibility and decreases with the fluid viscosity.

Considering a concentrated source of overpressure, the analytic solution of Equation (3) in a homogeneous medium can be written as [26]

$$p\left(x,t\right)=\frac{A}{\sqrt{4D\pi t}} exp \left[-\frac{{\left(x-x_0\right)}^2}{4Dt}\right]$$

(4)

where $A$ is the initial pore pressure perturbation.

We calculated the initial pressure perturbation by fixing the injection rate and assuming that the fluid diffusivity is negligible with respect to that of the pressure [3].

Since we want to simulate cases which are as general as possible, we hypothesized a continuous injection up to the observation time. This assumption can also be found in other works of fluid injection simulation (see, for instance, [3,27]). If the injection rate history is known precisely, the proposed analysis can be easily adapted accordingly.

EP retrieved empirical relations between the effective pressure and the seismic velocities in sandstone from laboratory measurements. We used the equation related to the P-wave velocity model:

$$V_p=5.77-6.94\phi -1.73\sqrt{C}+0.446\left(P_e-e^{-16.7P_e}\right)$$

(5)

where $\phi$ is the porosity, $C$ is the clay content, i.e., the clay fraction in soils, and $P_e$ is the differential pressure. In a recent work of [28], the authors used these relations in order to link the velocity variations around the injection well to the effective pressure perturbation caused by the water injection.

In detail, starting from the P-wave velocity structure of Venisti et al. [14], we used the EP relations to retrieve the initial effective pressure field. We then modeled the propagation of the 3D pore pressure perturbation into the medium by considering a continuous injection up to the three observation times (i.e., after 1 day, 1 week, and 1 month from the start of the continuous water injections). The retrieved pore pressure perturbation field was added to the initial effective pressure field. Finally, we calculated the perturbed P-wave velocity model by using the EP relations.

By considering the well logs in the Adriatic area, available through the ViDEPI Project (at http://unmig.sviluppoeconomico.gov.it/videpi/pozzi/pozzi.asp), we modeled the target volume as mainly composed of sand banks with clay intercalations. Typical porosities of these formations are in the range of 0.2–0.5 (20–50%), depending on the compaction (i.e., depth) and the clay content [29]. The densities are in the range of 2200–2600 kg/m³ [30]. The shear modulus was set at 35 MPa [30], and the Poisson modulus is in the range 0.1–0.5 [30]. These model parameters are summarized in Table 1.

In order to explore a variety of scenarios, we used two values for the injection rate (see

Table 2. Details regarding the variable parameters in the rock modeling.

PP Perturbation Model		Phase Catalogue Construction		Observation Time
Diffusivity (m2/s)	Injection Rate (m3/day)	Noise Level (m/s)	Sea Bottom Interface Characteristics	After 1 Day
0.1 very low	500 low	10-8 Min. (MIN)	0.54 Soft (SF)	After 1 Week
1 high	1500 moderate	10-7 Max. (MAX)	0.26 Hard (HD)	After 1 Month

), one of 500 m³/day, typical of wastewater disposal operations in Italy [31], and the other of 1500 m³/day, considered as a moderate value for a wastewater disposal operation [3]. We also used two different hydraulic diffusivity values (see Table 2), one representing a low (0.1 m²/s) and the other a high (1 m²/s) diffusivity medium [3]. In Table 2, we summarize the parameters used to build 12 different perturbed velocity models, i.e., four combinations of the model parameters (i.e., injection rate and diffusivity) for three observation times. This combination of two injection rate values and two diffusivity values allowed us to construct four different velocity models perturbed by fluid injection. The temporal evolution of the pore pressure perturbation front, in terms of extension and intensity, was strongly related to the chosen injection rate and diffusivity values. Indeed, under the same injection rate, a higher diffusivity implies a more broadly extended and less intense velocity anomaly. With the same diffusivity, a higher rate implies, instead, a more intense velocity anomaly. As an example, we show in Figure 4 the pore pressure perturbation (Figure 4a–c) and the consequent perturbed velocity model (Figure 4d–f) around the injection well (located in the center), in the case of the rate of 500 m³/day and diffusivity of 1 m²/s.

2.3. D Tomographic Technique

The complete catalogue of P-phase arrival times was created by considering wave propagation in the perturbed 3D velocity models. Since in real data acquisition the waves generated by a shot may not be recorded at all stations due to attenuation phenomena and noise level at the recording site, we selected the travel-time dataset for the analysis by considering the attenuation law for the area under study (see Section 2.1). Moreover, we added an uncertainty to the P-wave selected arrival times and we assigned an inversion weight based on the estimated S/N ratio, according to the conversion rates reported in

Table 3. Conversion table from the signal to noise ratio to the measurement uncertainties on the arrival times and inversion weights. The symbol “/” indicates that the correspondent data are not considered.

S/N Ratio	Error on Arrival Time (s)	Inversion Weight
>50	0.01	1
20–50	0.05	0.8
10–20	0.1	0.6
5–10	0.2	0.3
2–5	0.5	0.1
<2	/	/

. The uncertainties are randomly extracted from a uniform distribution of interval length equal to two times the corresponding S/N ratio.

In order to cover different scenarios, this procedure has been performed for the four different combinations of noise level and soil characteristics, allowing the generation of four datasets (i.e., SF MIN, SF MAX, HD MIN, HD MAX) for each epoch. Figures of uncertainties and weights distributions for each catalogue can be found in S2 of SM.

The 4D tomography was performed by applying the 3D tomography to the datasets built for the three observation times for each combination of noise and soil characteristics.

To perform the 3D seismic tomography, we used a method based on accurate finite difference travel-time computation and a simultaneous inversion of body wave travel times to infer velocity and source location parameters [32]. This method is well proven and used in many applications in tectonic regions and volcanic and geothermal areas using data from passive (e.g., [33,34,35]) and active surveys (e.g., [36,37]). In the case of datasets from active surveys, the source location parameters are not inverted and fixed to the known positions. The inversion approach is based on an iterative scheme, in which the linearized delay time inversion is performed. First arrival travel times are computed through a finite difference solution of the eikonal equation [38] in a fine grid of 0.5 × 0.5 × 0.5 km³. Travel times for each event–receiver pair are recalculated by numerical integration of the slowness on the inversion grid field along the rays traced in the finite difference travel-time field. The parameters are inverted using the least squares root (LSQR) method of Paige and Sanders [39]. The velocity model is parametrized by a nodal representation defined by a tridimensional grid. The spacing of the nodes can be different in the three dimensions. The optimal grid mesh has been chosen taking into account the dimension of the investigated volume, the source–station geometry, and the shape of heterogeneities, and corresponds to 1 × 1 × 1 km³.

Since we used a linearized approach for the inversion, the choice of the initial velocity model is crucial. In order to explore and quantify the issue of the dependence of the final solution on the starting velocity model, we generated 200 random initial velocity models (avoiding velocity inversion) in a range of 1 km around the reference (Figure 5a), then we inverted the noised dataset starting from each one of these initial models. The considered scenario is the case with the rate of 500 m³/day and diffusivity of 1 m²/s, after one month of continuous injection using the SF MIN catalogue.

Figure 5b shows that even though different initial models provided substantially different rms time residuals, the data misfit of all final solutions reached about the same value (close to ~0.018 s). We obtained a similar result using the reference initial model (red star in Figure 5b).

Selecting only the tomographic results for the models with a final rms less than 0.0175 s, we estimated the standard deviation of the velocity distribution at each node of the velocity model. The representation of the results and the standard deviation at the nodes affected by the anomaly provide an estimation of the sensibility of the final result to the initial velocity model (Figure 5c).

3. Results

The target of the synthetic tomographic study is the accurate reconstruction of velocity anomalies relative to the background model, their extent over time, and the location of the anomaly peak. The reliable retrieval of these features allows us to monitor the condition of the volume hosting the reservoir and can be used to manage the injection operations. We analyzed the tomographic results in terms of the parameters defined above.

In Figure 6, we show an example of retrieved tomographic models in the case of the rate of 500 m³/day and diffusivity of 0.1 m²/s, after one month of continuous injection. In the horizontal slices, the white regions represent the areas not covered by seismic rays. Moreover, in Figure 6, we show a comparison between the tomographic models retrieved with two different catalogues (i.e., SF MIN and HD MAX). The different intensities of the retrieved anomalies (see the light blue sphere at 2, 3, and 4 km of depth in the center of the model) in the two tomographic images are strongly correlated with the dataset quality and seabed characteristics. However, the shape and the location of the anomaly are well retrieved in both cases.

In order to quantitatively evaluate the tomography performances for all the considered cases (see Table 1), we analyzed the retrieved models by extracting the 1D anomaly trend (i.e., the difference between the retrieved models and the background model), centered at the injection point, along the three spatial directions (i.e., NS, EW, and Z). In Figure 7, we report these 1D anomaly trends (gray lines) in comparison with the true 1D anomaly trends (red lines, i.e., the difference between the true perturbed model and the background model).

The panels in Figure 7a are related to the model perturbed by an injection with a rate of 500 m³/day in a medium with diffusivity of 0.1 m²/s. The tomographic images reproduce the anomaly intensity over time well in the cases of the minimum noise level (solid gray lines). By considering the cases with the maximum noise level (dashed gray lines), the anomaly location and extension are well reconstructed over time in the SF case (dashed dark gray lines), while the retrieved intensity is about 2/3 of the true value, not well reconstructed but still detectable with respect to the background model. In the HD case (dashed light gray), the anomaly features are not well reconstructed.

By considering the model perturbed by the same injection rate with a higher hydraulic diffusivity, i.e., 1 m²/s (Figure 7b), the anomaly is well reconstructed in the cases of the minimum noise level (solid gray lines), and the anomaly location and extent are well defined in the cases of the maximum noise level (dashed gray lines). Moreover, in these cases, after one week of continuous injection, the intensity of the retrieved anomaly is about half of the true value. After one month, the retrieved intensity is about 2/3 of the true value. In these cases, the anomaly can be considered well reproduced even with the maximum noise level.

The panels in Figure 7c are related to the model perturbed by an injection rate of 1500 m³/day in a medium with diffusivity of 0.1 m²/s. The tomographic images reconstruct the anomaly intensity and location over time well in the cases of the minimum noise level (solid gray lines), while the anomaly extension is overestimated after one week of continuous injection. After one month of continuous injection, the intensity of the retrieved anomaly is about ¾ of the true value, so well reconstructed. By considering the cases with the maximum noise level (dashed gray lines), the anomaly location is well determined over time, while the retrieved intensity is about ½ of the true value, not well reconstructed but still detectable with respect to the background model.

By considering the model perturbed by the same injection rate with a higher hydraulic diffusivity, i.e., 1 m²/s (Figure 7d), the anomaly is well reconstructed in the cases of the minimum noise level (solid gray lines) after one day and one week, and the anomaly location and extension are well defined for all the other cases. In the cases of the minimum noise level, after one month, the reconstructed anomaly intensity is about half of the true value, so not well reconstructed. By considering the cases with the maximum noise level (dashed gray lines), the anomaly location and extension are well reconstructed over time, while the retrieved intensity is about half of the true value, not well reconstructed but still detectable with respect to the background model.

A schematic summary of tomographic study capability in reconstructing at least 2/3 of the anomaly intensity can be found in

Table 4. Tomographic method capability to reconstruct at least 2/3 of the true velocity anomaly intensity, for each combination of simulation parameters.

	Rate = 500 m3/Day								Rate = 1500 m3/Day
	D = 0.1 m2/s				D = 1 m2/s				D = 0.1 m2/s				D = 1 m2/s
	SF MIN	SF MAX	HD MIN	HDMAX	SF MIN	SF MAX	HD MIN	HD MAX	SF MIN	SF MAX	HD MIN	HD MAX	SF MIN	SF MAX	HD MIN	HD MAX
Day	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘
Week	✔	✔	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘
Month	✔	✔	✔	✘	✔	✘	✔	✘	✔	✘	✔	✘	✘	✘	✘	✘

.

4. Discussion

Our analysis relies on the following aspects: (1) we modeled the physical process of pore pressure diffusion in the investigated volume and the consequent perturbation of the velocity model; (2) in the active survey set-up, we accounted for an ad hoc attenuation law in order to build realistic arrival-time catalogues; (3) we simulated different scenarios by considering different parameters of the medium (crustal structure, hydraulic diffusivity, and noise level at stations) and injection (injection rate) typical of common cases in offshore exploitation activities; and (4) we evaluated the results in terms of the retrieved features of the velocity anomaly.

We fixed the source–receiver configuration in order to properly resolve the volume of a few kilometers embedding the injection reservoir. Figure 4 shows that this volume of interest, up to 6 km of depth, is well sampled by seismic rays, confirming the validity of the chosen setting. Therefore, the optimized design of the source–station geometry of the active surveys is crucial to properly image the volume of interest [40]. This is only possible through accurate ray paths reconstruction in the volume and, therefore, through reliable knowledge of the crustal velocity model.

The modeled velocity changes, induced by the fluid injection, range between 5% and 20%, depending on the medium and injection conditions. These values are in agreement with other works on the 4D seismic tomography application in exploitation areas [7,8,9,41].

We found that, in the proper parameter condition, a conventional tomographic technique, such as the one that inverts first P arrival times, is sufficient for accurate tracking of the pore pressure front expansion using the simulated dense active source–station geometry.

From the results we infer that the noise level at the OBS stations is a crucial parameter in this kind of analysis. In fact, the anomaly is completely identified and characterized with the minimum noise level in all simulated cases. In the case of the maximum noise level, although the location and extension of the anomaly are well retrieved, the failure of the accurate identification of the anomaly intensity can lead to a misinterpretation of the reservoir conditions. As an example, referring to the case with an injection rate of 500 m³/day and diffusivity of 0.1 m²/s, the real anomaly (red lines in Figure 7c) grows in intensity (i.e., about 1.7 km/s after one day, 1.9 km/s after one week, and 2.5 km/s after one month), while the retrieved readings in the HD case (dashed light gray lines in Figure 7c) slightly increase from 0.5, after one day, to 0.7 km/s, after one month. This observation would induce the misinterpretation that the pore pressure perturbation does not cause an instability that propagates into the injection volume and, consequently, that the reservoir state is stable in contrast to the real expansion of the perturbation. Therefore, the choice of installing sea-bottom borehole sensors, or deploying small-aperture arrays in order to increase the S/N ratio, could be relevant for this kind of monitoring technique.

However, there are cases in which the expansion of the anomaly is well reconstructed even with maximum noise level, i.e., both cases with diffusivity of 1 m²/s after one week (gray dashed line in Figure 7b–d, central panels) and one month (gray dashed lines in Figure 7b–d, left panels). Moreover, in the case of a rate of 1500 m³/day and diffusivity of 1 m²/s after one month (SF and HD, gray dashed lines in Figure 7d, left panels), the anomaly intensity is not well recovered even with the minimum noise level.

As to the different soil types, i.e., SF and HD, the differences in terms of the retrieved intensity of the anomaly are relevant in the cases of the maximum noise level (dashed dark and light gray lines in Figure 7). The retrieved intensity is higher in the case of soft soil (SF), from 10% to 30%, with respect to hard soil (HD).

Through a bootstrap analysis, we quantified that the initial 1D velocity model could affect the amplitude of the retrieved final velocity anomaly with a deviation between 10% and 25% (Figure 5c).

The discussed results lead us to affirm that, beside the noise level, the other simulation parameters characterizing the propagation medium and the injection operation are also crucial in the reliable evaluation of this feasibility study. As to, in particular, the host medium, information about velocity and attenuation parameters is crucial in order to build reliable attenuation curves, which, in turn, affect the catalogue quality. In areas affected by industrial operations, the crustal structure of the exploited volume is well characterized, so the elastic/anelastic parameters of the unperturbed medium are available.

The use of air gun sources in marine geophysical surveys raises the concern of public opinion worldwide (see, for instance, the critical review in [42]). The main criticism of this technique is associated with the impact on the marine ecosystem. Nevertheless, in recent years, several air gun surveys have been performed by oil and gas companies and research institutes for exploration and/or research purposes. In the present simulation, we fixed the number and volumes of air gun sources according to the characteristics of prospecting granted in recent years in order to set up a realistic, feasible experiment.

Usually, the space-time evolution of induced seismicity is modeled in order to recognize the pore pressure front propagation into the medium embedding the reservoir (see, for instance, [43,44]). Our methodology has the advantage that the state of the volume containing the reservoir can be monitored even before the triggering of induced seismicity. In addition, the use of active seismic data ensures the same station–source layout in the different time epochs and so the same resolution of the 4D tomographic images. This is a critical condition for the application of time-lapse seismic tomography [45].

Another interesting point of the presented methodology is that the availability of 3D smooth velocity models continuously updated is important in locating with extreme accuracy and characterizing, in terms of source parameters, the natural/induced seismicity with the aim of monitoring the time-dependent hazard in industrial exploitation areas. Indeed, for seismic events near fluid injection wells, knowing their accurate location and the pressurized volume could be crucial to discriminating if an event was injection-induced or not, a question of high social importance that is inherently difficult to answer [1,2].

Although our study is based on what is dictated by the Italian guidelines for the monitoring of seismicity, ground deformation, and pore pressure in oil–gas exploitation areas where fluid injection is operated, we believe that the proposed strategy for detecting and tracking medium velocity changes can be of general interest in mitigating the seismic hazard of injection-induced earthquakes which can possibly be caused by industrial activities. Indeed, in the performed simulations, we considered medium and injection parameters that cover most common cases of offshore and onshore exploitation activities. Therefore, our results can be used as first-order operative indications to properly design active surveys in order to use the tomographic investigation as a monitoring tool of the volume containing the reservoir.

Whenever the operational and environmental conditions are found to be different from those of the performed simulation, this work provides for a complete methodology, which includes the physical modeling of pore pressure perturbation propagation and the reliable evaluation of the attenuation law in offshore cases, aimed at the evaluation of the feasibility of reservoir condition monitoring with seismic time-lapse tomography and consequent active survey projecting.

5. Conclusions

Our results indicate that 4D travel-time velocity tomography based on active surveys can be suitable for tracking the pore pressure front diffusion in the volume hosting the reservoir under the condition that:

• The source–receiver configuration of the active surveys is properly designed through considering the volume to be resolved (the few kilometers embedding the reservoir), the operational (injection rate) and environmental (crustal structure and medium hydraulic diffusivity) conditions, and the noise level at the stations;
• The noise level at the station is comparable to or lower than 10⁻⁸ m/s in cases of both a low (500 m³/day) and moderate (1500 m³/day) injection rate and low (0.1 m²/s) and high (1 m²/s) diffusivity. Alternatively, the noise level at the station is about 10⁻⁶ m/s in the case of high diffusivity (1 m²/s).

Moreover, we believe that:

• Despite our case study is about an area of the Italian offshore, the proposed methodology, which allows for the tracking of medium velocity changes even before the occurrence of induced earthquakes, can be of general interest to the mitigation of seismic hazards in exploited areas;
• Since in the simulated scenarios we considered medium and injection parameters that cover most common cases of exploitation operations, our result can be used as first-order operative indications to properly design active surveys in order to use the tomographic investigation as a monitoring tool of reservoir condition.

In conclusion, this work provides a complete methodology, which includes the physical modeling of pore pressure perturbation propagation and the reliable evaluation of the attenuation law in offshore cases, aimed at the evaluation of the feasibility of reservoir condition monitoring with seismic time-lapse tomography and consequent active survey projecting.