Gravity Calculations from 3D
Geologic Models: A New Approach

Graham Brew, Dynamic Graphics Inc.

Summary

Described here is a new approach for calculating the gravity response of geologic models. This scheme has two major advantages over many previous methods. Firstly, the models and calculations are all fully three-dimensional—hence very complex models and arbitrary density distributions are permissible. Secondly, the models are built as part of a complete workflow using geologically-intuitive steps, and industry standard input data. Building of the models and population with density data is, in this case, performed within the EarthVision software.

The modeled density distributions are decomposed into regular grids, and the gravity contribution of each grid cell is calculated and summed—using well-established formulae—to produce the final anomaly. This approach is very computer-intensive, and currently the method is best suited to integration within a broader earth modeling workflow wherein gravity models are but one of a suite of outputs.

1. Objective

The goal of this project is to offer new features beyond what are often used in gravity modeling of earth models. Hence, the new approach should satisfy the following criteria:

Fully 3D geology: Construct models using geologically-intuitive processes, while maintaining full 3D flexibility⁄no along-strike restrictions.
Fully 3D density distributions: Models with any density distribution should be permissible. No requirements to use constants, gradients etc.
Integrated with earth modeling: Models will be built within a broader earth modeling workflow, as is common in the hydrocarbon industry. Tools from this workflow will be used in the gravity anomaly investigation.
Leverages modern visualization tools: Built within 3D software, the gravity analysis should utilize available 3D visualization capabilities.

2. Introduction to EarthVision

The core input to the gravity calculation will be a 3D grid of density values that was constructed within a geological⁄structural earth model. The population of the density grid is discussed in paragraph 4, and the gravity calculation in paragraph 5. This paragraph describes the software used to construct the geologic earth model—EarthVision from Dynamic Graphics, Inc. The attractions of using EarthVision include fully 3D workflows, powerful visualization capabilities, and customization options. Plus, EarthVision models are widely used in industry for volume calculations, well plans etc.

Fully 3D Workflows

Data can be input to EarthVision in most industry standard formats—and custom inputs can be generated if a format is not supported. EarthVision uses these data to construct fully three-dimensional models. These models can consist of almost any conceivable geologic geometry, no matter how complex. For example, overturned surfaces, thrust faults, isolated salt bodies, and channels are routinely modeled. Fault and horizon surfaces can be located anywhere, and at any orientation, within these 3D models with no arbitrary constraints on fault framework construction.

Furthermore, modeling building within EarthVision is accomplished using geologically-intuitive operations (e.g., faulting, deposition, erosion) that are familiar to any earth scientist. Hence, gravity models are built based on true structural concepts, rather than via arbitrarily shaped polyhedra. Models are also easily and quickly changed through logical user-interfaces. Perhaps most importantly, these models are easy to build, even by relatively inexperienced users.

Likewise, property models can have totally arbitrary geometry. Hence density distributions within a model are not artificially constrained to be constants or gradients—literally any density distribution is possible. Moreover, the powerful property modeling tools in the software can be used in the construction of the property distribution (e.g., gridding property distributions in pre-faulted space). Finally faults, horizons, and properties all considered simultaneously for building self-consistent 3D earth models.

Visualization

EarthVision includes many powerful tools for the editing, visualization, and interrogation of data, models, and anomalies. The diagram below shows just some of the potential for visualizing multiple datasets. Clearly, this can be very useful when comparing observed⁄calculated gravity anomalies, comparing models with anomalies in 3D, and when iterating the models to fit the observed gravity response. These ideas are further illustrated in later paragraphs.

Customization⁄Scripting

EarthVision software is built on a modular basis. Inputs and outputs to the modules can be handled by a user-defined scripting language. Moreover, outputs from EarthVision can be taken and manipulated in almost any conceivable way by these scripts. This puts tremendous power and flexibility in the hands of the user. For the case of gravity modeling, this means taking the density grid computed in EarthVision, calculating the gravity response of the density model, and sending the output back to EarthVision for visualization.

3. Modeling Workflow

A salt dome example illustrates a typical modeling workflow. Portions of this sequence are often iterative as new data and interpretations are added.

