Introduction

In the Molodensky method of determination of the figure of the Earth, the quasigeoid is evaluated from the data given on the surface of the Earth. The method is purely based on external gravity field parameters without any hypothesis about density distribution inside topographic masses (Heiskanen and Moritz 1967, chap. 8, Molodensky et al. 1962). Within this approach, the geoid is given by means of numerical addition of a geoid-quasigeoid separation term which can be determined by using the gravity anomalies, i.e. Bouguer as well as free air anomalies with elevation data. The physical and mathematical nature of geoid-quasigeoid separation term has been already investigated with good accuracy by different authors (Pick et al. 1973; Sjöberg 1995; Rapp 1997; Barlik 2000). A more recent contribution has been made by Tenzer et al. (2006). It relies upon the mean gravity disturbance formula and uses normal heights. Nevertheless, the geoid-quasigeoid separation term can be determined with sufficient accuracy through the use of gravity anomalies and orthometric heights. The expression for its determination has the following form, see Sjöberg (1995)

$$ {{\text{N}}_P} - {\zeta _{\text{P}}} \approx \frac{{\Delta {{\text{g}}_B}}}{{\overline \gamma \left( \Omega \right)}}{\text{H}}\left( \Omega \right) + \frac{{{{\text{H}}^2}\left( \Omega \right)}}{{2\overline \gamma }}\left( {\frac{{\partial \Delta {{\text{g}}^F}}}{{\partial {\text{H}}}}} \right) + \left( {{{\text{O}}^3}} \right) $$
(1)

where, Δg B is the complete Bouguer anomaly including the Bouguer plate effect plus terrain correction and \( \overline \gamma \) is the average theoretical gravity, M P the geoid height, H orthometric height, ζP the height anomaly at the corresponding point P at the Earth surface and (O 3) represents sum of the remaining higher order terms which are rather small and negligible. The 1st term on right side of Eq. 1 can be determined easily as all parameters are simple to compute. The 2nd term however, needs special solution for the vertical gradient of free air anomaly, i.e. \( \frac{{\partial \Delta {g^F}}}{{\partial H}} \), which can be estimated from the computation of a singular integral in Eq. 2 of the next section.

In the geodetic literature one finds that calculation of the geoid-quasigeoid separation is based on Eq. 1 derived by Sjöberg (1995) in the Helmert condensation formulation. This formula is based on the assumption that the internal mass distribution fulfills the hypotheses of constant topographic density and isostatic mass compensation according to Helmert's second condensation model. However, the effect of the variable density distribution and the deviation of Helmert's compensation model from reality has been considered to be constant here in this study. Though it has significant contribution which can be treated in a separate work.

Singular integrals, associated with the reciprocal distance, often appear in the different integral computations in physical geodesy, for example, geoidal heights, deflections of vertical, topographic corrections and vertical gravity anomaly gradient computations. These integrals become singular when the distances as well as the height differences between the computation and running points approach to zero. In this case these integrals need special treatment at the computation point and/or in the innermost zone. Numerous treatments were proposed in the past decades; among them include template computations, piecewise polynomials, spline interpolations, fast Fourier transform, finite element method, or a direct approach by means of Newton–Cotes integration formulas, e.g. Bian and Sun (1994), Bian (1997) and Klees and Lehmann (1998).

The approach adopted here for the computation of the singular integral (that yields the vertical gravity anomaly gradient) is based on its numerical solution using Newton–Cotes formulas (Bian 1997) with extended integration area and order (n = 6) of the integration in planar approximation. It should be mentioned, that integral expressions can be written with good accuracy in the planar approximation if these procedures are applied in the innermost area (Bian 1997). Furthermore, it is known fact that integral for the vertical gravity anomaly gradient has a fast decreasing character for the far zone due to an 1/l 3 term as compared to other integrals like e.g. the Stokes and the Vening Meinesz integral, see Heiskanen and Moritz (1967). Therefore, the effect of the innermost area can be considered as enough due to its relatively small contribution to the geoid-quasigeoid separation. In addition, numerical investigations have been made to estimate the effect of different integration radii on the geoid-quasigeoid separation in the innermost zone which is one of the main purposes of this work.

