Exploration Geophysics - Volume 37, Issue 1, 2006

We measured the changes in P-wave velocity that occur when injecting CO2 in gaseous, liquid, and supercritical phases into water-saturated anisotropic sandstones. P-wave velocities were measured in two cylindrical samples of Tako Sandstone, drilled along directions normal and parallel to the bedding plane, using a piezo-electric transducer array system. The velocity changes caused by CO2 injection are typically -6% on average, with maximum values about -16% for the case of supercritical CO2 injection. P-wave velocity tomograms obtained by the differential arrival-time method clearly show that CO2 migration behaviour is more complex when CO2 flows normal to the bedding plane than when it flows parallel to bedding. We also found that the differences in P-wave velocity images were associated both with the CO2 phases and with heterogeneity of pore distribution in the rocks. Seismic images showed that the highest velocity reduction occurred for supercritical CO2 injection, compared with gaseous or liquid CO2 injection. This result may justify the use of the seismic method for CO2 monitoring in geological sequestration.

The first Japanese pilot-scale CO2 sequestration project has been undertaken in an onshore saline aquifer, near Nagaoka in Niigata prefecture, and time-lapse well logs were carried out in observation wells to detect the arrival of injected CO2 and to evaluate CO2 saturation in the reservoir. CO2 was injected into a thin permeable zone at the depth of 1110 m at a rate of 20–40 tonnes per day. The total amount of injected CO2 was 10 400 tonnes, during the injection period from July 2003 to January 2005. The pilot-scale demonstration allowed an improved understanding of the CO2 movement in a porous sandstone reservoir, by conducting time-lapse geophysical well logs at three observation wells. Comparison between neutron well logging before and after the insertion of fibreglass casing in observation well OB-2 showed good agreement within the target formation, and the higher concentration of shale volume in the reservoir results in a bigger difference between the two well logging results. CO2 breakthrough was identified by induction, sonic, and neutron logs. By sonic logging, we confirmed P-wave velocity reduction that agreed fairly well with a laboratory measurement on drilled core samples from the Nagaoka site. We successfully matched the history changes of sonic P-wave velocity and estimated CO2 saturation after breakthrough in two observation wells out of three. The sonic-velocity history matching result suggested that the sweep efficiency was about 40%. Small effects of CO2 saturation on resistivity resulted in small changes in induction logs when the reservoir was partially saturated. We also found that CO2 saturation in the CO2-bearing zone responded to suspension of CO2 injection.

Japan’s first pilot-scale CO2 sequestration experiment has been conducted in Nagaoka, where 10 400 t of CO2 have been injected in an onshore aquifer at a depth of about 1100 m. Among various measurements conducted at the site for monitoring the injected CO2, we conducted time-lapse crosswell seismic tomography between two observation wells to determine the distribution of CO2 in the aquifer by the change of P-wave velocities. This paper reports the results of the crosswell seismic tomography conducted at the site.

The crosswell seismic tomography measurements were carried out three times; once before the injection as a baseline survey, and twice during the injection as monitoring surveys. The velocity tomograms resulting from the monitoring surveys were compared to the baseline survey tomogram, and velocity difference tomograms were generated. The velocity difference tomograms showed that velocity had decreased in a part of the aquifer around the injection well, where the injected CO2 was supposed to be distributed. We also found that the area in which velocity had decreased was expanding in the formation up-dip direction, as increasing amounts of CO2 were injected. The maximum velocity reductions observed were 3.0% after 3200 t of CO2 had been injected, and 3.5% after injection of 6200 t of CO2.

Although seismic tomography could map the area of velocity decrease due to CO2 injection, we observed some contradictions with the results of time-lapse sonic logging, and with the geological condition of the cap rock. To investigate these contradictions, we conducted numerical experiments simulating the test site. As a result, we found that part of the velocity distribution displayed in the tomograms was affected by artefacts or ghosts caused by the source-receiver geometry for the crosswell tomography in this particular site. The maximum velocity decrease obtained by tomography (3.5%) was much smaller than that observed by sonic logging (more than 20%). The numerical experiment results showed that only 5.5% velocity reduction might be observed, although the model was given a 20% velocity reduction zone. Judging from this result, the actual velocity reduction can be more than 3.5%, the value we obtained from the field data reconstruction. Further studies are needed to obtain more accurate velocity values that are comparable to those obtained by sonic logging.

