Application of DoubleDifference Seismic Tomography to Carbon Sequestration Monitoring at the Aneth Oil Field, Utah,” to

Double difference seismic tomography was performed using travel time data from a carbon sequestration site at the Aneth oil field in southeast Utah as part of a Department of Energy initiative on monitoring, verification, and accounting (MVA) of sequestered CO2. A total of 1211 seismic events were recorded from a borehole array consisting of 23 geophones. Artificial velocity models were created to determine the likelihood of detecting a CO2 plume with an unfavorable event and receiver arrangement. In tests involving artificially modeled ray paths through a velocity model, ideal event and receiver arrangements clearly show velocity reductions. When incorporating the unfavorable event and station locations from the Aneth Unit into synthetic models, the ability to detect velocity reductions is greatly diminished. Using the actual, recorded travel times, the Aneth Unit results show differences between a synthetic baseline model and the travel times obtained in the field, but the differences do not clearly indicate a region of injected CO2. MVA accuracy and precision may be improved through the use of a receiver array that provides more comprehensive ray path coverage, and a more detailed baseline velocity model.


Introduction
Increasing concentrations of carbon dioxide (CO 2 ) in the atmosphere have become a significant international concern in recent decades, as they are believed to contribute to climate change.The majority of anthropogenic CO 2 emissions come from the burning of fossil fuels, which remain a considerable portion of the world's energy portfolio.One method of reducing anthropogenic CO 2 emissions is through geologic carbon sequestration.Geologic carbon sequestration is the storage of captured CO 2 in geologic formations, such as unmineable coal seams [1], depleted petroleum reservoirs [2], or deep saline aquifers [3].A major engineering challenge preventing the widespread adoption of geologic carbon sequestration is the ability to monitor, verify, and account (MVA) for the injected CO 2 .One tool with a potential application for performing these three tasks is double difference seismic tomography.
Many tools exist for monitoring sequestered CO 2 , including well testing, pressure monitoring, tracers and chemical sampling, surface and bore-hole seismics, and electromagnetic and geomechanical meters [4].These methods are not mutually exclusive and may be used in conjunction with one another to better characterize CO 2 activity underground, such as the monitoring system at In Salah which incorporates seismics, tiltmeters, tracers, and wireline logging, among other monitoring methods [5].The tomographic monitoring performed in this study differentiates itself from other common seismic monitoring systems, such as the vertical seismic profile and crosswell methods used to monitor CO 2 injection at Frio [6], in that it uses passive seismic events already present nearby or inside the reservoir to image CO 2 movements.
Generating an image, or tomogram, of an object by examining its reaction to the probing energy from an external source is the foundation of tomography [7].The external source, in the case of seismic tomography, is a seismic event.Fluid substitution within a rock mass can be imaged by measuring these seismic signals.The in situ rocks may be saturated with fluids, such as brine, which have a specific, inherent P wave velocity.A fluid being displaced by another fluid with different acoustic properties will affect the velocity of the medium.This velocity change, due to fluids with different elastic properties displacing each other, may be detected through tomography.This technique can be applied to monitoring injected CO 2 at carbon sequestration sites because CO 2 has different elastic properties than the in situ fluid [8].
Double difference seismic tomography [9] is a form of tomography in which it is assumed that similarly located seismic events will generate seismic waves that propagate along a similar path.Absolute and differential travel times are then used to create a three-dimensional velocity model through iterative reconstruction.The double difference equation is given as: where is the arrival time of the seismic wave, emanating from location i, arriving at station k, and is the arrival time of the seismic wave, emanating from location j, arriving at station k.The difference between the calculated and observed arrival time residual is , or the double difference [9].The double difference tomography software used in this application is tomoDD, which is based on the double difference location algorithm, hypoDD, written by Waldhauser [10].
The ability to use tomography in monitoring, verifying, and accounting for sequestered CO 2 in the Aneth oil field will be demonstrated through analysis of time lapse images.Event to receiver travel time data was provided along with known receiver locations and estimated event locations.In an attempt to verify the velocity changes associated with CO 2 injection, time lapse images will be compared to results from artificially created plumes.Combining these methods may allow for the determination of CO 2 plume extents, changes in CO 2 concentrations, as well as monitoring of CO 2 leakage from the reservoir.