The main contribution to the geoid-quasigeoid separation term originates from the Bouguer effect that includes Bouguer anomaly plus terrain correction as determined from observed gravity data. Other studies include global geopotential model dependent terms as suggested and evaluated by Rapp (1997) and then by Nahavandchi (2002). In this case the gradients of height anomaly are evaluated at earth surface with respect to the earth’s radius and the normal gravity in the form of correction terms, see Rapp (1997) and Nahavandchi (1998), which are then added to the gravity data based geoid-quasigeoid separation. This study will concentrate on the free air anomaly dependent component which has an effect on the order of centimeter in low- to mid-range elevations. Therefore, it requires to be estimated with reasonable accuracy for accurate determination of the total geoid-quasigeoid separation term.

The height datum of Pakistan is based on the orthometric height system. Since Pakistan has a diversity of terrain distribution due to vast expanse of land comprising both plain lands of southern to mid elevation central and northern provinces and then to very high Himalayan mountain ranges. The quantification of the geoid-quasigeoid separation is essential due to the reason that most geodetic boundary value problems give quasigeoid as their final solution with the exception of the pure Helmert condensation. This fact lies in the way of handling topography in these methods, e.g. Molodensky's method with residual terrain modelling (RTM) and combined RTM/Helmert schemes (for details see e.g. Omang and Forsberg 2000). This study focuses on the 2nd term of Eq. 1 dependent on vertical gradient of free air anomaly in the innermost zone. “Brief problem statement and method of solution” describes briefly, the problem and solution of the integral of vertical gravity gradient in planar approximation, “Solutions for the different areas of integration” deals with the numerical solutions for different orders using Newton–Cotes integration, “Data processing methodology and analysis of the results” highlights the analysis of test results and “Conclusions and recommendations” concludes with some recommendations.

Brief problem statement and method of solution

The planar/flat Earth approximation of the vertical gravity anomaly gradient can be expressed as

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{1}{{2\pi }}\iint {_{xy}\frac{{\Delta g\left( {{\rm x}, y} \right) - \Delta {g_0}}}{{{r^3}}}{\rm d} xdy} - \frac{2}{R}\Delta {g^F} $$
(2)

where Δg 0 is the free air gravity anomaly at the computation point and \( r = \sqrt {{{x^2} + {y^2}}} \), with x and y as tangent planar coordinates of the moving point, see Heiskanen and Moritz (1967). The indirect effect (2/R) Δg F term will not be written further and its effect will be included numerically, though it has insignificant numerical value in mid-range elevations. For the gravity anomaly to be used for numerical integration, it is supposed that Δg(x, y) is twice differentiable at the computation point as

$$ \Delta g\left( {x,y} \right) = \Delta {g_0} + {g_x}x + {g_y}y + \frac{1}{{2!}}\left( {{g_{xx}}{x^2} + 2{g_{xy}}xy + {g_{yy}}{y^2}} \right) + \ldots, $$
(3)
$$ {\hbox{where}}\quad {g_x} = \frac{{\partial \Delta g}}{{\partial x}}\;{\hbox{and}}\;{g_y} = \frac{{\partial \Delta g}}{{\partial y}}, $$
(4)

see Heiskanen and Moritz(1967). Then the 1st part of R.H.S. of Eq. 2 can be rearranged as (Bian 1997),

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{1}{{2\pi }}\iint {\frac{{\Delta g\left( {x,y} \right) - \Delta {g_0} - {g_x}x - {g_y}y}}{{{r^3}}}}{\rm d} \sigma, $$
(5)

since \( \iint {\frac{{x{\text{d}}\sigma }}{{{r^3}}}} = 0 \) and \( \iint {\frac{{y{\text{d}}\sigma }}{{{r^3}}}} = 0 \), due to symmetry of the area of integration.

In order to remove the singularity of Eq. 2, Bain and Dong (1991) and Bian and Sun (1994) suggested the rearrangement and change of variable. The integration area in Eq. 5 dq = dxdy so that \( \sigma \equiv \left\{ {\left( {x,y} \right) \in {R^2}; - \infty < x < \infty, - \infty < y < \infty } \right\} \) is divided into two parts \( {\sigma_1} \equiv \left\{ {\left( {x,y} \right) \in {R^2};\left| y \right| < \left| x \right|} \right\} \) and \( {\sigma_2} \equiv \left\{ {\left( {x,y} \right) \in {R^2};\left| x \right| < \left| y \right|} \right\} \). Introducing x = x, y=kx for σ 1 and x = ky and y = y for σ 2, the final form of Eq. 5 becomes

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{1}{{2\pi }}\int\limits_{ - \infty }^\infty {\int\limits_{ - 1}^1 {\frac{{\Delta g\left( {x,kx} \right) - \Delta {g_0} - {g_x}x - {g_y}kx}}{{{x^2}{{\left( {1 + {k^2}} \right)}^{3/2}}}}} dxdk} + \frac{1}{{2\pi }}\int\limits_{ - \infty }^\infty {\int\limits_{ - 1}^1 {\frac{{\Delta g\left( {ky,y} \right) - \Delta {g_0} - {g_x}ky - {g_y}y}}{{{y^2}{{\left( {1 + {k^2}} \right)}^{3/2}}}}dydk} } $$
(6)