In 2006, the Cooperative Research Centre for Greenhouse Gas Technologies (CO2CRC) plans to undertake (subject to receiving the necessary approvals) a Pilot program for CO2 storage within a depleted gas reservoir. The Otway Basin Pilot Program (OBPP) aims to demonstrate that subsurface CO2 storage is both economically and environmentally sustainable in Australia. This will be the first CO2 storage program in the world to utilise a depleted gas reservoir and, hence, the experience gained will be a valuable addition to the range of international CO2 storage programs that are underway or being planned.

A key component of the OBPP is the design of an appropriate geophysical monitoring strategy that will allow the subsurface migration of the CO2 plume to be tracked and to verify that containment has been successful. This paper presents the results from modelling the predicted gravity response to CO2 injection into the Otway Basin reservoir, where the goal was to determine minimum volumes of CO2 that may be detectable using non-seismic geophysical techniques. Modelling results indicate that gravity measurements at 10 m spacing within the existing observation well and the planned CO2 injection well would provide excellent vertical resolution, even for the smallest CO2 volume modelled (10 000 tonnes), but resolving the lateral extent of the plume would not be possible without additional wells at closer spacing.

In this study, the resolution capabilities of electrical resistivity tomography (ERT) in the monitoring of CO2 injection are investigated. The pole-pole and bipole-bipole electrode configuration types are used between two uncased boreholes straddling the CO2 plume. Forward responses for an initial pre-injection model and three models for subsequent stages of CO2 injection are calculated for the two different electrode configuration types, noise is added and the theoretical data are inverted with both L1- and L2-norm optimisation.

The results show that CO2 volumes over a certain threshold can be detected with confidence. The L1-norm proved superior to the L2-norm in most instances. Normalisation of the inverted models with the pre-injection inverse model gives good images of the regions of changing resistivity, and an integrated measure of the total change in resistivity proves to be a valid measure of the total injected volume.

The risk of fault reactivation in the Gippsland Basin was calculated using the FAST (Fault Analysis Seal Technology) technique, which determines fault reactivation risk by estimating the increase in pore pressure required to cause reactivation within the present-day stress field. The stress regime in the Gippsland Basin is on the boundary between strike-slip and reverse faulting: maximum horizontal stress (~ 40.5 MPa/km) > vertical stress (21 MPa/km) ~ minimum horizontal stress (20 MPa/km). Pore pressure is hydrostatic above the Campanian Volcanics of the Golden Beach Subgroup. The NW-SE maximum horizontal stress orientation (139°N) determined herein is broadly consistent with previous estimates, and verifies a NW-SE maximum horizontal stress orientation in the Gippsland Basin. Fault reactivation risk in the Gippsland Basin was calculated using two fault strength scenarios; cohesionless faults (C = 0; μ = 0.65) and healed faults (C = 5.4; μ = 0.78). The orientations of faults with relatively high and relatively low reactivation potential are almost identical for healed and cohesionless fault strength scenarios. High-angle faults striking NE-SW are unlikely to reactivate in the current stress regime. High-angle faults oriented SSE-NNW and ENE-WSW have the highest fault reactivation risk. Additionally, low-angle faults (thrust faults) striking NE-SW have a relatively high risk of reactivation. The highest reactivation risk for optimally oriented faults corresponds to an estimated pore pressure increase (Delta- P) of 3.8 MPa (~548 psi) for cohesionless faults and 15.6 MPa (~2262 psi) for healed faults. The absolute values of pore pressure increase obtained from fault reactivation analysis presented in this paper are subject to large errors because of uncertainties in the geomechanical model (in situ stress and rock strength data). In particular, the maximum horizontal stress magnitude and fault strength data are poorly constrained. Therefore, fault reactivation analysis cannot be used to directly measure the maximum allowable pore pressure increase within a reservoir. We argue that fault reactivation analysis of this type can only be used for assessing the relative risk of fault reactivation and not to determine the maximum allowable pore pressure increase a fault can withstand prior to reactivation.

