Considerations on the further improvements of regional geoid modeling over the Baltic countries/Baltijos saliu regioninio geoido modeliavimo tobulinimas.

Ellmann, Artu

1. Introduction and motivation of the study

The geoid is an equipotential surface of the Earth's gravity field that coincides approximately with the mean sea level. Precise knowledge of the geoid contributes to geosciences and in solving of many engineering tasks. It supports many professional, economic and scientific activities and applications, such as navigation, mapping and surveying for the construction and maintenance of nationwide communications. The geoid is instrumental in geodetic infrastructure, as the topographic heights and the depths of the seas are reckoned from it. High-resolution geoid model, in particular, enables the user in many cases to replace the traditional height determination techniques, such as levelling, by faster and costeffective GPS measurements. Also the determination of the variations of the ocean currents and the interpretation of seismic disturbances benefit from the knowledge of this important reference surface.

Several recent space technologies have improved our knowledge of the global gravity field and Earth's topography. However, the space-borne gravity data is not only limited by its accuracy but also by the spatial resolution. For instance, the ongoing satellite gravimetric mission GRACE (Gravity Recovery and Climate Experiment) has resolved the long-wavelength component of the global geoid with an accuracy of a few cm, whilst the spatial resolution of such information is limited to about 200 km.

Even though the first satellite gradiometry mission GOCE (Gravity Field and Steady-State Ocean Circulation Explorer, launched by the European Space Agency in March 2009) will be capable to further enhance the intermediate wavelength information of the gravity field, but only up to the 65 km spatial resolution. Further improvements to the knowledge of the Earth's gravity field at shorter wavelengths should still come from the use of terrestrial surveys and satellite altimetry (over the oceans). These worldwide data-sets are being used by authorized research centers to develop models of the Earth's gravity field. Such Earth Geopotential Models (EGM) comprise a set of spherical-harmonic coefficients, which are obtained from the spectral analysis of the Earth's geopotential. Importantly, the spatial resolution of EGM-s is directly linked to the expansion degree of the Fourier series (the higher the degree the better the spatial resolution). Earlier EGM-s were developed only up to maximum degree as of 360. Note that a study by Ellmann, Jurgenson (2008) evaluated the quality of the EGM96 (Lemoine et al. 1998) and a GRACE-based EIGEN-GL04c (Forste et al. 2006) over the Baltic countries.

The resolution of a new combined EGM08 (Pavlis et al. 2008) is 5' arc-minutes (corresponding to 9 km, i.e. to the spectral degree of 2160). A study by Ellmann et al. (2009) evaluates the performance of the EGM08 model over the Baltic Sea region with emphasis to Estonia, Latvia and Lithuania. In particular, the EGM08-derived height anomalies were compared with an existing regional geoid model BALTgeoid-04 (Ellmann 2004; see also Section 4.2 of this paper). The detected discrepancies range within [+ or -] 0.3 m with a mean of -0.02 m, whereas the standard deviation (STD) of the discrepancies amounts to 0.08 m. The largest discrepancies occur in areas where only a few gravity data points were available either for the regional geoid modeling or at the EGM08 compilation, or both. The EGM08 model was also validated with respect to GPS-levelling data. After removal of the vertical offset (~ 0.5 m) the STD of detected discrepancies is 0.06 m. For more details see the original publication Ellmann et al. (2009).

Thus, for many applications the resolution and accuracy of the EGM08 may not be sufficient. For solving a large variety of engineering tasks a high-resolution (2-3 km) regional geoid model with an 1 cm accuracy is required. Obviously, due to tremendous computational burden and the voids of terrestrial gravity data it is unrealistic to develop such an ultra-high-degree spectral model of the global geoid. Therefore, the usage of the local terrestrial data is still requested for the high-resolution regional geoid modeling.

Alternatively, a geoid model can also be computed by using the Stokes integral formula from the global coverage of gravity anomalies (Stokes 1849). However, this method still remains impractical, due to the lack or limited access to the worldwide terrestrial gravity data. Therefore, regional improvements of the global geoid models can be obtained by modifying the original Stokes formula. This method was first proposed by Molodenskii et al. (1960) in the end of the 1950-ies. Their proposal coincided thus with the launch of the first artificial satellite--"Sputnik". A modified- Stokes formula combines local terrestrial gravity anomalies and the EGM-derived long-wavelength component (i.e. the Aglobal trend", the most reliable source of which is the satellite tracking data) of the geoid in a truncated Stokes's integral.

Since some recent studies (Ellmann, Jurgenson 2008; Ellmann et al. 2009) have already evaluated the performance of the global geopotential models in the Baltic Sea area, then it is appropriate to investigate other challenges in further improvements of the geoid modeling in the Baltic Sea region.

This study is described in six sections. The introduction is followed by a general review on modifications of Stokes's formula. Section 3 tackles earlier geoid modeling works over the Baltic countries. Section 4 describes the similarities and differences of the stochastic and deterministic modification methods. The emphasis is given to the selection of the most important geoid modeling parameters. The results of a numerical study are discussed as well. The terrestrial data evaluation results in Estonia are discussed in Section 5. A brief summary concludes the paper.

2. Modifications of Stokes's formula

2.1. General

Several different modification methods have been proposed in the geodetic literature over the past half-century. For computing a geoid estimator N a generalised Stokes modification scheme uses a modified Stokes's function and residual gravity anomaly in truncated integral (cf. Vanicek, Sjoberg 1991):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (1)

where R is the mean Earth radius, y is the normal gravity at the computation point, [psi] is the geocentric angle between the computation point and the integration element, d[sigma] is an infinitesimal surface element of the unit sphere [sigma], the integration area [[sigma].sub.0] is limited to some spatial domain (say, a spherical cap with radius [[psi].sub.0]) around the computation point. The modified Stokes function [S.sup.L]([psi]) is expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], (2)

where the modification parameters sk are selected by different criteria. The original Stokes function is denoted S([psi]), which can be represented via Legendre polynomials [P.sub.n](cos[psi]) as follows (cf. Heiskanen, Moritz 1967: 29)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. (3)

Thus Eq. (2) becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. (4)