Bian (1997) has computed the integral in Eq. 6 by means of unit quadrature structure. In the present work, as a generalized case, the actual grid interval of length “a” km has been implemented and used for the solution invoking Newton–Cotes formulas up to n = 6. The solution has been attempted for only even orders of integration since these formulas with even numbers of strips give zero error for polynomials, see Williams (1972). The general form of the Newton–Cotes formula can be written as

$$ \int\limits_a^b {f(x){\hbox{d}}x = {A_0}h} \sum\limits_{i = 0}^n {{w_i}{f_i} + {A_1}{h^{k + 1}}{f^k}\left( \chi \right)} $$
(7)

where n is the number of strips, \( h = \left( {b - a} \right)/n \) and χ is some value in the interval (a, b). The 2nd term on the right side of Eq. 7 is the amount of expected truncation error. The grid interval “a” is taken with square grid to solve for all of the following intervals of the innermost area

  1. 1.

    −a < x < a, −a < y < a (n = 2),

  2. 2.

    −2a < x < 2a, −2a < y < 2a (n = 4),

  3. 3.

    −3a < x < 3a, −3a < y < 3a (n = 6),

as shown in the Fig. 1.

Fig. 1
figure 1

The innermost area used for the computation and the moving window (used for n = 6)

Full size image

Solutions for the different areas of integration

The complete derivations for different areas of integration have been carried out for the solution of the vertical gravity gradient integral with constant grid interval “a”. The solution is derived for all three cases of the innermost area, i.e. from −a < x < a, −a < y < a to −3a < x < 3a, −3a < y < 3a to be applicable for different orders of integration with actual grid interval, which is taken constant in both directions. One can extend this integration for different grid dimension in x and y directions very easily. There could be three cases in this study which require independent and separate mathematical treatment.

Solution for the innermost area −a < x < a, −a < y < a

To find the solution of the innermost area of −a < x < a, −a < y < a, Simpson's rule for n = 2 has been used as a starting point. The solution for the Eq. 6 using Simpson's rule twice i.e. in the x-direction and then y-direction for the innermost area with grid interval “a” (Bian 1997)

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{1}{{3a\;\pi }}\left[ {\begin{array}{*{20}{c}} {\int\limits_{ - 1}^1 {\frac{{\Delta g\left( { - a, - ak} \right) + 2\left( {{g_{xx}} + {g_{yy}}{k^2}} \right) + \Delta g\left( {a,ak} \right)}}{{{{(1 + {k^2})}^{3/2}}}}{\hbox{d}}k} + } \hfill \\{\int\limits_{ - 1}^1 {\frac{{\Delta g\left( { - ak, - a} \right) + 2\left( {{g_{xx}}{k^2} + {g_{yy}}} \right) + \Delta g\left( {ak,a} \right)}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}{\hbox{d}}k} } \hfill \\\end{array} } \right] $$
(8)

The middle terms of both the above integrals are obtained from the following identities using limiting values in terms of Eqs. 3 and 4 we have

$$ {\lim_{x \to 0}}\frac{{\Delta g\left( {x,kx} \right) - \Delta {g_0} - {g_x}x - {g_y}kx}}{{{x^2}}} = \frac{{{g_{xx}}}}{2} + {g_{xy}}k + \frac{{{g_{yy}}}}{2}{k^2} < \infty $$
(9)
$$ {\lim_{y \to 0}}\frac{{\Delta g\left( {ky,y} \right) - \Delta {g_0} - {g_x}ky - {g_y}y}}{{{y^2}}} = \frac{{{g_{yy}}}}{2} + {g_{xy}}k + \frac{{{g_{xx}}}}{2}{k^2} < \infty $$
(10)

Where, g xy k goes to zero when integrated for limits −1 ≤ x ≤ 1and −1 ≤ y ≤ 1 being an odd function of “k”.

Now, after performing integration using the Simpson's rule at x = 0, ±1 and y = 0, ±1, and considering relationships for the two horizontal derivatives

$$ {g_{xx}} = \left[ {\Delta g\left( { - a,0} \right) + \Delta g\left( {a,0} \right) - 2\Delta g\left( {0,0} \right)} \right]/{a^2} $$
(11)

