**Abstract:** GRAVSOFT is a package of FORTRAN programs for gravity field modelling. The
specific use of the programs for the determination of the empirical covariance function and its
analytic modelling and geoid determination using least-squares collocation is described. An
example using data from New Mexico, USA, illustrates the use of the method.

**1.Introduction.**

The purpose of these notes are to provide a guide to gravity field modelling, and especially to
geoid determination, using least-squares collocation (LSC). The use of the GRAVSOFT
package of FORTRAN programs will be explained. The readers are supposed to be familiar
with LSC as described in Moritz (1980). However the theory will be reviewed briefly in order
to fix the terminology.

So, for example when using the term **geoid**, we will mean the quasi-geoid, i.e. the so-called
height anomaly evaluated at the surface of the Earth.

A general description of the GRAVSOFT programs are given in Tscherning et al (1994),
which is updated regularly when changes to the programs have been made. It is available on
computer-readable form as a part of the GRAVSOFT-files.

The general methodology for (regional or local) gravity field modelling is as follows:

A: Transform all data to a global geodetic datum (GRS80/WGS84)

B: Use the remove-restore method. Here the effect of a spherical harmonic

expansion and of the topography is removed from the data and subsequently

added to the result. This is used for all gravity modelling methods

including LSC. This will produce what we will call **residual** data.

C: Determine at statistical model (a covariance function) for the residual data in

the region in question.

D: Make a homogeneous selection of the data to be used for geoid determination

using rule-of-thumbs for the required data density, check for gross-errors

(make contour map of data), verify error estimates of data,

E: Determine using LSC a (regional) gravity field approximation. Compute estimates

of the geoid heights and their errors.

F: If the error is too large, and more data is available add new data and

repeat E.

G: Check model, by comparison with data not used to obtain the model.

The whole process can be carried through using the GRAVSOFT programs GEOCOL, EMPCOV, TC (TC1), TCGRID, COVFIT, SELECT and GEOIP. GRAVSOFT includes data from New Mexico, USA., which can be used to test the programs and procedures.

They have here been used to illustrate the use of the programs. The programs have been used
for geoid determination by means of least-squares collocation in several countries and regions
(Andreu & Simo, 1992, Arabelos, 1989, Ayhan, 1993, Benciolini et al. 1984, Catalao &
Sevilla, 1994, Gerrard, 1990, Gil et al. 1993, Sevilla & Rodriques, 1993, Tscherning, 1985,
Tscherning & Forsberg, 1986).

**2. Theory.**

The anomalous gravity potential, T, is equal to the difference between the gravity potential W and the so-called normal potential U, T = W-U. T is a harmonic function, and may as such be expanded in solid spherical harmonics

and GM is the product of the gravitational constant and the mass of the Earth, and

The functions Y_{nm}(P) are orthogonal base-functions in a Hilbert space with an isotropic inner-product, harmonic down to a so-called Bjerhammar-sphere totally enclosed in the Earth. T
will not necessarily be an element of such a space, but may be approximated arbitrarily well
with such functions. This justifies the equality sign in equation (1).