Apparently, a truncation bias may occur due to neglecting the high-frequency (n > M) contribution of gravity anomalies located outside the integration domain (i.e. [[psi].sub.0] < [psi] [less than or equal to] n), see Eq. (1). The primary objective of the kernel SL([psi]) modification is to reduce the truncation bias to a level, which is acceptable for modern geodetic applications. For this the low degrees (2 [less than or equal to] n [less that or equal to] L) of the original Stokes function are modified (or simply removed), implying, in general, that across the integration domain [parallel][S.sup.L]([psi])[parallel] < [parallel]S([psi])[parallel], for an illustration see also Ellmann (2004: Fig. 3.3). Essentially, modification methods differ from each other by the selection of the modification parameters sk in Eq. (2). For instance, in the Wong, Gore (1969) modification approach the sk coefficients are a priori fixed to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], (5)

which is equivalent to the case when the summation in Eq. (3) starts from L + 1. In other words, the modified Stokes function [S.sup.L]([psi]) tapers off more rapidly than S([psi]), thus the contribution of distant gravity anomalies is expected to become manageably small.

The estimator Eq. (1) employs the high degree residual gravity anomalies, which are obtained by the subtraction the long-wavelength contribution of gravity from the complete anomaly [DELTA][uu]. It is understood thus, that the gravity anomaly can be expanded into a series of Laplace harmonics, i.e. [DELTA]g = [[infinity].summation over (n = 2)][DELTA][g.sub.n]. Apparently, due to the existence of various errors the terrestrial gravity anomalies [DELTA][uu] and the harmonics [DELTA][[uu].sub.n] are only estimates of their true values. The harmonics [DELTA][[uu].sub.n] can be calculated from an EGM by a standard formula (cf. Heiskanen, Mo ritz 1967: 89)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (6)

where a is the equatorial radius of the used EGM, r is the geocentric radius of the computation point, GM is the adopted gravitational constant, the coefficients [C.sub.nm] are the fully normalised harmonic coefficients of the disturbing potential and [Y.sub.nm] are the fully normalized spherical harmonics of degree n and order m (cf. Heiskanen, Moritz 1967: 31).

Since the low degree gravity field (up to degree M) is removed from the Stokes integration, then these effects are compensated (i.e. "restored") by the second part of Eq. (1). The latter is nothing but the 'pure' long wavelength contribution of the geoidal height, cf. Heiskanen, Moritz (1967: Sec. 2-17). This method is commonly called a remove-compute-restore (r-c-r) technique and is frequently used in practical geoid computations nowadays.

Modified Stokes formula has also been used in earlier regional geoid modelings in the Baltic Sea region. With a few exceptions the computations were carried out by the Nordic Geodetic Commission (NKG) methodology and software (see Section 3). It is of interest to compare the NKG approach with the generalized geoid determination method by Eq. (1). A NKG geoid model is obtained as quasigeoidal heights (i.e. models of height anomalies), which thereafter can be converted into a geoid model (e.g., by a simplified approach in Heiskanen, Moritz 1967: 327). Indeed, the NKG computational methodology has many similarities (though not exactly the same, see the discussion in Sjoberg, Egren (2002)) with Vincent, Marsh (1974) modification scheme.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (7)

Note that the non-modified Stokes function is utilized here in conjunction with the residual gravity anomaly in the integral. This modification method follows implicitly the generalised Vanicek, Sjoberg (1991) scheme, Eq. (1), with parameters [s.sub.k] = 0 in the integral kernel, cf. Eq. (2). Note that the largest geoid error usually occurs due to the contribution of first degrees (i.e. possible systematic biases between different gravity datums) of the terrestrial data. Vanicek, Featherstone (1998) found that that the unmodified kernel S([psi]) allows low-frequency terrestrial gravity data errors to pass, almost undiminished, into geoid models. The [S.sup.L]([psi]) integration kernels attenuate these errors to a larger extent, but not completely, however. Consequently, even when the residual anomalies are employed in the Stokes integral it is impossible to eliminate (at least in this way), the long wavelength errors in the geoid solution. Also a study by Agren, Sjoberg (2004) concludes that the simple r-c-r scheme with unmodified kernel S([psi]) is sensitive to long-wavelength errors in gravity anomalies. A rather explicit review on the problems of the traditional r-c-r schemes can be found in Sjoberg, Agren (2002) and Sjoberg (2005).

2.2. Complete anomaly versus residual anomaly

The concept of the geoid determination by the r-c-r technique implies that low-frequency gravity signals are removed from the Stokes integration. Importantly, as shown in Sjoberg, Hunegnaw (2000: Eq. 3), the general geoid estimator Eq. (1) can be expressed such that the complete (i.e. the low-degree n [less than or equal to] M harmonics of [DELTA][[uu].sub.n] are included) gravity anomaly instead the residual anomaly is exploited in the integral. According to Sjoberg, Hunegnaw (2000) the geoid estimator Eq. (1) is theoretically equivalent to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], (8)

where the truncation coefficients [Q.sup.L.sub.n] are calculated as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. (9)

Note that the [s.sub.k] coefficient set is the same as used in [S.sup.L]([psi]), the coefficients [Q.sub.n] and [e.sub.nk] are functions of the integration cap radius. They are usually computed using some recursive algorithms.

Comparing Eq. (1) to Eq. (8) we see no particular advantage of reducing [DELTA][uu] in Eq. (1) to Eq. (8) that uses the complete gravity anomaly. One would intuitively expect that any numerical error in the integration becomes smaller due to the use of reduced gravity anomalies. However, studies by Sjoberg, Agren (2002) and Sjoberg (2005) demonstrate that the r-c-r result is as sensitive to various biases as is the case when Stokes's formula is used with complete anomaly as the integral argument.

Accordingly, it can be shown that the NKG modification method, Eq. (7), is theoretically identical to the following modification method using the complete gravity anomaly (cf. Sjoberg 2005).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. (10)

Shortly summarizing--a large variety of different modification schemes exist. Even though the theoretical principles of modifications of Stokes's formula could be rather straightforward, one still needs to make a number of decisions while setting up the actual computations. At this, obviously, the main constraints are the quality and resolution of the available gravimetric data. As already noted, this contribution focuses on geoid computation problems which may occur over the Baltic countries. It is thus appropriate to review earlier geoid modeling efforts over the same target area.