and

$$ {g_{yy}} = \left[ {\Delta g\left( {0, - a} \right) + \Delta g\left( {0,a} \right) - 2\Delta g\left( {0,0} \right)} \right]/{a^2} $$
(12)

We arrive at the final solution as

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{1}{{18a\pi \sqrt {2} }}\left[ {\Delta g\left( { - a,a} \right) + \Delta g\left( {a, - a} \right) + \Delta g\left( { - a, - a} \right) + \Delta g\left( {a,a} \right) - 4\Delta g\left( {0,0} \right)} \right] + \left( {\frac{2}{{3a\pi }}\ln \left( {1 + \sqrt {2} } \right) + \frac{2}{{9a\pi }}} \right)\left[ {\Delta g\left( {a,0} \right) + \Delta g\left( { - a,0} \right) + \Delta g\left( {0,a} \right) + \Delta g\left( {0, - a} \right) - 4\Delta g\left( {0,0} \right)} \right] $$
(13)

The solution for the grid interval “a” as derived in Heiskanen and Moritz (1967) for planar earth approximation is given by

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{{{s_0}}}{{4{a^2}}}\left[ {\Delta g\left( {a,0} \right) + \Delta g\left( { - a,0} \right) + \Delta g\left( {0,a} \right) + \Delta g\left( {0, - a} \right) - 4\Delta g\left( {0,0} \right)} \right] $$
(14)

where so is the spatial distance from the computation point in the xy-plane, tangent to the earth surface. On comparing the Eqs. 13 and 14, it can be seen that the main difference in both the equations is the first term in Eq. 13. Here, this difference comes from the effect of points lying at the corners of the square grid for the integration area. This additional contribution can be seen clearly from the comparison of results, derived from Eqs. 13 and 14 and shown in Table 2 in “Solutions for the different areas of integration.”

Solution for the innermost area, −2a < x < 2a, −2a < y < 2a

The Newton–Cotes formula for n = 4 has been used for this solution for the integration area bounded by −2a < x < 2a, −2a < y < 2a. After applying the above integration limits in Eq. 8 and arranging the terms with same coefficients we get.

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{a}{{45\pi }}\int\limits_{ - 1}^1 {\left\{ {\begin{array}{*{20}{c}} {\frac{7}{{4{a^2}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 2a, - 2ak} \right) + \Delta g\left( {2a,2ak} \right) + } \hfill \\{\Delta g\left( { - 2ak, - 2a} \right) + \Delta g\left( {2ak,2a} \right) - 4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right)} \hfill \\+ \hfill \\{\frac{{32}}{{{a^2}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - a, - ak} \right) + \Delta g\left( {a,ak} \right) + \Delta g\left( { - ak, - a} \right) + } \hfill \\{\Delta g\left( {ak,a} \right) - 4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right)} \hfill \\\end{array} } \right\}\frac{{{\text{d}}k}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}} + \frac{{2a}}{{15\pi }}\int\limits_{ - 1}^1 {\left( {{g_{xx}} + {g_{yy}}} \right)\frac{{{\text{d}}k}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}} $$
(15)

Now, we apply the Newton–Cotes integration on the 1st term of 1st integral at x = 0, ±1/2, ±1 and y = 0, ±1/2, ±1 and 2nd term of the 1st integral at x = 0, ±1 and y = 0, ±1. The 2nd integral term on right side of Eq. 15 can be solved analytically for the limits of x = 0, ±1 and y = 0, ±1, so we arrive at the final solution after arranging the similar terms,

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{a}{{45\pi }}\left[ {\begin{array}{*{20}{c}} {\frac{1}{{3{a^2}}}\left( {36\ln \left( {1 + \sqrt {2} } \right) + 128} \right)\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - a,0} \right) + \Delta g\left( {a,0} \right) + \Delta g\left( {0, - a} \right) + \Delta g\left( {0,a} \right)} \hfill \\{ - 4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{\frac{7}{{180{a^2}}}\left\{ {\begin{array}{*{20}{c}} {\frac{7}{{\sqrt {2} }}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 2a,2a} \right) + \Delta g\left( {2a, - 2a} \right) + \Delta g\left( {2a,2a} \right) + \Delta g\left( { - 2a, - 2a} \right)} \hfill \\{ - 4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{12\left( {\Delta g\left( { - 2a,0} \right) + \Delta g\left( {2a,0} \right) + \Delta g\left( {0,2a} \right) + \Delta g\left( {0, - 2a} \right) - 4\Delta g\left( {0,0} \right)} \right) + } \hfill \\{\frac{{256}}{{\sqrt {{125}} }}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 2a,a} \right) + \Delta g\left( {2a, - a} \right) + \Delta g\left( {a, - 2a} \right) + \Delta g\left( { - a,2a} \right) + } \hfill \\{\Delta g\left( { - 2a, - a} \right) + \Delta g\left( {2a,a} \right) + \Delta g\left( { - a, - 2a} \right) + \Delta g\left( {a,2a} \right)} \hfill \\{ - 8\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right)} \hfill \\\end{array} } \right\}} \hfill \\+ \hfill \\{\frac{{32}}{{3{a^2}}}\left( {\Delta g\left( { - a,a} \right) + \Delta g\left( {a, - a} \right) + \Delta g\left( { - a, - a} \right) + \Delta g\left( {a,a} \right) - 8\Delta g\left( {0,0} \right)} \right)/\sqrt {2} } \hfill \\\end{array} } \right] $$
(16)

