EVALUATION OF TERRAIN CORRECTIONS THROUGH FFT AND CLASSICAL INTEGRATION IN TWO SELECTED AREAS OF THE ANDES AND THEIR IMPACT ON GEOIDAL HEIGHTS

As part of the regional geoid modeling, terrain corrections were computed in Tierra del Fuego island and in the west side of the province of Mendoza. The first place is located in the southernmost region of Argentina and Mendoza is in the center-west of this country. Considering both study areas, elevations range from 0 m to 6500 m. The classical integration of prism contribution and the 2-D FFT technique were used to estimate terrain corrections. This study aims at discussing the results obtained by both approaches and their applicability considering their advantages and disadvantages according to the regions under investigation. The analysis allowed us to conclude that classical integration has a better performance than FFT methods, especially in the highest regions where terrain corrections can be overestimated in more than 20 mGals by FFT. Both techniques described show similar results in flat areas. Finally, the effect that the error of terrain corrections computation has on geoidal heights is also discussed and numerically tested. It is proved that an error in gravity anomalies of 20 mGals may cause up to 2 m geoid height error. Evaluation of terrain corrections through FFT and ... Bol. Ciênc. Geod., sec. Artigos, Curitiba, v. 19, n 3, p.407-419, jul-set, 2013. 4 0 8