3. Earlier geoid modeling works in the Baltic Sea area

During the last two decades the geoid determination for the whole Nordic region has been carried out within the framework of the Nordic Geodetic Commission. Several NKG geoid models were delivered, see e.g. reference list in (Forsberg 2001). In 1990-ies the NKG geoid models were extended to the Baltic countries. Access to new gravity data from formerly classified sources, and release of EGM96, resulted in achieving better than a dm-accuracy for the regional NKG96 geoid model (Forsberg et al. 1997). After the Baltic Sea Airborne Gravity Campaign 1999 another geoid model (Forsberg 2001) was computed using exactly the same computational setup as for the NKG96 model. The national geoid solutions, either for Estonia, Latvia or Lithuania, were published (Vermeer 1994; Kaminskis, Forsberg 1997; Forsberg 1998; Jurgenson 2001, 2003; Ellmann 2001, 2002, 2004). The Baltic countries are also covered by the European gravimetric (quasi)geoid project (i. e. the EGG97 geoid model and its successor EGG07 (Denker, Torge 1998).

With a few exceptions (e.g. Vermeer 1994; Denker, Torge 1998; Ellmann 2001, 2002, 2004, 2005) the computations were carried out by the NKG methodology and software, whereas EGM96 is utilized in the 1997-2003 computations. It should be noted that the most recent (NKG-2004) model, Forsberg et al. (2004), was also using a GRACE based EGM. The accuracy of the models is evaluated to be 3-7 cm (with a spatial resolution of 3 km). For more details see the original references. To our present knowledge no newer individual geoid solutions either for Estonia, Latvia or Lithuania have been computed from 2004 onwards. Note that the used geoid determination methods are very different, so are the used EGM-s and last but not least, the terrestrial data are incomparable (for more details see Ellmann 2002, 2005; Jurgenson 2003). Unfortunately, the latest geoid models in our region originate in 2004, when the GRACE data processing was only at the initial stage.

Thus, considering the current gravimetric data situation in the Baltic countries and their vicinity, it is feasible to employ such a modification method which performs the best in the case of limited data.

4. Deterministic and stochastic modification methods

4.1. General review

Let us review briefly the main error sources in the gravimetric geoid determination and ways of mitigating their undesired effect.

Recall that the coefficients of geopotential models are obtained from satellite tracking data, while the higher degrees may be obtained by combining the satellite data with terrestrial gravity information. Both datasets contain noise, which unavoidably propagates into the computed geoid undulations. One should also consider the erroneous terrestrial gravity data within the integration domain. Note that the gravity data for integration are usually presented as surface blocks. An additional error occurs thus due to loss of short wavelength gravity information (called discretization error) when estimating the mean (gridded) anomalies [DELTA][uu] from point gravity data.

Minimisation of the geoid estimator errors is the main objective of the modification procedure. The modification methods proposed in geodetic literature can be divided into two distinct classes: deterministic and stochastic approaches. The deterministic approaches principally aim at reducing the effect of the neglected remote zone ([sigma] - [[sigma].sub.0]) making use of a set of low-degree geopotential coefficients. No attempt is made to reduce the errors of the geopotential coefficients and terrestrial data, although the errors of both datasets are contributing to the total error budget. In other words, the modification parameters sn of the deterministic methods are invariant to the two error sources. The most prominent deterministic approaches are Molodenskii et al. (1960), Wong, Gore (1969), Meissl (1971), Vincent, Marsh (1974), Heck, Gruninger (1987), Vanicek, Kleusberg (1987), Vanicek, Sjoberg (1991).

In contrast, the least squares modification methods proposed by Sjoberg (1984, 1991, 2003) allow minimization of the truncation bias, the influence of erroneous gravity data and geopotential coefficients in the least squares (hereafter referred to as LS) sense. Basically, for minimising the errors in geoid modeling the stochastic methods aim at an optimal combination of the data sources (and their error estimates) by adopting a priori or empirical stochastic models.

4.2. A numerical evaluation of modification methods

A study by Ellmann (2004) assessed numerically six different modification methods (three deterministic and three stochastic) by computing a number of high-resolution geoid models over the Baltic countries. The deterministic methods used in this study are Wong, Gore (1969), Vincent, Marsh (1974), and Vanicek, Kleusberg (1987). The stochastic methods are biased, unbiased and optimum least squares modifications by Sjoberg (1984, 1991, 2003). The standard setup for the deterministic methods is to use the residual gravity anomaly in the integral, see Section 2.1. For the sake of comparison these deterministic methods are expressed such that the complete surface gravity anomaly, instead of the residual anomaly, is used in Stokes's integral, cf. Section 2.2. Five methods utilise the modified Stokes function SL([psi]) in Eq. (8), whereas the original (non-modified) Stokes function is used by one modification method (Vincent, Marsh 1974). Accordingly, the simple modification scheme by Eq. (10) is adopted in the Ellmann (2004) study as a very rough representation of the NKG approach. The principles and results of the Ellmann (2004) study are reviewed in the sequel.

4.2.1. Selection of geoid modeling parameters

The Ellmann (2004) experiment employed the EGM96 and the two first GRACE models GGM01s and GGM01c (Tapley et al. 2004) in the practical computations. It should be noted that in the geodetic literature the EGM contribution in the modified Stokes formula is often referred to as "reference field" or "reference model". This may yield a (misleading!) interpretation that the used global geopotential model is errorless. However, a caution is needed with this kind of "wishful thinking". Any EGM should be treated as an ordinary data-source, which thus unavoidably contains errors.

Also the choice of the limit M in Eq. (8) is directly related to the quality of a EGM to be used. In practice, due to restricted access to terrestrial data (or/and the necessity to increase the computational efficiency) the integration cap is often limited to a few hundred kilometers. This implies that a relatively high M should countermeasure this limitation. On the other hand, the EGM error grows with increasing degree, which provides a rationale for a compromise value of M. Traditionally, due to poor accuracy of early EGM-s a rather small modification degree was favored in the computations of many geoid models in the past (e.g. M = 20, see e.g. Vanicek, Kleusberg 1987). Consequently, in these computations even the intermediate EGM wavelength information (suspected to be too erroneous) is avoided.

Recall that a high-resolution EGM is determined from a combination of satellite data and terrestrial gravity data. This combination implies that the two datasets could be correlated in Eqs. (8) and (10). If this correlation appears (e.g. in the case of the full expansion of EGM96 the same data may be used twice in the geoid computation), rigorously, this feature could be accounted for by adding some corrective terms into both formulae. However, it is desirable excluding these cumbersome terms (e.g., see Sjoberg 1991) by selecting appropriate modification limits (more details are given below). If one utilises the "satellite-only" harmonics, this correlation is completely prevented, of course.

Importantly, the space technology advancements have significantly improved the accuracy of recent EGM-s, which allows the user to safely increase the modification degree up to 100 or even beyond. Note that the GGM01s is a "satellite-only" geopotential model developed up to degree/order 120/120, whereas the GGM01c is a combined geopotential model developed to degree/ order 200/200. Unlike the past geopotential models, the GGM01 is highly accurate and homogeneous also for the intermediate spectrum. For the final geoid model it was aimed at selecting the modification limit M such that the correlation between the EGM-derived and terrestrial datasets is entirely prevented.

Further. It was commonly believed that a good quality terrestrial data-set may somewhat compensate the shortages of the past geopotential models. This provides a rationale to exploit vast data areas, wherever possible in the past geoid computations. This goal is also aimed at the NKG geoid modeling works. In practical NKG computations the complete ([n.sub.max] = 360) expansion of the used EGM models in conjunction with a very large integration domain (where available) was often utilised for regional geoid modeling. Hereafter we refer to it as the simple modification method. In this respect it is interesting to review the results of Vanicek, Featherstone (1998). They concluded that even the error-free gravity data, when used in a limited spatial domain, can never completely correct the errors of geopotential models.

Let us focus on some aspects of the simple modification method. Assume for a moment that available terrestrial data is completely errorless. In such a case, according to a study by Sjoberg, Agren (2002), in order to reduce the truncation bias in the NKG approach to a cm level, the integration cap must exceed 10[degrees] (approx. 1100 km around each computation point).

In the context of the Baltic geoid modeling this means that good quality and dense terrestrial gravity data is needed even beyond the Moscow meridian. Recall, however, that the Baltic countries are located at the eastern edge of the European Gravimetric Geoid (EGG) and NKG geoid modeling projects. Even though some international gravity data-bases possess indirect gravity field information (i.e. gridded anomaly values) over eastern part of Russia and Belarus, but access to the pointdata records is still restricted. For instance, these gridded data were not made available for the Ellmann (2004) study, see Fig. 1. Further, to our present knowledge these data can be obsolete and are not verified against the possible presence of systematic biases with respect to modern (absolute gravity based) data-sets. In other words, over this significant area the quality of the gravity data remains largely unknown. Therefore, the geoid determination in the Baltic countries may suffer in the possible data shortages to the east from the target area.

[FIGURE 1 OMITTED]

One should admit, however, that the terrestrial data is never without errors and these errors will inevitably be present in the geoid estimation process by the simple modification method. Ellmann (2004) shows the larger the integration radius the larger the influence of terrestrial errors to the geoid modeling. Hence, conversely to the earlier assumption, the integration cap size should be in reasonable balance when applying the simple modification method.

It should be noted, that the limited extension of terrestrial data is the most serious restriction for the Ellmann (2004) study. Due to the limited data access the integration cap radius was chosen 2[degrees] (corresponding to 220 km) around each computation point. Note that the upper modification limit L in Eq. (4) is arbitrary and generally it is not equal to the series expansion M in Eq. (1). For instance, the choice of the upper limit L may be related to the integration cap radius [[psi].sub.0]. According to de Witte (1966) the truncation bias tends to minimum when the integration cap radius is extended to the zerocrossings of integration kernel (modified or non-modified Stokes function, i.e. [S.sup.L]([psi]) or S([psi])). Instead, Heck and Gruninger (1987) propose to place a constraint on the values that can be chosen for the modification, i.e. either on the parameter L or [[psi].sub.0], such that [S.sup.L]([psi]) ceases to zero at the edge of the integration cap. Accordingly, we want the modified Stokes function SL([psi]) to become zero at the edge of the integration cap. Therefore the choice [[psi].sub.0] = 2[degrees] is the basis for determining the L for the WongGore and Vanicek-Kleusberg methods. For the VanicekKleusberg modification method this condition is fulfilled with L = 67. Note also, that the limit M should preferably not exceed the "satellite-only" harmonics (i.e. M [less than or equal to] 95) of GGM01s. Moreover, since the immediate area of interest is gravimetrically well surveyed (see Fig. 1), the terrestrial data are probably a better source for the medium and short wavelength geoid information. This suggests a smaller value for M. We select L = M = 67 for our computations. In other words, the modification degree of the kernel SL([psi]) is the same as the upper limit M of the geopotential harmonics to be used. This choice is also supported by a circumstance that the error degree variances of the GGM01s harmonics for n [less than or equal to] 67 are smaller than the degree variances with n [less than or equal to] 20 (often preferred for regional geoid computations in the past) of any previous geopotential model. The kernel [S.sup.L]([psi]) with the WongGore coefficients, 2/(n - 1), becomes zero (at [[psi].sub.0] = 2[degrees]) when L = 31 in Eq. (2). For the sake of the comparison, we choose M = 67 for both deterministic approaches.

It should be noted that a typical feature of the LS parameters sn is that the LS parameters "force" [S.sup.L]([psi]) to almost zero at the edge of any pre-selected integration cap. For consistency of comparison we select the same limits also for the LS modifications as for the VanicekKleusberg deterministic method, i.e. L = M = 67. This choice allows us to take full advantage of the "satelliteonly" GGM01s, whereas there is strictly no correlation with terrestrial data. On the other hand, if the terrestrial gravity data in the area of interest is poor, there is no reason to abandon erroneous high-degree harmonics of the used EGM. Remember, that the gain from high-degree harmonics may be more rewarding than possible damage. A relevant matter is how to reduce the errors and find a correct balance (e.g. weights) between different data sources (i.e. EGM and terrestrial gravity anomaly). The LS modifications are designed for aiming at this goal. Besides, Ellmann (2001) concluded, that the LS procedure is able to adjust the data in a way that the geoid model becomes rather insensitive to the maximum degree of modification, because it matches the different types of data in an optimum way. Conversely, the limit L should be selected very carefully for the deterministic methods. Hence, one needs to consider many (often rather contro-versial) arguments when choosing the limits L and M.

Summarising, the limits L = M = 67 will be used everywhere in the computations for the Vanicek-Kleusberg and LS modification methods, whereas M = 67 and L = 31 are adopted for the Wong-Gore method. For the simple modification method see a note below.

4.2.2. The results

Different modification methods yielded corresponding gravimetric geoid models for the Baltic countries.

The quality of these models was assessed from the comparisons with GPS-levelling data. In particular, three sets (one from each country) of GPS-levelling points were used for an independent evaluation of the computed geoid models. Four transformation parameters between the gravimetric geoid models and the GPS-levelling data were defined and thereafter a polynomial fit (Ellmann 2004; Eq. 23) was applied. The achieved accuracy was more or less the same for modification methods with SL([psi]), but formally the unbiased LS modification method yielded the best statistics for the post-fit residuals. The corresponding geoid model is referred to as BALTgeoid-04, see Fig. 2. The best RMS error of the GPS-levelling post-fit residuals were as follows: 5.3 cm for the joint Baltic geoid model and 2.8, 5.6 and 4.2 cm for Estonia, Latvia and Lithuania, respectively. It seems that the accuracy of the tested modification methods (with [S.sup.L]([psi])) is at least the same level as the accuracy of the used control points.

The Ellmann (2004) study tested the simple modification {see Eq. (10), equivalent to Vincent, Marsh (1974) modification} in conjunction with three geopotential models--GGM01s, GGM01c and EGM96. The computation results are validated by the same sets of GPS-levelling points. Our attempts to exploit M = 67 and M = 120 (the latter is the maximum degree of GGM01s) and M = 200 (maximum degree of GGM01c) did not provide satisfactory results, i.e. the statistics of the models utilising Eq. (10) were worse (generally, a smaller M yielded the larger RMS error) than the ones for the five methods with [S.sup.L]([psi]). As it was already discussed, this could be mainly due to large truncation bias associated with the simple modification method, whereas in the methods utilising SL([psi]) this bias is efficiently mitigated. [Perhaps this problem for the methods with S([psi]) is less critical when data from more extended areas would be used in the integration.] Only the complete expansion of EGM96 (i.e. M = 360) provides comparable post-fit ac curacies with other results. It can be concluded thus, that the GRACE-only models do not qualify for the precise modeling by Eq. (10). Possibly, even future geopotential models from GOCE data alone (resolution up to degree 270 is expected, www.esa.int) also would not be satisfactory for the simple modification method.

[FIGURE 2 OMITTED]

Even though the results of the Ellmann (2004) study are generally superior to the NKG joint geoid models, there are several circumstances preventing direct comparison of the results. These are:

(i) The fast Fourier transform (FFT) method is utilised in the computations of the NKG geoid models. Most often all the data-points from rather large (rectangular-shaped) region are involved for computing each geoidal height. In contrast, in the Ellmann (2004) study the data is strictly limited to a spherical cap (with a pre-defined radius) around each computation point.

(ii) The NKG software uses residual gravity anomalies for the integration, whereas the Ellmann (2004) study used the complete gravity anomaly

(iii) Regional NKG solutions (such as Forsberg 2001, Forsberg et al. 2004) include a large Russian dataset, which was not available for the Ellmann (2004) study.

(iv) The NKG computations yield height anomalies, whereas the outcomes of the Ellmann (2004) study are geoidal heights.

Apparently, an appropriate geopotential model is also essential to determine the regional gravimetric geoid model accurately. It should be noted that the numerical discrepancies between various geoid models using SL([psi]) and GGM01s remain within [+ or -] 9 cm. At the same time the differences between the GGM01s and EGM96 based deterministic geoid models range in the target area from -6 to +17 cm. This implies that the deviations between the contemporary geopotential models are more crucial than the differences between the tested modification methods. It was also detected that the GRACE-based based regional geoid models agree with the 1977 Baltic Height System better than the EGM96 based geoid models. Note that Ellmann et al. (2009) results indicate that the EGM08 based regional models may possess a great potential for further refinement of the regional geoid models. As expected, the discrepancies between the upcoming geopotential models will be reasonably small, which prompts for more careful selection of geoid modeling methods and parameters. It is expected thus that the modifications of Stokes's will not lose its actuality for the years to come. Preferably, several different modification methods should be simultaneously tested at regional geoid determination in similar manner as in Ellmann (2004).

4.2.3. Further reading

It has been almost five years since 2004, when the aforementioned results were presented and defended in PhD Dissertation at the Royal Institute of Technology (KTH), Stockholm. The Dissertation was titled as "T/re geoid for the Baltic countries determined by the least squares modification of Stokes' formula" and it was published as a KTH internal publication 04:013 of Trita-INFRA series (ISSN 1651-0216), see also KTH's Geodesy Report 1000 Series, Dissertation No 1061. The main quintessence of the study was later published in Journal of Geodesy (Ell mann 2005).

Undoubtedly, over the past years geoid modeling methods have gradually developed. The advances of the satellite technology have already been mentioned in the Introduction. However, the detected geoid modeling problems in the Baltic countries have largely remained the same (for more details see Section 5).

Recently, the author was given an opportunity to reiterate the matters in question in a separate book. The title of this new book is Ellmann A. (2009): Modified Stokes's Formula for Regional Geoid Modeling: Deterministic and Stochastic Modifications of Stokes's Formula for Computing an Improved Geoid Model over the Baltic Countries, which was published by the VDM Verlag. This book (ISBN 9783639128192) can be accessed/ordered via online bookstores (such as Amazon.com) and libraries.

The overall structure and the core of the new book is intentionally left the same as the original PhD Dissertation. Minor mistakes have been corrected and also changes in wording of the text have been made. Accordingly, this book can be considered as a new, updated and corrected edition of the PhD Dissertation. The main intention of this book is to present comprehensive guidelines for the application of different modification methods that can be utilized in any given region worldwide. Throughout the book, it is assumed that the reader already knows the basics of physical geodesy. This book can be used as complementary text to graduate level courses in the discipline of physical geodesy. Researchers primarily interested in regional geoid modeling may also be interested in this book.

Besides the obvious theoretical challenges of geoid modeling (discussed explicitly in Ellmann 2004 and 2009) is affected by many practical issues. For instance, there are evidences that the quality of the terrestrial gravity data may distort the geoid modeling outcomes in the Baltic Sea region. Recall, that most of the data within the land masses of the Baltic countries have been collected before the 1990-ies. Conversely, the modern gravity networks are established decades after the historic gravity surveys. A practical case study is reviewed in next section, which validates Estonian terrestrial gravity data.

5. Terrestrial gravity data evaluation in Estonia

Clearly, perfection of any geoid modeling method is diminished or even meaningless with insufficient data quality and coverage. Therefore a great deal of attention should be paid to the reconciliation of the terrestrial gravity data. The treatment of the data collected with different methods and equipment, during several decades by different nations and specifications, requires careful study before their use in the geoid computation. Therefore, all the undesired systematic biases need to be detected and eliminated, followed by the conversion into the common gravimetric datum.

Ellmann et al. (in press, 2009) study focuses on the quality assessment of Estonian terrestrial gravity data, which have been collected by different institutions over many decades. The oldest gravity survey points were originally based on the 1930 realization of the international Potsdam gravity system. In the 1950-ies these (along with newer gravity data) were converted into a new (1955) realization of the Potsdam system in Estonia. In 1970-ies the worldwide gravity system IGSN71 was introduced also in Estonia. Further, the contemporary Estonian gravity system is currently based on a nationwide set of absolute gravity measurements. The gravity network is being developed and maintained by the Estonian Land Board (ELB). Attempts have been made to convert the historic gravity survey results to the current gravity system.

Accordingly, Ellmann et al. (in print, 2009) investigated the links between the contemporary gravity system and the following two datasets: (i) the 1949-1958 gravity survey by the Institute of Geology of the Estonian Academy of Sciences; (ii) the 1967-2007 gravity surveys of the Geological Survey of Estonia (GSE).

The ELB gravity network and recent survey points are used as control points in this study. The agreement of the datasets (i) and (ii) with the control points (altogether 424 points) was determined empirically by using the gravity survey results for predicting the simple Bouguer anomalies at the locations of the control points.

The 1949-58 survey consists of 4000+ gravity data points. The detected discrepancies between the observed and predicted Bouguer anomalies of the 1949-1958 gravity survey at the locations of the control points range from -4.5 to +3.8 mGal, with the RMS error of the discrepancies as of 1.38 mGal. The mean of the detected discrepancies is -0.53 mGal. Unfortunately, the discrepancies between the two datasets are not random at all. The largest systematic discrepancies can be observed in South-Estonia, where an average bias can approach up to 3-4 mGal. This is clearly inadmissible for precise geoid modeling. In other parts of Estonia, the discrepancies are less, but still worrisomely exceeding 1 mGal.

It should be noted that the 1949-1958 gravity survey is the only available set of gravity data covering the whole Estonia with suitable density for geoid modeling. Therefore, this data-set has been employed in all the earlier Estonian and NKG and EGG geoid modeling works.

The main purpose of the GSE gravity surveys was the geological mapping of the crystalline basement. The interval between the tracks was 1 km in average, whereas the along-track station separation was 250 m. At present the gravity database of the GSE consists of 126 609 survey points, which were used in comparisons. The detected discrepancies between the observed and predicted Bouguer anomalies at the locations of the control points range from -1.9 to +1.9 mGal, with a mean of -0.06 mGal. The RMS error of the discrepancies is 0.33 mGal. The detected discrepancies are more or less random, meaning that a reasonable agreement between the control values and the GSE gravity survey results was achieved.

The emphasis of this study was given for assessing the suitability of the existing gravity data to ensure a 1 cm geoid modeling accuracy over Estonia and its surroundings. The accuracy of the GSE gravity data seems to meet this requirement. Unfortunately, the GSE surveys are not covering the whole of Estonia.

Naturally, within the GSE survey areas the 1949-58 gravity survey data can be down-weighted or simply removed from the geoid computations. Over the rest of Estonia, however, the usage of the 1949-58 survey results seem still to be unavoidable.

An alternative to the 1949-58 survey data would be using gravity values from global geopotential models. For instance, the spatial resolution, 9 km, of the EGM08 (Pavlis et al. 2008) is quite comparable with the density of the 1949-58 survey points. The performance of the EGM08 model over the Baltic countries was evaluated by Ellmann et al. (2009). One of the main conclusions of their study was that the EGM08 derived gravity quantities agree reasonably well with the terrestrial survey data in Estonia. Apparently, the 1949-58 gravity survey data have entirely been utilised in the compilation of the EGM08. This data set has been made available for several international gravity databases since the 1990-ies. Conversely, the modern gravity network points and the results of new surveys were not accessible at the EGM08 compilation. For detecting the discrepancies between the contemporary gravity datum, and the EGM08 derived gravity field the free-air anomalies were computed at the locations of the control points (the same as used in the Ellmann et al. (2009) study) in Estonia. The statistics of the detected discrepancies between the newly measured and EGM08-derived gravity anomalies resembles the discrepancies between the control points and the 194958 survey data. Also the distribution of the discrepancies is rather similar. In other words, the 1949-58 survey data errors have been absorbed into the EGM08 highfrequency spectrum. Therefore the use of the EGM08 model cannot provide better results, than the 1949 - 58 survey data. Another implication is that such a distorted EGM will distort the resulting regional geoid model, if the modification limits are not selected properly.

All in all, it seems that for accurate geoid modeling the 1949-58 gravity survey results need to be replaced by new field measurements. Certainly, this is a quite burdensome task, requiring lots of effort and well coordinated actions. However, this is needed for the sake of the consistency of the national gravity datasets. For this the primary attention should be paid to the most critical regions, which are outlined in Ellmann et al. (in print, 2009). Additionally, the gravity field over major water bodies, such as the Peipsi lake and the Vortsjarv lake, also the Riga Gulf, need to be specified. Our further studies will be devoted for achieving this goal in Estonia. The detected discrepancies will be studied and ultimately resolved in future gravity field and geoid modeling works.

Most likely the problems revealed in Ellmann et al. (2009) can also be identified in Latvia and Lithuania as well. Thus it would be very useful to join (or at least coordinate) efforts to eliminate terrestrial data shortages in our region.

6. Summary and conclusions

This study reviewed efforts for geoid modeling over the Baltic countries with an emphasis on the comparison of the practical outcomes of different modifications of Stokes's formula. In general, all the modification methods combining the modified Stokes function SL([psi]) and an intermediate expansion of the geopotential models provide a reasonable accuracy even for rather small cap sizes. Conversely, to ensure cm accuracy in the geoid modeling with unmodified S([psi]), the employment of a high degree expansion of an EGM and data from a large integration domain are necessary. This yields that the "satellite-only" geopotential models in conjunction with unmodified S([psi]) are not suitable for computations of high-resolution regional geoid modeling.

It is also concluded that the methods using SL([psi]) mitigate the terrestrial data errors more efficiently than the modification methods using unmodified Stokes's function S([psi]) in the integral.

Recent technological and theoretical advances have created preconditions to achieve 1-cm accuracy for geoid model. However, in order to achieve a high resolution for the regional geoid models the terrestrial data should not be totally abandoned even in the upcoming GOCE era. Disappointingly, the deficiencies in historical terrestrial gravity data may corrupt this objective. For in many cases the shortages of the gravity data cannot be completely removed, say by mathematical correction methods.

Consequently, a careful and versatile reconciliation of all available geodetic data in the Baltic Sea region is desired. This can only be done within the frame of international cooperation, desirably resulting in modernization of the gravity databases and vertical networks of the Baltic countries. Also the need for the new gravity datacollection has become quite obvious.

But also other geoid modeling related subjects need to be properly treated. Recall that geoid determination by Stokes's formula holds rigorously only on a spherical boundary, assuming also the masses outside the geoid to be absent. Consequently, the original and modified Stokes formulae should also comprise some correction terms accounting for the Earth's ellipticity and the contribution of topographic and atmospheric masses. For instance in the Ellmann (2004) study all these effects are accounted for by means of additive corrections. However, all these geoid modeling related subjects need to be further investigated in the context of the Baltic Sea region.

The first GOCE results are expected to be released in the beginning of 2010. Therefore it is now a good opportunity to revise theoretical and practical aspects for future geoid modeling works in the Baltic countries.

doi: 10.3846/gc.2010.01

Acknowledgements

This study is a part of activities undertaken within the framework of the Estonian Science Foundation grant ETF7356 "Application of space technologies to improve geoid and gravity field models".

Received -6 10 2009, accepted 05 10 2010

References

Agren, J.; Sjoberg, L. E. 2004. Comparison of some methods for modifying Stokes' formula in the GOCE era, in Proceedings CD of 2nd International GOCE user workshop "GOCE, The Geoid and Oceanograpy", hold at ESA-ESRIN, Frascati, Italy, March 08-10, 2004.

de Witte, L. 1966. Truncation errors in the Stokes and Vening Meinesz formulae for different order spherical harmonic gravity terms. Aerospace rep TR-669 (SG 230-35)-5.

Denker, H.; Torge, W. 1998. The European gravimetric quasigeoid EGG97--an IAG supported continental enterprise, in Forsberg, R.; Feissel, M.; Dietrich, R. (Eds.). Geodesy on the Move. IAG Symp. Series vol. 119, Springer Verlag, 249-254.

Denker, H., et al. 2007. The Development of the European Gravimetric Geoid Model EGG07 Observing our Changing Earth, in Proceedings of the 2007 IAG General Assembly, Perugia, Italy, 2-13 July, 2007, International Association of Geodesy Symposia, Vol. 133 (M. Sideris, Ed.), 177-185, Springer.

Ellmann, A. 2001. Least squares modification of Stokes formula with applications to the Estonian geoid. Royal Institute of Technology, Division of Geodesy, Report No. 1056.

Ellmann, A. 2002. An improved gravity anomaly grid and a geoid model for Estonia, in Proceedings of the Estonian Academy of Sciences. Geology 51: 199-214.

Ellmann, A. 2004. Te geoid for the Baltic countries determined by the least squares modification of Stokes' formula. Royal Institute of Technology, Department of Infrastructure, Geodesy Report No. 1061, Stockholm.

Ellmann, A. 2005. Two deterministic and three stochastic modifications of Stokes's formula: a case study for the Baltic countries, Journal of Geodesy 79: 11-23. doi:10.1007/s00190-005-0438-1

Ellmann, A. 2009. Modified Stokes's Formula for Regional Geoid Modeling: Deterministic and Stochastic Modifications of Stokes's Formula for Computing an Improved Geoid Model over the Baltic Countries, VDM Verlag.

Ellmann, A.; Jurgenson, H. 2008. Evaluation of a GRACE-based combined geopotential model over the Baltic countries, Geodesy and Cartography 34: 35-44. doi:10.3846/1392-1541.2008.34.35-44

Ellmann, A.; Kaminskis, J.; Parseliunas, E.; Jurgenson, H.; Oja, T. 2009. Evaluation results of the Earth Gravitational Model EGM08 over the Baltic countries, Newton's Bulletin 4: 110-121.

Ellmann, A.; All, T.; Oja, T. (in print, 2009). Toward unification of terrestial gravity data sets in Estonia, Estonian Journal of Earth Sciences 58.

Forsberg, R. 1998. Geoid tailoring to GPS--with example of a 1-cm geoid of Denmark, in Vermeer, M.; Adam, J. (Eds.). Second Continental Workshop on the Geoid in Europe. Budapest, 10-14 March, 1998. Reports of the Finnish Geodetic Institute 98: 4, 191-197.

Forsberg, R. 2001. Development of a Nordic cm-geoid with basics of geoid determination, in Harsson, B. G. (Ed.). Nordic Geodesy towards the 21st century. Lecture notes for Autumn School, Nordic Geodetic Commission, 67-88, Statens kartverk, Geodetic Publication 2001: 1.

Forsberg, R.; Kaminskis, J.; Solheim, D. 1997. Geoid for the Nordic and Baltic Region from Gravimetry and Satellite Altimetry, in Segawa, J.; Fujimoto, H. and Okubo, S. (Eds.). Gravity, Geoid and Marine Geodesy. IAG Symposium Series volume 117, 540-547, Springer Verlag.

Forsberg, R.; Strykowski, G.; Solheim, D. 2004. NKG-2004 ge oid of the Nordic and Baltic area, Poster submitted to the unrefereed proceedings on CD-ROM from the International Association of Geodesy Conference "Gravity, Geoid and Satellite Gravity Missions", Aug 30-Sep 3, 2004, Porto, Portugal.

Forste, C., et al. 2006. A mean global gravity field model from the combination of satellite mission and altimetry/gravimetry surface data: EIGEN-GL04c, Geophysical Research Abstracts 8: 03462.

Heck, B.; Gruninger, W. 1987. Modification of Stokes's integral formula by combining two classical approaches, in Proceedings of the XIX General Assembly of the International Union of Geodesy and Geophysics, Vol. 2. Vancouver, Canada, 309-337.

Heiskanen, W. A.; Moritz, H. 1967. Physical Geodesy. W. H. Freeman and Company. San Francisco.

Jurgenson, H. 2001. Gravity measurements and calculations of the geoid model for Estonia, in Ingemann, V. and Madsen, F. (Eds.). Proceedings from Seminar of the Danish & Baltic sector Programme, Riga March 2-3, 1999, 94-110, Copenhagen.

Jurgenson, H. 2003. Determination of Estonian Precise Geoid. PhD Thesis, Estonian Agricultural University, Tartu.

Kaminskis, J.; Forsberg, R. 1997. A new detailed geoid for Latvia, in Segawa, J.; Fujimoto, H.; Okubo, S. (Eds.). Gravity, Geoid and Marine geodesy. IAG Symposium Series, volume 117, 621-628, Springer Verlag.

Lemoine, F.G., et al. 1998. The development of the joint NASA, GSFC and NIMA Geopotential Model EGM96. NASA/TP 1998-206861, NASA Goddard Space Flight Center, Greenbelt Maryland, USA.

Meissl, P. 1971. Preparation for the numerical evaluation of second order Molodenskii- type formulas. Dept. of Geodetic Science, OSU Report No. 163, Columbus, Ohio.

Molodenskii, M. S.; Eremeev, V. F.; Yurkina, M. I. 1960. Methods for study of the external gravitational field and figure of the Earth (translated from Russian in 1962). Israel program for scientific translations, Jerusalem.

Pavlis, N. K.; Holmes, S. A.; Kenyon, S. C.; Factor, J. K. 2008. An Earth Gravitational Model to Degree 2160: EGM2008. Presented at the 2008 General Assembly of the European Geosciences Union, Vienna, Austria, April 13-18, 2008.

Sjoberg, L. E. 1984. Least squares modification of Stokes and Vening-Meinesz formulas by accounting for errors of truncation, potential coefficients and gravity data. The Dept. of Geodesy report No. 27, University of Uppsala.

Sjoberg, L. E. 1991. Refined least-squares modification of Stokes formula, Manuscripta Geodaetica 16: 367-375.

Sjoberg, L. E. 2003. A general model of modifying Stokes' formula and its least-squares solution, Journal of Geodesy 77: 459-464. doi:10.1007/s00190-003-0346-1

Sjoberg, L. E. 2005. A discussion on the approximations made in the practical implementation of the remove-computerestore technique in regional geoid modelling, Journal of Geodesy 78: 645-653. doi:10.1007/s00190-004-0430-1

Sjoberg, L. E.; Hunegnaw, A. 2000. Some modifications of Stokes formula that account for truncation and potential coefficient errors, Journal Geodesy 74: 232-238. doi:10.1007/s001900050281

Sjoberg, L.E.; Agren, J. 2002. Some Problems in the Theory used for the NKG Geoid Model, in Poutanen, M. and Suurmaki, H. (Eds.). Proceedings of the 14th General Meeting of the Nordic Geodetic Commission, 121-129. Kirkonummi/ Helsinki.

Stokes, G.G. 1849. On the variations of gravity on the surface of the Earth, Transactions of the Cambridge Philosophical Society 8: 672-695.

Tapley, B. D.; Bettadpur, S.; Watkins, M. M.; Reigber, Ch. 2004. The gravity recovery and climate experiment: mission overview and early results, Geophysical Research Letters 31, L09607. doi:10.1029/2004GL019920

Vanicek, P.; Kleusberg, A. 1987. The Canadian geoid-Stokesian approach, Manuscripta Geodaetica 12: 86-98.

Vanicek, P.; Sjoberg, L. E. 1991. Reformulation of Stokes's theory for higher than second-degree reference field and modification of integration kernels, Journal of Geophysical Research 96(B4): 6529-6339. doi:10.1029/90JB02782

Vanicek, P.; Featherstone, W. E. 1998. Performance of three types of Stokes's kernel in the combined solution for the geoid, Journal of Geodesy 72: 684-697. doi:10.1007/s001900050209

Vermeer, M. 1994. A fast delivery GPS-gravimetric geoid for Estonia. Reports of Finnish Geodetic Institute, Helsinki.

Vincent, S.; Marsh, J. 1974. Gravimetric global geoid, in Veis, G. (ed.). Proc. of Int. Symp. On the use of artifical satellites for geodesy and geodynamics. National Technical University, Athens, Greece.

Wong, L.; Gore, R. 1969. Accuracy of geoid heights from the modified Stokes kernels, Geophysical Journal Royal Astronomy Society 18: 81-91.

Artu Ellmann

Tallinn University of Technology, Dept of Civil Engineering, Ehitajate rd. 5, 19086 Tallinn, Estonia E-mail: artu.ellmann@ttu.ee

Artu ELLMANN, Prof., Tallinn University of Technology, Dept of Transportation, Ph. +372 6202 603, Fax +372 6202 601.

He holds a PhD (since 2004) degree from the Royal Institute of Technology (KTH) in Stockholm. 2004-2006, a research fellow at the Geodetic Research Laboratory, University of New Brunswick, Canada. Over 30 presentations in international conferences, author/co-author of over 60 publications in referred Journals and Conference Proceedings, technical reports, magazine articles. National correspondent to the International Association of Geodesy (IAG).

Research interests: physical geodesy, gravity field and geoid modelling in particular, national geodetic networks.