The solutions of Eqs. 13 and 15 are similar to those computed by Bian (1997), where, only unit grid interval was implemented.

Solution for the inner most area −3a < x < 3a, −3a < y < 3a

The solution procedure has been developed in the similar manner, as above, by using the Newton–Cotes formula (for n = 6) for the above mentioned innermost area. After applying the required integration limits in Eq. 8 and arranging the terms with same coefficients, we get,

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{a}{{280\pi }}\left[ {\int\limits_{ - 1}^1 {\left\{ {\begin{array}{*{20}{c}} {\frac{{41}}{{9{a^2}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 3a, - 3ak} \right) + \Delta g\left( {3a,3ak} \right) + } \hfill \\{\Delta g\left( { - 3ak, - 3a} \right) + \Delta g\left( {3ak,3a} \right)} \hfill \\{ - \left. {4\Delta g\left( {0,0} \right)} \right)} \hfill \\\end{array} } \right) + } \hfill \\{\frac{{54}}{{{a^2}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 2a, - 2ak} \right) + \Delta g\left( {2a,2ak} \right) + } \hfill \\{\Delta g\left( { - 2ak, - 2a} \right) + \Delta g\left( {2ak,2a} \right) - } \hfill \\{4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{\frac{{27}}{{{a^2}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - a, - ak} \right) + \Delta g\left( {a,ak} \right) + } \hfill \\{\Delta g\left( { - ak, - a} \right) + \Delta g\left( {ak,a} \right) - } \hfill \\{4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right)} \hfill \\\end{array} } \right\}\frac{{{\text{d}}k}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}} } \right] + \frac{{93a}}{{140\pi }}\int\limits_{ - 1}^1 {\left( {{g_{xx}} + {g_{yy}}} \right)\frac{{{\text{d}}k}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}} $$
(17)

Again, we can solve the 1st integral term in Eq. 17 by applying the Newton–Cotes formula on the 1st term at x = 0, ±1/3, ±2/3, ±1 and y = 0, ±1/3, ±2/3, ±1, the 2nd term at x = 0, ±1/2, ±1 and y = 0, ±1/2, ±1, respectively. The 2nd integral term on the right side of Eq. 17 has been solved analytically. For the sake of simplicity at x = 0, ±1 and y = 0, ±1. The whole solution can be split into different terms as shown in the following relationship.

$$ \frac{{\partial \Delta g}}{{\partial H}} = \frac{1}{{2\pi }}\frac{a}{{140}}\left( {T1 + T2 + T3 + T4} \right) $$
(18)

where