Data input: Initial input data consisted of well data and depth-migrated seismic. Time-migrated seismic data, together with a velocity model, could also be used and depth converted within EarthVision.
Data interpretation: Seismic and well data are interpreted either within EarthVision, or interpreted picks are imported from another package.
Fault surface modeling: Fault picks are gridded and surfaces created. For the case of dying faults, fault tip-lines can be limited by polygons, automatically or manually (see figure below).
Fault tree construction: Fault relationships and dependencies are built into a fault ″tree,″ either automatically or manually.
Horizon surface modeling: Horizon picks and well tops are gridded and surfaces created in the context of the fault framework. Hence, the resulting model is geometrically consistent and geologically sensible, even where input data are sparse.
Property modeling: Stratigraphic zones can be populated with property grids (e.g., density). This property information can easily be reconstructed into a pre-faulting configuration for more accurate modeling.
3D model construction: Faults, horizons, and properties are automatically built into a self-consistent 3D solid earth model.

4. Populating with Density Data

Density information is included in the 3D modeling workflow as a ″property″ of the geological layers. In the case of the salt dome example, these data came from well logs. Over twenty wells with density logs covering at least part of the hole were obtained, as shown below.

During property modeling, the input data are considered within the structural framework and hence can be reconstructed into unfaulted space before gridding—often yielding more desirable results. Also, several stratigraphic zones can be ″grouped″ during gridding, or constants specified for zones—see the salt dome example. Geostatisical tools for property prediction are also available.

Once the model is populated with density data, the property is decomposed into a regular 3D grid. The grid resolution is user-defined. This grid is the input to gravity calculation.

5. Calculating the Gravity Response

Once the density grid is constructed the gravity calculation can begin. The gravity anomaly at each point of observation is the sum of the individual anomaly contributions of all the grid cells (cuboids) in the model, as measured at the observation point. The contribution of a single cuboid, after Plouff (1976), is:

g = Gρ Σ_i=1,2 Σ_j=1,2 Σ_k=1,2 s [z_k tan^-1 (x_iy_j ⁄ z_kR_ijk) - x_i ln(R_ijk + y_j) - y_j ln(R_ijk + x_i)]

x₁ is the distance from the point of observation to the near edge of the cuboid along the x-axis,
x₂ is the distance to the far side of the cuboid;
y_1,2 and z_1,2 are similar measurements along the y and z axes.
ρ = density
G = universal gravitational constant
R_ijk = √ (x²_i + y²_j + z²_k)
s = s_is_js_k; s₁ = -1; s₂ = 1

An EarthVision script has been written that reads the density grid and various associated variables. The script then calculates the gravity contribution of each grid cell within a larger loop that considers each point of observation. The actual computation algorithms are a modification of those presented by Blakely (1995). A snapshot of the user interface to the gravity routine is shown here. Other inputs to the routine are the height at which to make the observations (currently a single value across the entire model), any bulk shift, and the ″background″ density—see paragraph 7. The algorithm only calculates relative gravity anomalies, not absolute values. Observed values (to compare against the calculated anomaly) will most likely be Bouguer anomalies or isostatic residuals.

Computation Time

Previously, elemental approaches such as this were rendered impractical due to prohibitively long compute times. Modern computing power, however, allows small models to calculate nearly instantly, and larger models usually complete in several minutes. Any delay is clearly unacceptable for ″interactive″ modeling (where changes in the model are simultaneously reflected in changes to the calculated anomaly). However, as stated above, the usage of this routine is currently intended as a parallel process to the model building, rather than as an interactive tool.

Testing

The gravity calculation routine was tested using a variety of simple polygonal cases, such as an infinite slab, small prisms, and other regular distributions. The anomalies were then compared to those computed directly from modification of the well-known Newtonian formulae (e.g., Telford et al., 1990). Furthermore, more complex shapes were tested at a variety of along strike locations and the results compared to publicly available implementations of 2.5D gravity computations (e.g., Roecker, 1997). The results were all found to agree with the other methods. Also see the later paragraph (number 8) for the Syria example of precision testing.

Notwithstanding the above, a key criteria in the approach discussed herein is using 3D density grid resolutions fine enough to capture the structural and property complexities within the geologic model.

6. Edge Effects

Any gravity modeling scheme must consider what ″exists″ at points beyond the model space. This is equally true in the current case where, in theory, the contributions of all points out to infinity in the x and y directions should be considered (since the anomalies are relative, consideration in the Z dimension can be limited to the thickness of the modeled volume). In practice, the consideration is limited to those volumes relatively close to the area of interest. Currently the routine adds an ″extended volume″—of background density—out to 10,000 times the model size in the x and y dimensions. In future versions of the routine, this added volume could have variable density, based on the chosen density sequence within the modeled space.