The space will have a so-called reproducing kernel, which is a function of two points in space P, Q. (The coordinates of Q will be distinguished from these of P with an '). The kernel will be

where is the spherical distance between P and Q, P_{n} the Legendre polynomials and _{n} are
positive constants, the so-called (potential) degree-variances.

It is well -known, that if the degree-variances are selected equal to simple polynomial
functions in the degree n multiplied by exponential expressions like q^{n}, where q < 1 then
K(P,Q) may be represented by a closed expression. The most simple example are where _{n}=q^{n},
where we get well-known expressions related to the reciprocal distance.

In a reproducing kernel Hilbert space (RKHS) one may determine approximations to the
elements from data for which the associated linear functionals are bounded. The relationship
between the data and T are expressed through functionals L_{i}, so that

where y_{i} is the i'th data element, L_{i} the functional, e_{i} the error, A_{i} a vector of dimension k and

X a vector of parameters also of dimension k. A_{i}^{T}X may for example express the effect of a
datum-shift or of a bias and tilt in the data.

For the data we will meet (cf. Table 1) we have

or mean values of theses quantities. (Mean values will in practice be represented a means of point-values). The point Q is a point with latitude and longitude equal to P, but having ellipsoidal height equal to the orthometric height of P.

---------------------------------------------------------------------------

Table 1. GRAVSOFT data from the New Mexico Area to be used in LSC geoid test.

Type Format Error File name

Free-air gravity number, , , height, g (mgal) 0.2 nmfa

Defl. vertical number, , , height, , (arcsec) 0.5 nmdfv

DTM Grid label, grid data (m) 25 m nmdtm

---------------------------------------------------------------------------

The linear functionals are then

shown here in spherical approximation. (This approximation is only used on the residual
quantities, see section 3).

An optimal approximation to T using error-free data in a geocentric system may then be
obtained using that the observations are given by, L_{i}(T) = y_{i}. The "optimal" solution is the
projection on the n-dimensional sub-space spanned by the so-called representers of the linear
functionals, L_{i}(K(P,Q)) = K(L_{i},Q). The projection is the intersection between the subspace and
the subspace which consist of functions which agree exactly with the observations, L_{i}(g)=y_{i}.
(See Figure 1).

Figure 1.
Collocation estimates the
"point" with shortest distance from zero agreeing with the observations, and is simultaneously the orthogonal projection of T on the subspace spanned by the basefunctions K(L_{i},P),
if T is an element of the space.

Then

with

Based on the gravity field model given by equation (7), any quantity associated with a linear functional, L, may then be computed

If the data contain noise, then the elements _{ij} of the variance-covariance matrix of the noise-vector is added to K(L_{i},L_{j}). The solution will then both minimalize the square of the norm of T
and the noise variance. If the noise is zero, the solution will agree exactly with the observations. This is the reason for the name collocation. Upper limits for the approximation error
may be calculated if the norm of T is known,

This is not very useful, since the upper limit generally is very pessimistic (see Tscherning,
1984).

If we want to minimalize the mean-square error, the reproducing kernel must be selected equal to the so-called empirical covariance function, COV(P,Q). This function is equal to the reproducing kernel given in eq. (2), and having the degree-variances equal to

Here the normal equation matrix

result in predictions

and error estimates

In the diagonal of the normal equation matrix we find the sum of the data variance and the
noise variance. We can say that we here have obtained a balance between the signal and the
noise.

The covariances are computed using the "law" of covariance propagation, i.e. COV(L_{i},L_{j}) =
L_{i}(L_{j}(COV(P,Q))), where COV(P,Q) is the basic "potential" covariance function. COV(P,Q) is
an isotropic reproducing kernel with the degree-variances given by eq. (10).

**Example.** If we want to derive the gravity anomaly covariance function for gravity anomalies
in two points P and Q then we must appy the functionals given in eq. (6) on K(P,Q). This is

The quantities COV(L,L), COV(L,L_{i}) and COV(L_{i},L_{j}) may all be evaluated by the sequence of
subroutines COVAX, COVBX and COVCX which form a part of the programs GEOCOL and
COVFIT for the functionals given in eq. (6). For further details see Tscherning(1976, 1993).

**3. The remove-restore method.**

The least-squares collocation solution is giving the minimum mean square error in a very
specific sense, namely as the mean over all data-configurations which by a rotation of the
Earth's center may be mapped into each other. So if this should work locally, we must make
all areas of the Earth look alike, see from the gravity field standpoint.

This is done by removing as much as we know, and later adding it. We obtain a field which is
statistically more homogeneous than before.

First we may remove the contribution T_{s} from a known spherical harmonic expansion like the
OSU91A field, complete to degree N=360 (Rapp et al.(1991)). Secondly we may remove the
effect of the local topography, T_{M}, using Residual Terrain Modelling (RTM, see Forsberg &
Tscherning, 1981). We will then be left with a **residual** field, with a smoothness in terms of
standard deviation of gravity anomalies between 50 % and 25 % less than the original
standard deviation.

The residual quantities are then

**Exercise 1. **

Compute residual gravity anomalies and deflections of the vertical using the OSU91A spherical harmonic expansion and the New Mexico DTM, cf. Table 1. The free-air gravity anomalies are shown in Fig. 2.

The program GEOCOL may be used to subtract the contribution from OSU91A. A job-file with the name jobosu91.nmfa is shown in the Appendix 1. An example of a run is shown in Appendix 2 . The coefficients of the OSU91A model are found in the file osu91a1f. The difference files should be denoted nmfa.osu91 and nmdfv.osu91. The gravity anomaly differences are shown in Fig. 3.

The RTM contribution may be computed and subtracted using the program tc1. First a reference terrain model must be constructed using the program TCGRID with the file nmdtm as basis. A jobfile to run tc1 named jobtc.nmfa is found in Appendix 1. The result should be stored in files with names nmfa.rd and nmdfv.rd, respectively. The residual gravity anomalies are shown in Fig. 4.

The degree-variances will be changed up to the maximal degree, N, of the spherical harmonic series, contingently up to a smaller value, if the series is not agreeing well with the local data (i.e. if no data in the area were used when the series were determined). The first of N new degree-variances will depend on the error of the coefficients of the series. We will here suppose that the degree-variances at least are proportional to the so-called error-degree-variances,

The scaling factor must therefore be determined from the data (in the program COVFIT, see
later).

**4. Covariance function estimation and representation.**

The covariance function to be used in LSC is equal to

where is the azimuth between P and Q and , are the coordinates of P. Note that this is a
global expression, and that it will only dependent on the radial distances r, r' of P and Q and of
the spherical distance between the points.

In practice it must be evaluated in a local area by taking a sum of products of the data grouped according to an interval i of spherical distance,

where is the interval length (also denoted the sampling interval size). Note that two
intervals may be merged, so that the sampling interval becomes the double. In the program
EMPCOV this may be done a number of times, as an aid in selecting the right size of the
sampling interval.

Hence

where in actual calculations the covariance will refer to the mean height of the area.

In practice data directly related to T will not be available, but we may have geoid heights, gravity anomalies and deflections of the vertical. Covariance functions of geoid heights and gravity anomalies are estimated like for T, but the representation as a Legendre-series, cf. eq. (4) will be different. In spherical approximation we have already derived

where R is the mean radius of the Earth.

**Exercise 2.**

Compute the empirical gravity anomaly covariance function using the program EMPCOV for the New Mexico area both for the anomalies minus OSU91A and for the anomalies from which also RTM-effects have been subtracted (input files nmfa.osu91 and nmfa.rd). A sample input file to EMPCOV is called empcov.nmfa, see Appendix 1. A sample run is shown in Appendix 3 . The estimated covariances are shown in Fig. 5.

Note the correlation distance _{1}, i.e. the distance where the covariance becomes equal to 50%
of the variance C_{0}.

What is the mean height to which the calculated covariances refer ?

We see here, that if we can find the gravity anomaly degree-variances, we also can find the
potential degree variances. However, we also see that we need to determine infinitely many
quantities in order to find the covariance function.

The solution to this problem is to use a so-called degree-variance model, i.e. a functional
dependence between the degree and the degree-variances. In the program COVFIT, three
different models (1, 2 and 3) may be used. The main difference is related to whether the
(potential) degree-variances go to zero like n^{-2}, n^{-3} or n^{-4}. The best model (see Tscherning &
Rapp, 1974) is of the type 2,

where R_{B} is the radius of the Bjerhammar-sphere, A is a constant in units of (m/s)^{4}, B an
integer. If a spherical harmonic series expansion is used, B is typically put equal to a small
number like 4, while in the original work it was put equal to 24, so that the low-degree
degree-variances could be modelled appropriately.

The actual modelling of the empirically determined values is done using the program
COVFIT. The factors , A and R_{B} need to be determined (the first index N+1 must be fixed).

The program makes an iterative non-linear adjustment for the Bjerhammar-sphere radius, and linear for the two other quantities (see Knudsen, 1987). Unfortunately the iteration may diverge (e.g. result in a Bjerhammar-sphere radius larger than R). This will normally occur, if the data has a very inhomogeneous statistical character. Therefore simple histograms are always produced together with the covariances (in EMPCOV) in order to check that the data distribution is reasonably symmetric, if not normal.

**Exercise 3.**

Compute using COVFIT an analytic representation for the covariance function. An example
of an input file is found in covfit.nmfa, see Appendix 1. An example of a run is shown in
Appendix 3
. Gravity error-degree-variances for the OSU91A coefficients are found in the file
edgv.osu91. The estimated and the fitted covariance values are shown in Fig. 5.

Subsequently use COVFIT to tabulate the auto-and cross-covariance functions for residual
gravity, deflections and geoid heights. (Job file covfit.nmta, see Appendix 1). (Values at the
altitude of 1700 m are shown in Table 2.)

Table 2. Table of model-covariances of height anomalies, gravity anomalies and deflections of
the vertical. Height 1700 m. Model degree-variances equal to A/((i-1)(i-2)(i+4)). Error degree-variances used from degree 2 to 360 with scale factor 0.207. Error degree variances from
OSU91A (file osu91.edg). Depth to Bjerhammar-sphere = -6061.00 m, variance of point
gravity anomalies at 0 height = 196.20 mgal^{2}, the factor A, divided by R_{E}^{2} is = 452.52 mgal^{2}

Covariances between quantities KP and KQ, Evaluated in P,Q, having spherical distance PSI
and an azimuth equal to zero. Values of KP or KQ used signify 1=height anomaly, 3=gravity
anomaly, 6= component of deflection of the vertical in the direction between P and Q.
Units: m, mgal and arcsec.

KP= 1 3 3 6 6 6

KQ= 1 3 1 6 3 1

PSI(ddmm.m) m

^{2}mgal^{2}m*mgal arcsec^{2}arcsec*mgal arcsec*m

0 0.00 0.06524 142.22765 2.22634 3.17297 0.00000 0.00000

0 3.00 0.06411 131.37725 2.12867 2.81328 -7.00904 -0.08223

0 6.00 0.06096 104.86753 1.86979 1.95443 -11.71086 -0.14717

0 9.00 0.05640 73.50534 1.52289 0.97988 -13.58973 -0.18664

0 12.00 0.05112 44.55151 1.15654 0.13173 -13.29683 -0.20119

0 15.00 0.04573 21.11034 0.81785 -0.49983 -11.66534 -0.19571

0 18.00 0.04069 3.91214 0.53348 -0.90540 -9.36200 -0.17628

0 21.00 0.03630 -7.42000 0.31480 -1.11005 -6.85612 -0.14870

0 24.00 0.03270 -13.73315 0.16276 -1.15150 -4.45955 -0.11789

0 27.00 0.02994 -16.04179 0.07157 -1.07052 -2.36722 -0.08771

0 30.00 0.02794 -15.37807 0.03139 -0.90642 -0.68736 -0.06091

0 33.00 0.02661 -12.71152 0.03039 -0.69476 0.53671 -0.03925

0 36.00 0.02577 -8.90302 0.05627 -0.46611 1.31012 -0.02360

0 39.00 0.02528 -4.67778 0.09734 -0.24525 1.67322 -0.01405

0 42.00 0.02496 -0.61127 0.14338 -0.05094 1.69044 -0.01013

0 45.00 0.02469 2.87494 0.18609 0.10414 1.44048 -0.01095

0 48.00 0.02434 5.50874 0.21941 0.21302 1.00791 -0.01533

0 51.00 0.02384 7.15599 0.23955 0.27394 0.47589 -0.02200

0 54.00 0.02315 7.80293 0.24487 0.28952 -0.07969 -0.02969

0 57.00 0.02224 7.53425 0.23562 0.26588 -0.59476 -0.03726

1 0.00 0.02115 6.50831 0.21361 0.21158 -1.02003 -0.04376

1 3.00 0.01990 4.93120 0.18173 0.13657 -1.32262 -0.04848

1 6.00 0.01855 3.03134 0.14354 0.05124 -1.48616 -0.05103

1 9.00 0.01717 1.03626 0.10285 -0.03458 -1.50975 -0.05124

1 12.00 0.01581 -0.84740 0.06328 -0.11229 -1.40588 -0.04924

---------------------------------------------------------------------------

The numerical evaluation of the expression for the covariance function is rather time-consuming because it involves the summation of a Legendre-series up to degree N, here equal
to 360. However, the covariances of geoid heights, gravity anomalies and deflections of the
vertical may be tabulated. This option is available in COVFIT and in GEOCOL. However, the
selection of the optimal table entries is complicated. The solution recommended for the
moment is trial and error. The tool here is COVFIT, which may compute the differences
between tabulated and correctly evaluated quantities.

**Exercise 4.**

Use COVFIT to compute differences between tabulated and analytically calculated values. For
tabulation as a function of height use three intervals between 1500 m and 2600 m. For the
tabulation in spherical distance use three equidistant intervals bounded by 0" to 60", 60" to
360" and 360" to 3600" and with 20 sub-intervals in each main interval. An input example is
provided in the jobfile covfit.nmtt, see Appendix 1.

The covariances may also be tabulated in 1-dimension, i.e. only as a function of spherical
distance. This is very useful if we deal with ocean data (at height zero).

**5. LSC geoid determination from residual data.**

We now have all the tools available for using LSC: residual data and a covariance model. The
rest is to establish the normal equations eq. (11), solve the equations, eq. (12) and compute
predictions and error estimates, eq. (12) and (13). This may be made using GEOCOL.

However, as realized from eq. (8) we have to solve a system of equations as large as the
number of observations. This is one of the key problems with using the LSC method. The
problem may be reduced by using means values of data in the border area. Also, if the
observations are clustered, we may not need all observations. Rules for the necessary data
density (d) as a function of the correlation length _{1 }of the covariance function are given in
Tscherning(1984, p. 330). Suppose we want to determine geoid height differences with an
error of 10 cm over 100 km. This corresponds to an error in deflections of the vertical of 0.2".
This is equivalent to that we must be able to interpolate gravity anomalies with a mean error
of 1.2 mgal. The rule-of-thumb is

**Exercise 5.**

Use the residual gravity variance C_{0}, and the correlation distance determined in exercise 3 for
the determination of the needed data spacing.

Then use the program SELECT for the selection of points as close a possible to the nodes of a
grid having the required data spacing, and which covers the area of interest. The area covered
should be larger than the area in which the geoid is to be computed. Data in a distance at least
equal to the distance for which gravity and geoid becomes less than 10 % correlated, cf. the
result of exercise 3. Denote this file nmfa.rd1.

When data have been selected (as described in exercise 5) it is recommended to prepare a
contour plot of the data. This will show whether the data should contain any gross-errors. LSC
may also be used for the detection of gross-errors, see Tscherning (1991).

An input file for the program GEOCOL must then be prepared, or the program may be run
interactively. However, in order to compute height-anomalies at terrain altitude, a file with
points consisting of number, latitude, longitude and altitude must be prepared. This may be
prepared using the program GEOIP, and input from a digital terrain model.

**Exercise 6. **

Prepare a file named nm.h covering the area bounded by 33.0^{o} and 34.0^{o} in latitude and -107.0^{o} and -106.0^{o} in longitude consisting of sequence number, latitude, longitude and height
given in a grid with 0.1 degree spacing. Use the program GEOIP with input from nmdtm. This
will produce a grid-file. This must be converted to a standard point data file (named nmh2)
using the program GLIST.

When using GEOCOL there must then be specified

the coordinate system used (GRS80),

the spherical harmonic expansion subtracted (and later to be added),

the constants defining the covariance model and contingently its tabulation

the input data files (nmfa.rd or nmfa.rd1 if a selected subset is used)

the files containing the points in which the predictions should be made (nm.h2).

Several options must be selected such as wether error-estimates should be computed or
whether we want statistics to be output. The program may also produce a so-called restart file.
This file is an ASCII-file which contains input to GEOCOL which enamles the calculation of
predictions only. But it has the advantage that it may be used on different computers.

**Exercise 7.**

Run the program GEOCOL (geocol11) with the selected gravity data for the prediction of geoid heights and their errors in the points given by nm.h2. The result should be output to a file named nm.geoid. Predict also residual deflections of the vertical (nmdfv.rd) and compare with the observed quantities. Make sure a restart file is produced.

A model input file is found in jobnmlsc, see Appendix 1. An example of a run where all data
in a sub-area are used is found in Appendix 5 .

When the LSC-solution has been made, the RTM contribution to the geoid must be determined. Here the program tc1 may be used with the file nm.h defining the points of computation. The LSC determined residual geoid heights and the associated error-estimates are shown in Fig. 6 and 7.

Figure 6.

**Exercise 8.**

Compute the RTM contribution to the geoid using tc1 and add the contribution to the output
file from exercise 7, nm.geoid.

If mean gravity anomalies, deflections or GPS/levelling determined geoid-heights were to be
used, they could easily have been added to the data. It would not be necessary to recalculate
the full set of normal-equations. Only the columns related to the new data need to be added.
Likewise, an obtained solution may be used to calculate such quantities and their error-estimates.

**Exercise 9.**

Compute a new solution with the same observations as in exercise 7, but add as observation one of the predicted residual geoid heights. Define the error to be 0.01 m. Recalculate the geoid heights and the error-estimates. Use the possibility for re-using the Cholesky-reduced (see Tscherning, 1984) normal-equations generated in exercise 7. Verify that the error-estimates, which now are equivalent to error-estimates of geoid height differences, have a magnitude smaller than the one specified in exercise 5. (Error-estimates corresponding to one observed geoid height are shown in Fig. 8).

The use of deflections and geoid heights (e.g. from satellite altimetry) may require that
parameters such as datum shifts and bias/tilts are determined. These possibilities are also
included in GEOCOL (see Tscherning, 1984).

**6. Conclusion.**

We have now went through all the steps from data to predicted geoid heights. The exercises
describes the use of gravity data only, but observed mean gravity anomalies, GPS/levelling
derived height anomalies as well as deflections could have been used as well.

The difficult steps in the application of LSC is the estimation of the covariance function and
subsequent selection of an analytic representation.

The flexibility of the method is very useful in many circumstances, and is one of the reasons
why the method has been applied in many countries. If the reference spherical harmonic
expansion is of good quality, only a limited amount of data outside the area of interest are
needed in order to obtain a good solution. But if this is not the case, data from a large border-area must be used so that a vast computational effort is needed to obtain a solution. This may
make it impossible to apply the method.

A way out is then to use the method only for the determination of gridded values, which then
may be used with Fourier transform techniques (Schwarz et al., 1991) or Fast Collocation
(Bottoni & Barzaghi, 1993).

**Acknowledgement: **This paper was prepared for the International School for the determination and Use of the Geoid, Milano, Oct., 1994. Thanks to the US National geodetic Survey for
the permission to use the New Mexico data. Also many colleagues have contributed to or
inspired the delvelpment of the programs used: D.Arabelos (Thessaloniki), T.Krarup
(Horsens), K.Poder (Oelstykke), H.Suenkel (Graz), R.H.Rapp (Columbus, Ohio).

**References:**

Andreu, M.A. and C.Simo: Determination del geoide UB91 a Catalunya. Institut Cartografic de Catalunya, Monografies tecniques, Num. 1, 1992.

Arabelos, D.: Gravity field approximation in the area of Greece with emphasis on local characteristics. Bulletin Geodesique, Vol. 63, pp. 6984, 1989.

Ayhan, M.E.: Geoid determination in Turkey (TG-91). Bulletin Geodesique, Vol. 67, pp. 10-22, 1993.

Benciolini, B., L.Mussio, F.Sanso, P.Gasperini and S.Zerbini: Geoid Computation in the Italian Area. Boll. di Geodesia e Sc. Affini, Vol. XLIII, n. 3, pp. 213-244, 1984.

Bottoni, G.P. & R.Barzaghi: Fast Collocation. Bulletin Geodesique,Vol. 67, pp. 119-126, 1993.

Catalao, J.C. & M.J.Sevilla: Geoid comoutation in the North-East Atlantic (Acores-Portugal). Mare Nostrum - Geomed Report no. 4, pp. 77-90, Thessaloniki, Sept. 1994.

Forsberg, R. and C.C.Tscherning: The use of Height Data in Gravity Field Approximation by Collocation. J.Geophys.Res., Vol. 86, No. B9, pp. 7843-7854, 1981.

Gerrard,S.: The geoid, GPS and levelling. Thesis, University of Nottingham, Institute of Engineering Surveying and Space geodesy, 1990.

Gil, A.J., M.J.Sevilla, G.Rodriguez-Caderot: Geoid determination in Central Spain from gravity and height data. Bulletin Geodesique, Vol. 67, pp. 41-50, 1993.

Knudsen, P.: Estimation and Modelling of the Local Empirical Covariance Function using gravity and satellite altimeter data. Bulletin Geodesique, Vol. 61, pp. 145-160, 1987.

Moritz, H.: Advanced Physical Geodesy. H.Wichmann Verlag, Karlsruhe, 1980.

Rapp, R.H., Y.M.Wang and N.K.Pavlis: The Ohio State 1991 Geopotential and Sea Surface Topography Harmonic Coefficient Models. Rep. of the Dep. of Geodetic Science and Surveying n. 410, The Ohio State University, Columbus, 1991.

Schwarz,K.P., M.G.Sideris and R.Forsberg: The use of FFT techniques in physical geodesy. Geophys. J. Int., Vol. 100, pp. 485-514, 1990.

Sevilla, M.J. & G. Rodriguez Velasco: Preliminary determination of a gravimetric geoid in Portugal. Mare Nostrum - Geomed Report no. 3, pp. 139-150, Milano, 1993.

Tscherning, C.C.: A FORTRAN IV Program for the Determination of the Anomalous Potential Using Stepwise Least Squares Collocation. Reports of the Department of Geodetic Science No. 212, The Ohio State University, Columbus, Ohio, 1974.

Tscherning, C.C.: Covariance Expressions for Second and Lower Order Derivatives of the Anomalous Potential. Reports of the Department of Geodetic Science No. 225, The Ohio State University, Columbus, Ohio, 1976.

Tscherning, C.C.: Geoid Modelling using Collocation in Scandinavia and Greenland. Marine Geodesy, Vol. 9, no. 1, pp. 1-16, 1985.

Tscherning, C.C.: Local Approximation of the Gravity Potential by Least Squares Collocation. In: K.P.Schwarz (Ed.): Proceedings of the International Summer School on Local Gravity Field Approximation, Beijing, China, Aug. 21 - Sept. 4, 1984. Publ. 60003, Univ. of Calgary, Calgary, Canada, pp. 277-362, 1985.

Tscherning, C.C.: The use of optimal estimation for gross-error detection in databases of spatially correlated data. Bulletin d'Information, no. 68, pp.79-89, Bureau Gravimetrique International, 1991.

Tscherning, C.C.: Computation of covariances of derivatives of the anomalous gravity potential in a rotated reference frame. Manuscripta Geodaetica, Vol. 18, no. 3, pp. 115-123, 1993.

Tscherning, C.C., P.Knudsen and R.Forsberg: Description of the GRAVSOFT package. Geophysical Institute, University of Copenhagen, Technical Report, 1991, 2. Ed. 1992, 3. Ed. 1993, 4. ed, 1994.

Tscherning,C.C. and R.Forsberg: Geoid determination in the Nordic countries from gravity and height data. Boll. di Geodesia e Sc. Aff, Vol. XLVI, pp. 21-43, 1986.

Tscherning, C.C. and R.H.Rapp: Closed Covariance Expressions for Gravity Anomalies, Geoid Undulations, and Deflections of the Vertical Implied by Anomaly Degree-Variance Models. Reports of the Department of Geodetic Science No. 208, The Ohio State University, Columbus, Ohio, 1974.

**Appendix 1.**

Here are shown examples of job-files used when solving the excercises. They should be
regarded as a help when running the programs interactivally.

**Job to subtract OSU91A contribution from gravity**. Execution of the job is shown in
appendix 2 . jobosu91.nmfa.

geocol11<<!

f

f T f F f t f

5

OSU91A

3.98600500E+14 6378137.0 -484.1655 360 f f T F

(2i3,2d19.12)

osu91a1f

F F t

-1 2 3 3 4 5 0 13 -1 0.00 t F F F F f t f F t

nmfa

25

nmfa.osu91

5.0

T F

**Job to compute and subtract RTM effects from gravity**. jobtc.nmfa.

tc1<<!

nmfa.osu91

nmdtm

nmdtm5

nmdtm30

nmfa.rd

-1 4 1 1 2.67

31.52 34.98 -107.98 -105.02

18.0 999.0 2.67

2 2

-1

**Job to compute empirical covariance function of gravity**. empcov.nmfa. Execution of job
is shown in appendix 3 .

empcov<<!

MN GRAVITY DATA - OSU91A -TC .

5.0 96 3 F T F

nm91.covt

4800 9 T 3 2 0 5.0 f

3 3 0

nmfa.rd

T

**Job to determine parameters for covariance kernel**. covfit.nmfa. Execution of job is shown
in appendix 4 .

covfit<<!

f

4 f f t f

2

4

-6.150 200.0 360 f t t

0 2 0.27 F

osu91.edg

10 0.5 0.2 0.1

1 1 000 000 T

18 3 3 1700 1700 1 1.0 t

177.8 177.8 40.0 42.0 0.2 100.0 103.0 0.2

nm91.covt

**Job to tabulate covariance functions**. covfit.nmta. Result is found in Table 2.

covfit<<!

t

1 f f t f

2

4

-6.061 196.2 360 f t t

0 2 0.207 F

osu91.edg

25 3 0 0 0.0 t

1 1 1700 1700 F 3 3 1700 1700 f 3 1 1700 1700 f

6 6 1700 1700 f 6 3 1700 1700 f 6 1 1700 1700 t

**Job to control tabulation error**. covfit.nmtt.

covfit<<!

t

3 f f t f

2

4

-6.061 196.2 360 f t f

0 2 0.207 F

osu91.edg

3 1600.0 2600.0 4

60.0 300.0 1750.0 4000.0

10 20 30 40

25 3 0 0 0.0 t

1 1 1700 1700 F 3 3 1700 1700 f 3 1 1700 1700 f

6 6 1700 1700 f 6 3 1700 1700 f 6 1 1700 1700 t

**Job to predict geoid from residual gravity** and compare predicted and observed deflections
of the vertical. jobnmlsc. Execution of job is found in
appendix 5 .

geocol11<<!

t

F T T F f F t

t f f f f

nmrestart

nmneq1

F T f

t

5

OSU91A

3.98600500E+14 6378137.0 -484.1655 36 f f T F

(2i3,2d19.12)

osu91a1f **(use local file name !!) **

2

4

-6.061 196.20 360 F f T f

0 4 0.207227

osu91.edg

t

1 2 3 3 4 7 0 -13 -1 1700.00 f F F t F f t t F t

nmfa.rd

25

5.0

32.5 34.5 -107.5 -105.5

0.2

t

T f

F F F

F t f

1 2 3 3 5 0 0 11 -1 1700.0 t f F F F F f T F T

nm.h2

27

nm.geoid

33.0 34.0 -107.0 -106.0

t

f

F t t

1 2 3 3 4 5 6 -5 -1 1700.0 f F F F F F T T F T

nmdfv.rd

20

0.2

33.0 34.0 -107.0 -106.0

0.6

-1.0

t

t

Updated 2000-05-11 by cct