$$ T1 = \frac{{41}}{{3,780{a^2}}}\left\{ {\begin{array}{*{20}{c}} {\frac{{41}}{{\sqrt {2} }}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 3a,3a} \right) + \Delta g\left( {3a, - 3a} \right) + } \hfill \\{\Delta g\left( { - 3a, - 3a} \right) + \Delta g\left( {3a,3a} \right) - } \hfill \\{4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{\frac{{216}}{{{{(13/9)}^{3/2}}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 3a,2a} \right) + \Delta g\left( {3a, - 2a} \right) + } \hfill \\{\Delta g\left( {2a, - 3a} \right) + \Delta g\left( { - 2a,3a} \right) + } \hfill \\{\Delta g\left( { - 3a, - 2a} \right) + \Delta g\left( {3a,2a} \right) + } \hfill \\{\Delta g\left( { - 2a, - 3a} \right) + \Delta g\left( {2a,3a} \right) - } \hfill \\{8\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{\frac{{27}}{{{{(10/9)}^{3/2}}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 3a,a} \right) + \Delta g\left( {3a, - a} \right) + } \hfill \\{\Delta g\left( {a, - 3a} \right) + \Delta g\left( { - a,3a} \right) + } \hfill \\{\Delta g\left( { - 3a, - a} \right) + \Delta g\left( {3a,a} \right) + } \hfill \\{\Delta g\left( { - a, - 3a} \right) + \Delta g\left( {a,3a} \right) - } \hfill \\{8\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{272\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 3a,0} \right) + \Delta g\left( {3a,0} \right) + \Delta g\left( {0, - 3a} \right)} \hfill \\{ + \Delta g\left( {0,3a} \right) - 4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right)} \hfill \\\end{array} } \right\} $$
(19)
$$ T2 = \frac{6}{{5{a^2}}}\left\{ {\begin{array}{*{20}{c}} {\frac{7}{{\sqrt {2} }}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 2a,2a} \right) + \Delta g\left( {2a, - 2a} \right) + } \hfill \\{\Delta g\left( { - 2a, - 2a} \right) + \Delta g\left( {2a,2a} \right) - } \hfill \\{8\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{\frac{{32}}{{\sqrt[3]{{1.25}}}}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 2a,a} \right) + \Delta g\left( {2a, - a} \right) + } \hfill \\{\Delta g\left( {a, - 2a} \right) + \Delta g\left( { - a,2a} \right) + } \hfill \\{\Delta g\left( { - 2a, - a} \right) + \Delta g\left( {2a,a} \right) + } \hfill \\{\Delta g\left( { - a, - 2a} \right) + \Delta g\left( {a,2a} \right) - } \hfill \\{8\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{12\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - 2a,0} \right) + \Delta g\left( {2a,0} \right) + } \hfill \\{\Delta g\left( {0, - 2a} \right) + \Delta g\left( {0,2a} \right) - } \hfill \\{4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right)} \hfill \\\end{array} } \right\} $$
(20)
$$ T3 = \frac{9}{{{a^2}}}\left\{ {\begin{array}{*{20}{c}} {\frac{1}{{\sqrt {2} }}\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - a,a} \right) + \Delta g\left( {a, - a} \right) + } \hfill \\{\Delta g\left( { - a, - a} \right) + \Delta g\left( {a,a} \right) - } \hfill \\{4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right) + } \hfill \\{4\left( {\begin{array}{*{20}{c}} {\Delta g\left( { - a,0} \right) + \Delta g\left( {a,0} \right) + } \hfill \\{\Delta g\left( {0, - a} \right) + \Delta g\left( {0,a} \right) - } \hfill \\{4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right)} \hfill \\\end{array} } \right\} $$
(21)

and

$$ T4 = \frac{{372\ln \left( {1 + \sqrt {2} } \right)}}{{{a^2}}}\left\{ {\begin{array}{*{20}{c}} {\Delta g\left( { - a,0} \right) + \Delta g\left( {a,0} \right) + } \hfill \\{\Delta g\left( {0, - a} \right) + \Delta g\left( {0,a} \right) - } \hfill \\{4\Delta g\left( {0,0} \right)} \hfill \\\end{array} } \right\} $$
(22)

Data processing methodology and analysis of the results

The observed gravity and orthometric height data was available with random distribution for this study. The elevation data used here was a part of the gravity surveys and obtained through levelling measurements. Therefore, an assumption has been made here that the use of orthometric height instead of normal height would have insignificant effect on the computed geoid-quasigeoid separation in the low to mid-range elevations. So, orthometric height has been used in the determination of theoretical normal gravity and other related gravity field parameters e.g. gravity anomalies and terrain correction etc. The terrain correction was computed with the help of SRTM30 DEM using GRAVSOFT (2005) package. The DEM was mainly used for fine and reference grids development as per requirements of. The statistics of the related parameters are given in Table 1. The gridding of the free air anomaly and elevation data was made with the Kriging routine of Golden Sotware Surfer.7. The grid interval for all parameters comprising elevation, free air anomaly, Bouguer anomaly and average theoretical normal gravity was consistently used as 0.030303° for the computation of the gravity anomaly gradient and onward corresponding geoid-quasigeoid separation term dependent on it. Equivalently, an average grid interval in terms of the spatial distance of 3.1 km (∼3.096 km corresponding to the integration radius for n = 2 in the study area) was implemented for computation of the vertical gravity anomaly gradient in planar approximation. The higher orders of integration, i.e., n = 4 and n = 6, correspond to the maximum integration radii of 6.2 and 9.3 km, respectively.