Three differing sandstones, two synthetic and one field sample, have been tested ultrasonically under a range of confining pressures and pore pressures representative of in-situ reservoir pressures. These sandstones include: a synthetic sandstone with calcite intergranular cement produced using the CSIRO Calcite In-situ Precipitation Process (CIPS); a synthetic sandstone with silica intergranular cement; and a core sample from the Otway Basin Waarre Formation, Boggy Creek 1 well, from the target lithology for a trial CO2 pilot project. Initial testing was carried on the cores at “room-dried” conditions, with confining pressures up to 65 MPa in steps of 5 MPa. All cores were then flooded with CO2, initially in the gas phase at 6 MPa, 22°C, then with liquid-phase CO2 at a temperature of 22°C and pressures from 7 MPa to 17 MPa in steps of 5 MPa. Confining pressures varied from 10 MPa to 65 MPa. Ultrasonic waveforms for both P- and S-waves were recorded at each effective pressure increment. Velocity versus effective pressure responses were calculated from the experimental data for both P- and S-waves. Attenuations (1/Qp) were calculated from the waveform data using spectral ratio methods. Theoretical calculations of velocity as a function of effective pressure for each sandstone were made using the CO2 pressure-density and CO2 bulk modulus-pressure phase diagrams and Gassmann effective medium theory.

Flooding the cores with gaseous phase CO2 produced negligible change in velocity-effective stress relationships compared to the dry state (air saturated). Flooding with liquid-phase CO2 at various pore pressures lowered velocities by approximately 8% on average compared to the air-saturated state. Attenuations increased with liquid-phase CO2 flooding compared to the air-saturated case. Experimental data agreed with the Gassmann calculations at high effective pressures. The “critical” effective pressure, at which agreement with theory occurred, varied with sandstone type. Discrepancies are thought to be due to differing micro-crack populations in the microstructure of each sandstone type. The agreement with theory at high effective pressures is significant and gives some confidence in predicting seismic behaviour under field conditions when CO2 is injected.

Injection of CO2 into underground saline formations, due to their large storage capacity, is probably the most promising approach for the reduction of CO2 emissions into the atmosphere. CO2 storage must be carefully planned and monitored to ensure that the CO2 is safely retained in the formation for periods of at least thousands of years. Seismic methods, particularly for offshore reservoirs, are the primary tool for monitoring the injection process and distribution of CO2 in the reservoir over time provided that reservoir properties are favourable. Seismic methods are equally essential for the characterisation of a potential trap, determining the reservoir properties, and estimating its capacity. Hence, an assessment of the change in seismic response to CO2 storage needs to be carried out at a very early stage. This must be revisited at later stages, to assess potential changes in seismic response arising from changes in fluid properties or mineral composition that may arise from chemical interactions between the host rock and the CO2. Thus, carefully structured modelling of the seismic response changes caused by injection of CO2 into a reservoir over time helps in the design of a long-term monitoring program. For that purpose we have developed a Graphical User Interface (GUI) driven rock physics simulator, designed to model both short and long-term 4D seismic responses to injected CO2. The application incorporates CO2 phase changes, local pressure and temperature changes, chemical reactions and mineral precipitation. By incorporating anisotropic Gassmann equations into the simulator, the seismic response of faults and fractures reactivated by CO2 can also be predicted.

We show field examples (potential CO2 sequestration sites offshore and onshore) where we have tested our rock physics simulator. 4D seismic responses are modelled to help design the monitoring program.