Aneth Site Characteristics
The Aneth Unit is located in the Pennsylvanian Paradox Basin of southeastern Utah, and is part of the Greater Aneth oil field, shown in Figure 1.The Aneth Unit is part of an ongoing carbon sequestration operation, conducted in conjunction with the Southwest Regional Partnership on Carbon Sequestration (SWP).The SWP in comprised of state and federal agencies, universities, electric utilities, non governmental organizations, coal and gas companies, and the Navajo Nation [11].The target reservoir for carbon sequestration is the Desert Creek carbonate, which overlies the Chimney Rock shale and underlies the Gothic shale.The Gothic shale acts as the primary caprock for the reservoir.In addition to the Gothic shale, the Desert Creek reservoir is overlain by approximately 244 m thick layer of tight carbonates, evaporites, and fine-grained siliciclastics.Above that is a 975 m thick layer of shales, siltstones, sandstones, and evaporites, providing multiple layers of confinement against CO 2 migration [13].The Desert Creek reservoir has an average vertical extent of approximately 15 m and porosity of 10.3%.The permeability of the reservoir ranges from 6 to 27 millidarcies, and has a water saturation of 23.3% [14].The reservoir geometry can be seen in Figure 2, for the top of the Desert Creek reservoir.Injection of CO 2 into the Desert Creek reservoir first began in August 2007.The planned rate of injection into the reservoir was 136,000 t of CO 2 per year from approximately 16 wells across the 5 km 2 test site [11].Seismic events were monitored through a vertical array of 23 geophones.The geophone array extended from 805 to 1704 m in depth.
A total of 1211 passive seismic events were recorded between 25 April 2008 and 16 March 2009 from a monitor well located in the western part of the Aneth unit.The event location method is described by Rutledge and Soma [15].The recorded events were spread unevenly throughout this monitoring period, as shown in Figure 3. Moment scalars of the events range from 1.46 × 10 7 to 1.44 × 10 10 Nm.P wave arrival times as well as seismic wave origin times were provided for each event to receiver pair.As shown in Figure 4, the seismic events form both a northern cluster (approximately 4% of all events) and a southern cluster (the remaining 96% of events), with one event outlying the rest near the monitor well.These events are arranged in a way that may suggest the existence of two NW to SE striking fault zones.The estimated locations for the events in the northern cluster are at the approximate depth of the Desert Creek reservoir, with an average depth of 1761 m and standard deviation of 21 m, while the estimated locations for the events in the southern cluster are noticeably below it, with an average depth of 2007 m and standard deviation of 22 m.The error associated with these event locations can be difficult to estimate, as the velocity of the medium through which the seismic waves travel is not precisely known.Due to the distance from the sources to the receivers, the source depths are highly dependent on the velocity model and degree of uncertainty [15].Attempts to correlate the fluid volume changes (oil, water, or CO 2 ) with microseismicity were inconclusive [15].

Aneth Site Modeling
In order to perform a time lapse analysis of the provided event data, separate time periods had to be defined and the number of events needed to be refined.While two clusters of events were provided, the northern and southern, there are only 44 events in the northern cluster.As this represents a small fraction of the total number of events, it was decided that the decrease in resolution created by expanding the studied area to include this cluster would be less beneficial than a higher resolution model centered on the majority of the events.As a result, the northern cluster, as well as the event located near the monitor well, were removed to allow for creation of a higher resolution velocity model surrounding the southern cluster and monitor well.The remaining 1166 events, in the southern cluster, were separated into four time periods, as given in Table 1.The time periods were selected based on a balance between maintaining similarly spaced time periods and creating as many time periods as possible while maintaining enough events to allow for similar characterization of the reservoir across all time periods.The background velocity model used throughout this application of double difference tomography is given in Figure 5.It ranges in velocity from 4340 to 5270 m/s and is divided into five layers.This velocity model was used for creation of synthetic travel times to compare to real travel times.This background velocity model is assumed to be representative of the site geology prior to any CO 2 injection.During tomographic inversion, a node spacing of 80, 50 and 35 m in the x, y and z directions are used.As it is unknown which nodes will be traversed by the wave during tomography, the nodes which experience a velocity change will be assumed to have ray path coverage while those that do not change are removed, as they exceed the model boundaries.