Table 1 Statistics of different parameters of gravity and elevation
Full size table

The solutions of Eqs. 13, 14, 16 and 18 have been implemented in the area bounded by longitudes 70°E to 73°E and latitudes 30°N to 33°N using routines developed in MATLAB (Appendix 1). This provides the basis to estimate the size and variability of the vertical gravity anomaly gradient and to compute the corresponding effect on the geoid-quasigeoid separation for different integration radii. Different combination solutions of the innermost zone with increasing integration radii were studied in planar approximation to assess the generally feasible integration radius for the best estimate of vertical gravity anomaly gradient.

The values of the vertical gravity anomaly gradient for different orders of integration are shown in Table 2. The vertical gravity anomaly gradient for the 2nd order planar approximations varies from −8.954 to −9.874 mGal/km in the study area. Here a significant difference (∼14%) is observed in the mean value computed by the planar approximation according to Eq. 13 given in Heiskanen and Moritz (1967) and the one derived from Eq. 14 using the Newton–Cotes integration with Simpson's rule for n = 2. The difference between two results emerges due to the contributions of points at the corners of the grid along the diagonals which are included in Eq. 13. However, it has negligible effect on the corresponding geoid-quasigeoid separation. These two formulas cover only the effects of the adjacent grid points and need to be extended to higher orders to include the effects of broader surrounding area for better estimates of the vertical gravity anomaly gradient. Therefore, same data was processed for the 4th and the 6th order Newton–Cotes integration formulas. The effect of an extended integration radius is evident from the results shown in Table 2, where the mean value has become double for the 4th/6th order as compared to the corresponding 2nd order computed value.

Table 2 Statistics of vertical gravity anomaly gradient (mGal/km) for different orders of integration
Full size table

An increasing trend can be seen in the values of the vertical gravity anomaly gradient with increasing integration radii, which reaches almost a saturation value for the orders of n = 4 and n = 6. These results also indicate that mean difference from the 4th to 6th order is much lower than that from 2nd to 4th order of integration. The mean deviation of the vertical gravity gradient for the 4th and the 6th order is only 1.3% whereas the difference of mean value of 4th order from the 2nd order, i.e. Eqs. 14 and 16, is significantly higher and considerable i.e. ∼47%. Therefore, it appears feasible to extend the integration to the 4th or 6th orders only in the innermost zone due to the obvious fact that vertical gravity anomaly gradient term has relatively small effect on the geoid-quasigeoid separation. The results shown in Table 3 are obvious proof of this fact.

Table 3 Statistics of the geoid-quasigeoid separation term (mm) dependent on the vertical gravity anomaly gradient part (2nd term of Eq. 1) determined by using different integration radii
Full size table

There is a similar pattern of the geoid-quasigeoid separation term dependent on vertical gravity anomaly gradient with increasing order of the integration due to the fact that the geoid-quasigeoid separation is linearly dependent on the vertical gravity anomaly gradient. It has a maximum absolute value on the order of 3 mm for the 2nd order and reaches 6.18 mm with absolute maximum for the 6th order of the integration. Its value can be considered almost negligible i.e. only ∼3% as compared to the Bouguer anomaly dependent absolute value of 0.206 m for this elevation range. These results indicate that the order of integration may be selected either n = 4 or n = 6 for better estimates of vertical gravity anomaly gradient and the onward geoid-quasigeoid separation term part dependent on it.

The estimated geoid-quasigeoid term dependent on the vertical gravity anomaly gradient based on the planar approximation shows that it has very minor effects and are on the order of a centimeter in mid-range elevation. This term has a linear relationship with the squared elevation. However, the effect of elevation is not so pronounced and increasing with elevation. This is due to the reason that the parameters like the trend of the free air anomaly and the derived vertical gravity anomaly gradient also have direct and linear relationship with geoid-quasigeoid separation. This fact is evident from Fig. 2 which indicates that most of the points have very small effects even for highest elevations of ∼1,500 m. The diversity of terrain elevation in Pakistan (from 0 to ∼8,000 m) can suggest that it could reach to sub decimetre level. These results have been confirmed in another study where this term has been found to range between −20 to 30 mm all over Pakistan (Sadiq et al. 2009).The major part of the geoid-quasigeoid separation term derives from the Bouguer anomaly term (∼97%) and it may well approximate it in low and mid range elevations. Therefore, in case of very high terrain elevations it is mandatory to include the vertical gravity anomaly gradient dependent term due to the reason that it also depends on the 2nd order height term, i.e. H2. The variance of gravity anomaly and its gradient with elevation can also be considered as deciding parameter for its inclusion in the total part.