Array observations of the vertical component of microtremors are frequently conducted to estimate a subsurface layered-earth structure on the assumption that microtremors consist predominantly of the fundamental mode Rayleigh waves. As a useful tool in the data collection, processing and analysis, the spatial autocorrelation (SPAC) method is widely used, which in practice requires a circle array consisting of M circumferential stations and one centre station (called “M-station circle array”, where M is the number of stations). The present paper considers the minimum number of stations required for a circle array for efficient data collection in terms of analytical efficacy and field effort.

This study first rearranges the theoretical background of the SPAC algorithm, in which the SPAC coefficient for a circle array with M infinite is solely expressed as the Bessel function, J0(rk) (r is the radius and k the wavenumber). Secondly, the SPAC coefficient including error terms independent of the microtremor energy field for an M-station circle array is analytically derived within a constraint for the wave direction across the array, and is numerically evaluated in respect of these error terms. The main results of the evaluation are: 1) that the 3-station circle array when compared with other 4-, 5-, and 9-station arrays is the most efficient and favourable for observation of microtremors if the SPAC coefficients are used up to a frequency at which the coefficient takes the first minimum value, and 2) that the Nyquist wavenumber is the most influential factor that determines the upper limit of the frequency range up to which the valid SPAC coefficient can be estimated.

The full-waveform inversion algorithm using normalised seismic wavefields can avoid potential inversion errors due to source estimation required in conventional full-waveform inversion methods. In this paper, we have modified the inversion scheme to install a weighted smoothness constraint for better resolution, and to implement a staged approach using normalised wavefields in order of increasing frequency instead of inverting all frequency components simultaneously. The newly developed scheme is verified by using a simple two-dimensional fault model. One of the most significant improvements is based on introducing weights in model parameters, which can be derived from integrated sensitivities. The model-parameter weighting matrix is effective in selectively relaxing the smoothness constraint and in reducing artefacts in the reconstructed image. Simultaneous multiplefrequency inversion can almost be replicated by multiple singlefrequency inversions. In particular, consecutively ordered singlefrequency inversion, in which lower frequencies are used first, is useful for computation efficiency.

In this study, we propose a joint inversion method, using genetic algorithms, to determine the shear-wave velocity structure of deep sedimentary layers from receiver functions and surfacewave phase velocity. Numerical experiments with synthetic data indicate that the proposed method can avoid the trade-off between shear-wave velocity and thickness that arises when inverting the receiver function only, and the uncertainty in deep structure from surface-wave phase velocity inversion alone. We apply the method to receiver functions obtained from earthquake records with epicentral distances of about 100 km, and Rayleigh-wave phase velocities obtained from a microtremor array survey in the Kanto Plain, Japan. The estimated subsurface structure is in good agreement with the previous results of seismic refraction surveys and deep borehole data.

Recently, rapid developments in computer hardware have enabled reverse-time migration to be applied to various production imaging problems. As a wave-equation technique using the twoway wave equation, reverse-time migration can handle not only multi-path arrivals but also steep dips and overturned reflections. However, reverse-time migration causes unwanted artefacts, which arise from the two-way characteristics of the hyperbolic wave equation. Zero-lag cross correlation with diving waves, head waves and back-scattered waves result in spurious artefacts. These strong artefacts have the common feature that the correlating forward and backward wavefields propagate in almost the opposite direction to each other at each correlation point. This is because the ray paths of the forward and backward wavefields are almost identical. In this paper, we present several tactics to avoid artefacts in shot-domain reverse-time migration. Simple muting of a shot gather before migration, or wavefront migration which performs correlation only within a time window following first arriving travel times, are useful in suppressing artefacts. Calculating the wave propagation direction from the Poynting vector gives rise to a new imaging condition, which can eliminate strong artefacts and can produce common image gathers in the reflection angle domain.

Urban conditions, such as existing underground facilities and ambient noise due to cultural activity, restrict the general application of conventional geophysical techniques. At a tunnelling site in an urban area along an existing railroad, we used the refraction microtremor (REMI) technique (Louie, 2001) as an alternative way to get geotechnical information. The REMI method uses ambient noise recorded by standard refraction equipment and a linear geophone array to derive a shear-wave velocity profile. In the inversion procedure, the Rayleigh wave dispersion curve is picked from a wavefield transformation, and iteratively modelled to get the S-wave velocity structure.