Testing of this approach (see paragraph 8) shows that there is no discontinuity between the gravity response over the model space and that in the extended volume.

7. Visualizing the Results

One of the great advantages to building and calculating the gravity model within a comprehensive 3D environment is the visualization power that can be used for analyzing and interrogating the result.

In this manner, input data, models, and output anomalies can be spatially compared and investigated.

A further advantage of using the broader modeling software is the availability of additional mapping tools. These include 2D map creation, isopach and cross-section generation, volumetics tools, and geostatistical processing power.

In the case of the salt dome example, gravity observations were not available. If they were, these too could be incorporated into the 3D visualization. Model adjustments could be made based on the comparison of the observed and calculated anomaly. Paragraph 8 shows an example where this was accomplished.

8. Regional Example from Syria

As a further example of the 3D gravity routine, an area from southern Syria has been selected for modeling. This location has a well-known gravity response that has previously been investigated using 2D profile methods.

The next two figures represent a first-order attempt to forward-model the observed gravity anomaly with a 3D geologic construction. Modeled area indicated with dashed box on map above.

The figure below shows a test of precision of the modeling process and gravity calculation. The densities of all the layers within the model were replaced by a single constant (2.75 Mg/m3). The ″background density″ also used this same number. Hence, the modeled geometry was effectively an infinite flat slab of constant density. The calculated anomaly shows departure from the expected constant of less than one part in 200 million.

9. Applications⁄Future Development

As currently implemented, this new approach is not a replacement for more traditional 2.5D profile-based gravity modeling. One limitation is that the 3D scheme is computationally too intensive for the interactive style of geometry iteration commonly used in forward modeling. Moreover, the distribution and quality of data needed to construct a 3D model are not usually available when 2.5D modeling is being considered (frequently the early stages of the exploration process).

However, the concept has been proved, and at least two areas of potential application are foreseen:

Gravity modeling to assist seismic analysis. Full 3D models are now a common product of the exploration and development workflow. However, despite huge advances in seismic technology, gravity data can still offer new insights in hard-to-image areas such are sub-salt and sub-basalt. Furthermore, potential fields can be used to better constrain the seismic velocity field in such areas.
Gravity calculation for regional models: Recent advances in software have lead to an increasing number of regional studies being accomplished in 3D. Often gravity data are available throughout the studied area, but are perhaps only utilized along profiles. The approach outlined here could easily be employed for validating such models in three-dimensions, and offering further geometric constraints where other data are lacking.

Future Developments

Speed: decrease compute time by automatically combining like-cells to make larger, faster-calculating prisms; parallel processing; code optimization.

Instantaneous Update: Code improvements to allow instantaneous, interactive manipulation of model to change anomaly. Other tools for more interactive workflow (dependent on speed improvements).

Spatially Variant Level of Observation: Include option to use a topographic (or other) surface, rather than a simple constant, as the observation height.

Inversion: Semi-automation to match calculated anomaly with observation.

Magnetics: add magnetic anomaly calculation to the routine.

References

Blakely, R.J., 1995. Potential theory in gravity and magnetic applications, Cambridge University Press, 441 pp.

Brew, G.E., Barazangi, M., Sawaf, T., & Al-Maleh, K., 2001. Tectonic and geologic evolution of Syria. Geoarabia, 6, 573-616.

Ponce, D.A., T.G. Hilterbrand, & R.C. Jachens, 2003. Gravity and magnetic expression of the San Leandro Gabbro with implications for the geometry and evolution of the Hayward Fault Zone, northern California. Bulletin Seismological Society of America, 93, 14-26. Also http://3d.wr.usgs.gov/docs/wgmt/3d/ressv2.html

Plouff, D. 1976. Gravity and magnetic field of polygonal prisms and application to magnetic terrain corrections. Geophysics, 41, 727-741.

Roecker, S, 1997. Xgrav: Interactive gravity modeling software, available via internet at: http://gretchen.geo.rpi.edu/software.html#sect2

Telford, W.M., L.P. Geldart, & R.E. Sheriff, 1990. Applied Geophysics, 2nd ed., Cambridge U. Press, 770 pp.

Acknowledgements

Thanks to all the staff at Dynamic Graphics, Inc. for their help and support during development. Thanks to Bob Jachens of the USGS for providing gravity data and models, and to Meghan Miller, formally of the USGS and DGI, for ideas and suggestions.

Gravity Calculations from 3DGeologic Models: A New Approach