The current terrain correction procedure corrects gravity measurements for the effect of terrain from a distance Rmin to a distance Rmax. Each gravity station is processed independently, and therefore corrections calculated for a particular station do not depend upon possible location errors of other stations. The correction procedure is designed to provide corrections which can be added to the "outer-zone" corrections calculated by existing terrain correction programs which are in common use in the United States, such as that of Plouff (1966). The distance Rmax is specified by the user during execution, and should be chosen to be large enough that the corrections calculated by the outer-zone correction program are accurate. The value selected for Rmax. therefore depends upon the smallest compartment size available to the outer-zone terrain correction program. Typically, that size is 15", 30", or 1' of arc, or roughly 450 to 1800 m. The distance Rmin is also supplied by the user; in areas of relatively low relief, values for Rmin of 5 m or so seem appropriate. In areas of very high local relief, the 30 m digitization cannot accurately represent the local terrain, and thus Rmin should be selected accordingly. A value for Rmin of 45 m [about the radius of Hammer's (1939) Zone C] seems to provide good results where local relief is large.
The correction procedure divides the circular area enclosed by Rmax into an inner zone and an outer zone; the radius separating the inner and outer zones is denoted Rmed. A surface is fit to the elevations between Rmin and Rmed; this surface is numerically integrated to calculate that portion of the terrain correction that is due to the terrain located in the interval Rmin <= R <= Rmed. The figure below illustrates the subdivision of areas used in the calculation.
The current version of the program uses the triangulation and interpolation procedures of Renka (1984) to interpolate elevations in this interval; earlier versions of the programs fit a multiquadric surface to the elevation data within the inner zone. The multiquadric surface approach, while quite accurate, was exceedingly slow and required substantial memory. Between Rmed and Rmax, the terrain effect is calculated by assuming that each elevation sample represents the elevation of a rectangular compartment 30 m on a side; a line element formula is used to calculate the effect of each such compartment. Terrain compartments lying between Rmed and Rmax that partially intersect the circular radius Rmed are treated in such a way that only the effect of that portion lying outside Rmed is added to the overall terrain effect. A similar procedure is applied to those compartments that partially intersect the outer radius Rmax
The use of multiquadric surfaces to fit the terrain data, while not the most efficient method from a computational standpoint, is a rather accurate method of representing topography (Hardy, 1971), and has been used for calculating terrain corrections for gravity stations by others (e.g, Krohn, 1976; Gettings, 1981). A comparison of several scattered data interpolation methods by Franke (1982) concludes that the multiquadric quation method, from the viewpoint of smoothness and accuracy, was one of the best scattered data interpolation methods of those that were tested. A disadvantage of using the multiquadric equation surface technique is that both time and memory requirements for fitting the multiquadric surface are proportional to N², where N represents the number of elevation samples used to construct the multiquadric surface; for the case under consideration, N is the number of DEM samples falling within or close to the innermost radius Rmin. Moreover, the evaluation of the surface at a single location requires summing the N terms defining the surface; as such a sum must be constructed many times for each innermost terrain correction, the evaluation phase of the calculation can take as long as the surface fitting phase. Finally, all the computations related to the construction and evaluation of the multiquadric equations should be performed in double precision on 32-bit machines.
Note that the numerical integration procedure used to calculate the effect of terrain lying between Rmin and Rmed essentially integrates the effect of a line element between the two radii. This integration is repeated every 5° to obtain the effect of all the terrain lying within the circular region. The radial portion of the integration is performed using an adaptive technique called QUAN8; this routine is described by Forsythe and others (1977, page 97). The technique is really the same as that of Gettings (1981); however, Gettings derived his method in an extremely complicated manner. An very simple derivation can show that the method used here is exactly equivalent.
A terrain surface is used close to the station location because the elevations provided by the DEM are samples and do not actually represent mean elevations of 30 m compartments. The use of a surface provides a terrain representation that should be much closer to reality than the use of compartments of constant mean elevation. At a certain distance from the station, the procedure of considering the elevation samples as representing the mean elevation of a compartment 30 m on a side should yield numerical results than are not distinguishable from the results that would be obtained by actually using mean elevations of compartments whose size would necessarily be larger than 30 m on a side. In fact, the method used is equivalent to numerically integrating the terrain effect (at distances greater than Rmed) using a rectangular rule. Because the compartment size is relatively small, use of the rectangular rule should be rather accurate. Typical values for Rmed and Rmax are 200 m and 2000 m. From a practical standpoint, the value of Rmed is restricted to relatively small values (usually less than 400 m), due to the increased memory and time required to fit and integrate the multiquadric surface to the elevation data. The value of Rmax can rather easily be increased to greater values, though, as no surface fitting is involved. However, little increase in accuracy of the calculation of the overall terrain effect is expected by increasing Rmax to values much greater than 2500 m.
The elevation of the gravity station is not directly used during the computation of the effect of terrain effects. Instead, the elevation at the horizontal location of the gravity station is calculated from the multiquadric representation of the terrain surface, and this calculated elevation is used in the computation of the effect of terrain. The actual elevation of the gravity station is not used at all.
The reason for not using the actual elevation of the gravity station for the computation of the terrain effect is simply that the DEM elevations often differ from the actual elevations by amounts of up to 7 m. The use of the actual elevation, rather than one consistent with the DEM terrain, effectively leads to gravity stations being located in deep holes or on very steep hills, whereas in fact such terrain features probably do not exist in the immediate vicinity of the gravity station. The departure of the DEM elevations from the actual elevations is a relatively local phenomenon, as the mean difference between the DEM and the actual terrain surface should be much closer to 0 than to 7 m.
Numerically, the procedures used for the calculation of the terrain effect are extremely accurate, especially when compared to terrain corrections calculated using template methods. However, the corrections calculated are no better than the terrain data that are used to represent the terrain about each gravity station. These data are typically at least as accurate as the elevation data present on standard USGS topographic maps, though, and thus it seems that the resulting terrain corrections will always be more accurate than those calculated using template methods, when the basic source of terrain data consists of USGS topographic maps.