Synthetic Experimental Plume
An unfavorable event and receiver arrangement will adversely affect the ability to precisely locate velocity changes in a rock mass.Low velocity zones being stretched, misplaced, or absent altogether in a synthetic test would suggest similar results may occur with the Aneth Unit data.It is necessary, therefore, to establish what degree of stretching or misplacement may be expected in the Aneth Unit results.To determine this, a travel time calculation code, using the fast marching ray tracing method, was used to calculate travel times from events to receivers for two synthetic tests.
As a first step towards determining the degree to which velocity changes can be accurately and precisely imaged through double difference tomography, synthetic tests are performed under both ideal conditions and those more closely matching the conditions of Aneth.To verify that detecting velocity changes in a layered velocity model under ideal conditions is possible, a scenario is created in which 1000 events are randomly distributed through a 1000 × 1000 × 1000 m velocity model with 8 receivers on the vertexes and 42 receivers distributed randomly throughout.This event and receiver arrangement ensures comprehensive ray path coverage of the entire velocity model.A cross section of the velocity model used for this test is shown in Figure 6a.The velocities were chosen to match the geology of the Aneth Unit.This synthetic test includes a simulation of travel times through the velocity model without a saturated layer, and another simulation of travel times with a saturated layer, from a depth 450 to 550 m.The velocity reduction magnitude due to CO 2 saturation as well as extent of affected area are unknown, so a simplified simulation involving a 10% (0.52 km/s) velocity reduction across the entire 450-550 m depth range is performed to test the method under controlled, idealized conditions.With the two different simulations, the results were compared to see if the area of change precisely matched the created area of change in the initial travel time model.The second scenario is similar to the first but more accurately models the conditions of the Aneth Unit site, shown in Figure 6b.Instead of using randomly dispersed events and ideal receiver locations, the 1166 events and 23 receiver locations from the Aneth Unit field data were used, in addition to the velocity model extents from the Aneth Unit site.This scenario also includes both an unsaturated and a saturated condition to model the effects of CO 2 saturation on the resultant tomograms.The layer of saturation corresponds to the extents of the Desert Creek reservoir, 1700 to 1783 m deep, and runs throughout the model extents in the lateral directions.The velocity model used for the initial tomographic inversion was set to a constant 5 km/s.The methods of tomographic inversion for the model simulating the Aneth conditions are the same as the methods described in Section 2.1.For the model simulating ideal conditions for tomography, the node spacing used for inversion is 25 m in the x, y and z directions.
To determine how reliable the results can be with this arrangement, several more synthetic tests are performed, using the checkerboard test, each with increasing similarity to the condition of the Aneth data.The following three synthetic tests contain velocities through a 6 by 6 checkerboard velocity model, with one set of checkers having a velocity of 4.5 km/s and the alternating checker having a velocity of 5 km/s.During tomographic inversion, each initial velocity model is held constant at 5 km/s, rather than starting the model with an accurate background.The first of these has a random distribution of 198 events in a cube used earlier, with 50 receivers distributed on the perimeter and throughout, ensuring comprehensive ray path coverage.The second test uses the Aneth receiver locations, but the number of events is limited to those in the first time period, 126, which is the smallest number of events used in a timer period, to determine if there is a significant difference between using all the events or a much smaller subset when they are all similarly located.The last test uses the Aneth event and receiver locations, limited to the 1166 events in the southern cluster.This test will highlight differences between using a small subset of the data.