INTRODUCTION
Terrain corrections (TCR) constitute an important part in the reduction scheme of gravimetric data, particularly in mountainous regions.They are a useful tool for smoothing gravity anomalies and their application allows the use of interpolation techniques in order to calculate a geoid model.
There is no single way to determine TCR.Broadly speaking, it is possible to distinguish two cases: methods which involve the use of prisms and estimate TCR through a classical integration and those employing FFT (Fast Fourier Transform) techniques.
Classical integration has the advantage of being straightforward to interpret while FFT algorithms are faster.The estimation of TCR can take hours or even days of processing in the case of the classical integration due to how the gravitational contribution of each prism is integrated.
The speed of computation depends mainly on the resolution of the Digital Elevation Model (DEM) employed, its extension and the capacity of the CPU processor.
In practice, the FFT is often given preference to the classical approach due to its higher integration speed.Regarding the computed values, the differences between both techniques are negligible in flat regions.However, in high mountainous areas they can exceed tenths of mGals.The reason for these discrepancies relies mainly on the steep topographic gradients which sometimes exceed 45° and cause the Fourier series to fail to converge (KIRBY and FEATHERSTONE, 2001).It must be clearly stated that this kind of topography can cause any method to give wrong results due to the rough terrain or the inadequate representation of the surface relief by the DEM.
In order to apply FFT despite the convergence problem, usually the resolution of the DEM is decreased leading to smaller topographic gradients.This smoothing can imply a resolution of 1' or less and can be considered a safe strategy in regions of moderate gradients.In mountainous areas, the denser the DEM, the better the representation of the highest frequencies.It is then necessary to reach a compromise between the best representation of the topography and the applicability of the TCR computation methods.
Two GRAVSOFT routines, namely TC and TCFOUR (TSCHERNING et al., 1992) were applied to quantify the effect of using classical integration or FFT for terrain correction estimations.
One focus lies on Tierra del Fuego island.In this region the elevation ranges from 0 to 2600 m.The second test area is situated in the province of Mendoza.
Here, the elevations range between 200 and 6500 m.Both regions represent parts of the Andean Cordilleras, one of the highest mountain ranges on Earth (Figure 1).
The following sections are devoted to the explanation of each method and their results.The practical relevance of the results obtained is demonstrated by a numerical estimation of the effect of terrain correction errors on geoid determinations for the area of Tierra del Fuego.TC estimates the effect of topography with respect to the Bouguer plate in a conventional manner.It integrates the gravitational effect produced by each prism constructed from the DEM (FORSBERG, 1984).
For a prism bounded by Δx= xp-x, Δy=yp-y, Δz=hp-z, the gravitational effect of terrain at a point P is given by equation ( 1 where ρ is the density, G is the gravitational constant and r is the spherical geocentric radius. (xp,yp,hp) and (x,y,z) represent the Cartesian coordinates of the computation point of height hp and the corresponding coordinates of mass elements respectively.
The total TCR is given by the sum of the terrain effect from all prisms.Due to the instabilities of these formulas over large distances, McMillan (1958) formulas are used in such case (equation 2).
(2) α, β and γ depend on the distance between the mass element and its distance to point P.For more details, see McMillan (1958).
One way to speed up the calculation is to use two DEM grids of different resolution: one high-resolution DEM covering the region close to the point of calculation and a less dense grid to take into account the gravitational effect caused by distant structures.For this purpose, two different radii must be defined.
When the proximity to the point of computation, P, is small enough (3x3 grid points), the calculation is performed using bilinear interpolation (TSCHERNING et al., 1992).
The described expressions are the formalism applied by the TC program to compute TCR.

Terrain Corrections Determined by FFT (TCFOUR)
Under a planar approximation, TCR at point P is written as equation (3).

∫ ∫∫
In the FFT approach, the solution is based on splitting equation (3) into three convolutions.
Distance l can be written as equation ( 5): If l is written as (5), integration of equation ( 1) along the z component, leads to: A series development of the term in square brackets shows that it converges only if Δh/s <1.This implies that gradients must not exceed 45°.
This development can be transformed into three convolutions in terms of h and h 2 , such as TCR is calculated by equation ( 8 with f and f 0 , given by equations ( 9) and (10).
f 0 is a singular integral which presents no problem at zero distance when discrete data is used (FORSBERG, 1985).8) can be solved applying a two-dimensional FFT, as is the case of the TCFOUR routine used in this work.
As in the calculation with prisms, two different DEM grids may be introduced here for the same reasons stated above.
TCFOUR is a suitable routine and FFT methods are widely in use due to the speed of FFT computation (FORSBERG, 1985;SIDERIS, 1985;SCHWARZ et al., 1990).

Topographic Analysis
Two topographic maps of the corresponding regions are shown in Figures 2a  and 2b.Both were derived from SRTM3 data (Shuttle Radar Topography Mission, FARR et al., 2007).
When computing TCR in Tierra del Fuego, a 20 km radius was used for the interior zone to take into consideration the gravitational effect of the closest topography.A 160 km radius according to the Hayford 167 km zone was applied for the external part (FORSBERG, 1984).Both radii were obtained after experimentations regarding the relationship between the minimum and maximum radius to achieve a required TCR accuracy of 0.1 mGal.For the inner grid a 3" resolution was chosen while for the external grid a 30" resolution was applied.
The radii for the inner and outer zones of Mendoza were the same as those for Tierra del Fuego.However, the resolution was 6" for the inner part and 12" for the external.As can be noticed, this resolution was less than that of the SRTM3.Although there was no significant difference when using 3" resolution for the inner grid, a 6" grid provided the best result.These values were adopted while applying the techniques described above.
The adopted resolutions for the inner and outer grids differ in both regions because it is a compromise between the most detailed representation of the topography and the minimization of the effect of 45° gradients.These gradients are larger in the area of Mendoza because of the steep topography in the whole region under investigation.Under this situation, a lower resolution for the inner grid is mandatory.On the contrary, the outer grid requires higher resolution than that applied to the corresponding grid of Tierra del Fuego.
The heights do not exceed 2600 m in Tierra del Fuego, which makes possible to use a resolution of 3" for the inner grid, even when this region can be elevated.Due to the fact that the highest mountains are located in the south-west part and the remaining heights are of 500 m or less, it is reasonable to use a 30" resolution for the outer grid.This number was determined through a few experiments concerning the variation of both radii.
The density value used in both regions for the TCR determination was 2.67 g/cm 3 according to Lodolo et al. (2007) and Blitzcow et al. (2009).
In Figures 3a and 3b terrain correction maps for both study regions are shown.In both cases, TCR were computed using the TC program; the classical approach.
Terrain corrections in Tierra del Fuego reach 60 mGals and 85 % are under 35 mGals.
Figure 2 -Topographic maps of both regions under investigation.The topography includes heights ranging from 0 to 6500 m.
In the case of Mendoza's province it is observed that terrain corrections exceed 120 mGals and 85 % are under 28 mGals for the selected region.Kirby and Featherstone (2001) concluded that TCR should not exceed 75 mGals.However, their results were obtained in regions like Hawaii, but not in the Andes.As it is seen here and in Blitzcow et al. (2009), TCR can exceed 100 mGals in the highest parts of the Andes.Figures 4a and 4b show the results of TCR determination for both regions in analogy to Figures 3a and 3b, but based on the application of the 2D-FFT.A remarkable increase in the maximum values is observed in both areas: the maximums values of TCR in Tierra del Fuego and Mendoza are 80 mGals and 175 mGals, respectively.
Figure 4 -Map of Terrain Corrections using TCFOUR.
In the maps of topographic gradients (Figures 5a and 5b), it is possible to appreciate that the steeper slopes, often exceeding 45°, coincide with the highest mountains.The correlation coefficients between TCR and the topographic gradients are of 60 % in both cases.Figure 6 confirms that the greatest differences (more than 50 mGals) between both methods correspond to the highest zones.These differences also show a correlation with the largest gradients.This condition may imply that FFT solution is not the most satisfactory one in high mountainous areas even if the resolution of the DEM is lowered.FFT overestimates TCR values in the steepest areas due to the convergence problem.
In flat regions the difference between both methods are negligible like in the center-east part of Tierra del Fuego.This means that if the computation area is extended and flat, FFT can be advantageous because of its speed.

THE EFFECT ON GEOID DETERMINATION
Terrain corrections, as stated before, are part of the procedure of geoid determination when gravity data is used.Especially in mountainous regions, they may have a significant impact on the gravity anomalies.
Assuming that the observed gravity anomalies have been accurately estimated (better than 0.1 mGals) there are still two error sources related to the DEM accuracy and the error in the estimation of topographic effects.Then, the focus should be the effect that terrain correction errors have on geoidal height determinations.
It was demonstrated that the differences between both methods for terrain correction determinations exceed 15 mGals in the southwest part of Tierra del Fuego island (Figure 6a).In the case of Mendoza, they are also larger than 20 mGals in the highest region (Figure 6b).In both places the largest differences between both techniques cover extended high mountain regions.Those discrepancies can be regarded as a systematic error of nearly 20 mGal in those areas.
In order to estimate the numerical effect that a 20 mGals error in mountainous zones has on geoid estimation, a test employing the Equivalent Source Technique (EST) was done.The EST method attempts to solve the inverse gravimetric problem: the distribution of masses below a given surface (geoid or topography) which are responsible for the observed external gravitational potential is nonunique.The EST technique consists in the determination of the location and mass of point masses, the configuration of which is consistent to the observed gravity anomalies, geoidal undulations and other observables of the Earth gravity field.For further details, see Dampney (1969), Del Cogliano (2006) and Guspí et al. (2004).
In this case, 654 Helmert (Faye) anomalies distributed in the province of Tierra del Fuego together with the remove-compute-restore method were used.These anomalies resulted from gravimetric reduction and interpolation of 512 observed gravimetric anomalies distributed in the province (Figure 7).20 mGals were added to produce a bias of that magnitude in part of these 654 gravity anomalies.As it was shown in the previous section, the differences between methods were around 20 mGals in the hightest places of Tierra del Fuego.In this way, the impact of erroneous TCR derived by the FFT approach is simulated.The affected gravity anomalies were located on the south-west side of the province, below 54°S and to the west of 67.6°S.This procedure provided a solution which was compared to an undisturbed gravimetric geoid, without applying the 20 mGals bias.The differences in geoid undulation are shown in Figure 8.The maximum error reaches 1 m in the highest places.Since the difference between both methods is about 0.5 mGals in flat areas (Figure 6a), their contribution would be up to 0.05 m in geoidal heights.A different result was obtained when 20 mGals were systematically added to all the 654 gravity anomaly values.In this case, the maximum effect on geoidal heights was 2 m (Figure 9).The results shown in Figure 9 confirm the results published by Vaniček and Martinec (1994).They indicated that a 0.01 mGals systematic error in gravity anomalies produces 1 mm effect in geoidal height.
It would be of great interest to carry out in Mendoza the same experiences as for Tierra del Fuego but the coverage of gravimetric observed data is not sufficient to such purpose.
In Mendoza, where large gradients are distributed all over the region, the error introduced by FFT based TCR can be regarded as systematic.This may imply an error of 2 m or more in geoid determination, depending on the error of the resulting gravity anomalies

DISCUSSION AND CONCLUSION
TCR were computed in two different regions of Argentina with two different techniques.
In both regions the differences between both methods correlate with the topography.
The fact that FFT overestimates TCR by tenths of mGals in regions of rough topography, makes classical integration the method which provides the most reliable results in those places.It should be mentioned that the difference between both methods can be neglected in flat areas.
In the case of Mendoza's province, TCR computed by FFT are in excess by more than 50 mGals, which can produce up to 5 m error in geoid determination and may affect geologic interpretations.
A numerical test done in Tierra del Fuego determined up to 1 m geoidal undulation effect caused by a 20 mGals error in gravity anomalies located in the highest places.This error also affected geoidal height determinations in lower areas, in a sensible way.
When the 20 mGals constitute a systematic error in all gravity anomalies, there is a systematic effect of 1m to 2m in geoidal heights too.
Considering that TCR participate in the smoothing procedure of gravity anomalies, when TCR are needed for high-accuracy geoid modelling in high mountains, the use of classical integration instead of FFT is recommended.

Figure 1 -
Figure 1 -Map of the geographic location of the study regions which are marked within rectangles.

Figure 3 -
Figure 3 -Terrain corrections computed with TC program.

Figure 5 -
Figure 5 -Map of topographic gradients in both regions.

Figure 7 -
Figure 7 -Distribution map of observed gravity data in Tierra del Fuego.

Figure 8 -
Figure 8 -Distribution of differences after comparing the geoid solutions caused by a disturbed and an undisturbed set of gravity anomalies.The gravity anomalies located in the south-west region were biased by 20 mGals in one of the cases.The maximum difference in geoidal heights reached 1m.

Figure 9 -
Figure 9 -Distribution of differences after comparing the geoid solutions due to a disturbed and an undisturbed set of gravity anomalies.The gravity anomalies in the whole region were perturbed by 20 mGals in one of the cases.The error in geoidal heights was near 2 m.