Fig. 2
figure 2

Vertical gravity anomaly gradient dependent part of the geoid-quasigeoid separation term (m) as a function of elevation with integration radius of 9.3 km for n = 6

Full size image

A comparative study for the change of both the vertical gravity anomaly gradient and its dependent part of geoid-quasigeoid separation with integration radii has also been conducted to predict the behaviour of these closely related parameters. For this purpose standard deviation was selected as a parameter due to the reason that it better estimates the range over which data are distributed around the mean value. The scaled standard deviations relative to the values of vertical gravity anomaly gradient determined using Eq. 14 and the corresponding term of the geoid-quasigeoid separation, i.e. \( \frac{{{H^2}\left( \Omega \right)}}{{2\overline \gamma }}\left( {\frac{{\partial \Delta {g^F}}}{{\partial H}}} \right) \), were plotted against different integration radii (Fig. 3). The rate of change of the vertical gravity anomaly gradient is comparatively higher than the corresponding geoid-qausigeoid separation and this is evident from the relatively faster saturation obtained for respective integration radii. This shows that vertical gravity anomaly gradient has relatively strong dependence on the integration radius as compared with its geoid-quasigeoid separation counterpart. On the other hand, it should be kept in mind that the integration radius also depends on the characteristics of terrain and characteristics of free air anomaly in a particular area. This characteristic of geoid-quasigeoid separation may also be associated to its dependence on the 2nd order height values and average theoretical normal gravity.

Fig. 3
figure 3

Plot of scaled standard deviation of the vertical gravity anomaly gradient and the geoid-quasigeoid separation with integration radii relative to the linear approximation by Heiskanen and Moritz (1967)

Full size image

In the sequel, one can conclude that both the vertical gravity anomaly gradient and its corresponding effect on the geoid-quasigeoid separation are significant even for the innermost area. The method proposed by Heiskanen and Moritz (1967) for vertical gravity anomaly gradient (Eq. 14) is too simple and approximate and need extended radius of integration (Bian 1997). It can be seen explicitly from the treatise developed above that effects of vertical gravity anomaly gradient and its corresponding geoid-quasigeoid separation term can be approximated in an innermost area with reasonable accuracy by using formulae with integration radii corresponding to the orders of n = 4 and n = 6. This fact may be attributed mainly to the fast decreasing characteristics of the integral for the solution of vertical gravity anomaly gradient with distance and partly to the accurate solution of Simpson’s rule for even orders of integration.

Conclusions and recommendations

The computation and comparative study of different relationships of vertical gravity anomaly gradient for different integration radii was carried out to asses and select an optimum radius of integration with planar approximation in the innermost area. The results lead to the following main conclusions.

  1. 1.

    The integration radii of 3.1 km (n = 2), 6.2 km (n = 4) and 9.3 km(n = 6) exhibit trend of saturation as per expectation due to inclusion of more surrounding effects. The results show that vertical gravity anomaly gradient reaches to a saturation value for n = 4 as its mean value differs only by 1.3% from next higher order integration radius mean value i.e. n = 6.

  2. 2.

    The characteristic of the geoid-quasigeoid separation term dependent on vertical gravity anomaly gradient has similar to the trend of the vertical gravity anomaly gradient as a function of the integration order.

  3. 3.

    The order of integration could be selected either n = 4 or n = 6 for better estimates of vertical gravity anomaly gradient and onward geoid-quasigeoid separation term dependent on it in the areas of low- to mid-range elevation.

  4. 4.

    The vertical gravity anomaly gradient dependent geoid-quasigeoid separation term has insignificant value in low- and mid-range elevation ranges and may be neglected for most practical purposes. It may have of the order of sub decimetre in the very high mountains of northern Pakistan.

  5. 5.

    The changes in the vertical gravity anomaly gradient are not significant enough (Table. 2) for n = 4 to n = 6, due to fast decreasing character of integral, hence the effect of outer zone can be neglected due to its small effect on geoid-quasigeoid separation. However, outer zone contribution may be evaluated in the sense to quantify its real effect on the geoid-quasigeoid separation term.

In the end, it is pertinent to compare these results with other methods such as those recently developed by Tenzer et al. (2006) in addition to the spherical approximation instead of planar and use of more rugged terrain elevation.