Synthetic Analysis
Distinct low and high velocity zones are expected after performing tomographic inversion.These zones can make it difficult to detect the area that changes between the different simulations.For this time-lapse analysis, the unsaturated condition is subtracted from the saturated condition to reveal areas of change between the two conditions.Figure 7 shows the difference in the first test, with the randomly dispersed events.
The low velocity zone is clearly visible in the 450 to 550 m range, with no significant velocity anomalies appearing anywhere else in the velocity model.In addition to the low velocity zone being accurately located, the change in velocity closely matches the actual velocity reduction of 0.52 km/s.
The second scenario, which more closely models the Aneth Unit geology and geometry is analyzed in the same manner as the first scenario.The unsaturated condition is subtracted from the saturated condition to reveal areas of change between the two conditions, shown in Figure 8.The velocity change for this scenario does not match the induced velocity change, as was the case in the first scenario.A small area of decreasing velocity is found in the 1700 to 1783 m range, but does not stretch across the entire extents of the model.In addition to the smaller affected area, the magnitude of velocity change is smaller than the 0.52 km/s induced velocity change for the layer when the travel times were generated.The poor results obtained from the second synthetic test may be a result of the poor event and receiver arrangement.A checkerboard test is used to determine how well alternating velocities can be resolved under several controlled scenarios.The first test, using an ideal event and receiver arrangement inside of a cubic velocity model is shown in Figure 9.The 6 by 6 checkerboard is well resolved with this arrangement except for the outermost checkers where there is little raypath coverage.
The second test, using the stations and events from the first time period is shown in Figure 10.There is no clear checkerboard pattern that appears in this tomogram.There is an alternation between high and low velocity, but it does not resemble the square pattern actually found in the initial velocity model.The 126 events used in the second test is significantly fewer than the total number of events available, however, due to the clustered nature of the events, they all experience a similar travel path.To determine whether including more events would improve the tomogram, resolving the checkerboard pattern more accurately, the same test is performed instead using all 1166 events, shown in Figure 11.The inclusion of more events does not improve the accuracy of the tomogram.Faint layers of higher and lower velocity are present, but they do not resemble a checkerboard pattern.

Aneth Unit Control Comparisons
The Desert Creek reservoir lies on the upper boundary of the 5.22 km/s to 4.34 km/s velocity boundary in the background velocity model, at 1783 m deep.This velocity change is significantly more pronounced than the expected velocity change.In order to remove the effect of the layering on the results, each time period for the time-lapse analysis will be compared to a synthetic analysis, using artificial travel times generated through the velocity model provided in Figure 5, and using the same events from the Aneth Unit analysis.By subtracting the Aneth tomograms from the synthetic tomograms, using the same event and receiver locations, anomalies resulting from the tomography should be reduced, leaving low and high velocity regions unaccounted for in the original geologic model.Both velocity models used for the initial tomographic inversion were set to a constant 5 km/s background velocity.
Each time period, compared to its synthetic counterpart, is shown in Figure 12.The tomograms are positioned at the midpoint between the stations and the centroid of the event cluster.Due to the arrangement of the events and receivers, significant stretching or misplacement of the velocity changes is expected.
The tomograms tend to show a change towards a higher velocity with decreasing depth, and a lower velocity with increasing depth.The low velocity regions are below the expected saturation layer, in the Desert Creek reservoir, but as shown in the synthetic tests, a low velocity zone is unlikely to appear in that location, even if fully saturated.

Discussion
Unlike with the first two synthetic tests, the actual initial condition for the Aneth Unit is unknown.The initial condition can be estimated by creating travel times through an assumed velocity model, but a one-dimensional velocity model may oversimplify the geology, creating many anomalous regions of low and high velocity that are not due to CO 2 saturation, but instead due to geologic conditions that were unknown or simplified.
There are two distinct zones of velocity change in the tomograms in Figure 12.The change to a lower velocity in the lower region may be caused by an increased CO 2 concentration in the Desert Creek reservoir.A decrease in velocity of that layer would have a proportionally larger effect on the average velocity of ray paths traveling to the lower receivers.The high velocity region may be caused by the velocity of that region being higher in reality than in the synthetic model.Neither conclusion can be made with confidence given the synthetic tests.
The synthetic test in Figures 7 and 9 shows the monitoring potential of tomography, while the synthetic test in Figures 8,10 and 11 shows the reality of tomography applied to a poor event and receiver arrangement.It may be possible to draw conclusions about the CO 2 injected into the Desert Creek reservoir with this event and receiver arrangement, but it would require a more detailed velocity model, and would likely be limited to quantifying the amount of CO 2 injected rather than determining the extent of its migration.The position of a plume in the Desert Creek reservoir, or elsewhere, cannot be determined with reliability as the low velocity region will be stretched between the events and receivers, given their existing arrangement.
Both the Desert Creek reservoir and the southern cluster of events are positioned near a velocity layer change, which makes differentiation between CO 2 induced velocity changes and stretching of existing velocity reductions difficult.Four velocity layer boundaries currently exist in the model, at depths of 1050, 1400, 1790 and 2050 m.There are no clear velocity anomalies localized to these boundaries, suggesting that significant change in velocity between layers is not adversely affecting the ability to locate a plume.