The REMI survey was carried out along the line of the planned railway tunnel. At this site vibrations from trains and cars provided strong seismic sources that allowed REMI to be very effective. The objective of the survey was to evaluate the rock mass rating (RMR), using shear-wave velocity information from REMI. First, the relation between uniaxial compressive strength, which is a component of the RMR, and shear-wave velocity from laboratory tests was studied to learn whether shear-wave velocity and RMR are closely related. Then Suspension PS (SPS) logging was performed in selected boreholes along the profile, in order to draw out the quantitative relation between the shear-wave velocity from SPS logging and the RMR determined from inspection of core from the same boreholes. In these tests, shear-wave velocity showed fairly good correlation with RMR. A good relation between shear-wave velocity from REMI and RMR could be obtained, so it is possible to estimate the RMR of the entire profile for use in design of the underground tunnel.

The main purpose of this study is to apply spatial analysis using semivariograms to small-loop electromagnetic survey data to assess the extent of seawater intrusion in an experimental watershed. To indicate the extent of seawater intrusion over the study area, vertical electrical soundings at 33 points and electrical conductivity logging in two wells were conducted. From the correlation between resistivities obtained by inversion and the depth of the aquifer at the two wells, the region of seawater intrusion was identified and demonstrated by electrical conductivity logging results obtained over two years. To measure the variation of apparent conductivity with depth, an electromagnetic survey in six frequency bands was adopted. Apparent conductivity mapping with spatial analysis using semivariograms is an effective technique for identifying the region of seawater intrusion at shallow depth.

Magnetic survey flights by helicopters are usually parallel to the topographic surface, with a nominal clearance, but especially in high-resolution surveys the altitudes at which observations are made may be too variable to be regarded as a smooth surface. We have developed a reduction procedure for such data using the method of equivalent sources, where surrounding sources are included to control edge effects, and data from points distributed randomly in three dimensions are directly modelled. Although the problem is generally underdetermined, the method of conjugate gradients can be used to find a minimum-norm solution. There is freedom to select the harmonic function that relates the magnetic anomaly with the source. When the upward continuation function operator is selected, the equivalent source is the magnetic anomaly itself. If we select as source a distribution of magnetic dipoles in the direction of the ambient magnetic field, we can easily derive reduction-to-pole anomalies by rotating the direction of the magnetic dipoles to vertical.

Helicopter-borne aeromagnetic surveys at two different times separated by three years were conducted to better understand the shallow subsurface structure of the Vulcano and Lipari volcanic complex, Aeolian Islands, southern Italy, and also to monitor the volcanic activity of the area. As there was no meaningful difference between the two magnetic datasets to imply an apparent change of the volcanic activity, the datasets were merged to produce an aeromagnetic map with wider coverage than was given by a single dataset. Apparent magnetisation intensity mapping was applied to terrain-corrected magnetic anomalies, and showed local magnetisation highs in and around Fossa Cone, suggesting heterogeneity of the cone. Magnetic modelling was conducted for three of those magnetisation highs.

Each model implied the presence of concealed volcanic products overlain by pyroclastic rocks from the Fossa crater. The model for the Fossa crater area suggests a buried trachytic lava flow on the southern edge of the present crater. The magnetic model at Forgia Vecchia suggests that phreatic cones can be interpreted as resulting from a concealed eruptive centre, with thick latitic lavas that fill up Fossa Caldera. However, the distribution of lavas seems to be limited to a smaller area than was expected from drilling results. This can be explained partly by alteration of the lavas by intense hydrothermal activity, as seen at geothermal areas close to Porto Levante. The magnetic model at the north-eastern Fossa Cone implies that thick lavas accumulated as another eruption centre in the early stage of the activity of Fossa. Recent geoelectric surveys showed high-resistivity zones in the areas of the last two magnetic models.