Conclusions
The purpose of this project is to contribute to determining the feasibility of using double difference seismic tomography for monitoring, verifying, and accounting for geologically sequestered CO 2 .The three components of MVA to be addressed are: • Monitoring both the location and impact of sequestered CO 2 ; • Verifying movement of CO 2 and ensure that sequestered CO 2 is not permeating the sealing mechanism, or migrating to an unsealed area; • Accounting for the amount of injected CO 2 , by comparing it to the amount of CO 2 estimated to be in place.
The potential for geologic carbon sequestration monitoring applications through tomography may exist.Overall velocity changes before and after injection can be used to make estimates of injection volume in a specific region.This would require a control model before injection or be entirely dependent on an accurate initial velocity model.
Verification of injected CO 2 into the Aneth Unit by means of tomography may not be possible with the existing event and receiver arrangement.Substantial stretching of the low velocity zone outside of the reservoir extents would make detection of any leakage or migration nearly impossible.
The amount of CO 2 injected into the Aneth Unit cannot be accounted for with the current event and receiver arrangement.One reason for this inability to account for the CO 2 is the large size of the plume.Determining the amount of CO 2 underground requires an accurate measure of the influenced region.In addition to this exaggerated size, the amount of CO 2 injected can also not be determined because the lateral extents of the observable region do not include the entire affected region.The horizontal extents of the low velocity zones would have to be known to make any determination of the amount of CO 2 that can be accounted for.The best way of achieving this plume boundary characterization is likely to implement a surface array of receivers that extend beyond the injection region or using more than one borehole for monitoring.
Perhaps the greatest limitation of seismic tomography is the difficulty in obtaining large numbers of omnidirectional ray paths that traverse the full area of interest.The further removed tomographic monitoring is from that idealized scenario, the more accuracy and precision of results will suffer.The feasibility of using passive seismic tomography will vary to great effect depending on the site geometry and adaptive ability to design a monitoring array that is able to provide comprehensive reservoir coverage.The results presented suggest that a single vertical array of boreholes is not ideal for monitoring considering the passive seismic event locations.With the existing event and receiver arrangement, double difference seismic tomography is a poor choice as a method of MVA for the Aneth Unit.Monitoring is possible to a limited extent, but precise and accurate locations for a CO 2 plume cannot be determined.A more comprehensive event-receiver arrangement is required for accurate tomographic imaging of the Desert Creek reservoir.

Figure 2 .
Figure 2. Depth of the top of the Desert Creek reservoir [12].

Figure 3 .
Figure 3. Frequency of events per month.

Figure 5 .
Figure 5. Aneth P wave sonic log from the salt water disposal well (converted from [15]) (left) and derived velocity model (right).

Figure 6 .
Figure 6.(a) Cross section of the initial velocity model for the first synthetic test using a random distribution of 1000 events.Eight receivers are placed on the vertexes of the cube and 42 are distributed randomly throughout the rest of the cubic velocity model.And (b) Cross section of the initial velocity model for the second synthetic test using the Aneth event and receiver locations.

Figure 7 .
Figure 7. Tomogram of the change in velocity for the first synthetic test, using the velocity model in Figure 6a.

Figure 8 .
Figure 8. Tomogram of the change in velocity for the second synthetic test, using the velocity model in Figure 6a at a northing slice of −290 m.

Figure 9 .
Figure 9. Tomogram showing the checkerboard velocity model with an ideal event and receiver arrangement.

Figure 10 .
Figure 10.Tomogram showing the checkerboard velocity model using the events and receivers from the first time period at a northing slice of −290 m.

Figure 11 .
Figure 11.Tomogram showing the checkerboard velocity model using all the events and receivers at a northing slice of −290 m.

Figure 12 .
Figure 12.Tomograms showing velocity change between Time 1 and Synthetic Time 1 (top left), Time 2 and Synthetic Time 2 (top right), Time 3 and Synthetic Time 3 (bottom left), and Time 4 and Synthetic Time 4 (bottom right).Gray coloring represents no ray path coverage, and all tomograms are at the midpoint between the receivers and the centroid of the event cluster, a northing slice of −290 m.

Table 1 .
Events contained in each time period.