1 Introduction

Mars is the fourth planet in the solar system. It is a small rocky planet, with a thin atmosphere, and cold surface. Despite these attributes, Mars is the most Earth-like planet in the solar system, and has following Earth-like features: It has mountains, valleys, near 24-h rotation period, tilted rotational axis and seasons. Due to those features, and its closeness to the Earth, Mars has been the subject of the most of planetary investigations and studies about the solar system plants other than the Earth. Geodesists are also interest in applying their expertise, i.e. geometry (shape) and gravity, to Mars. Development of shape models, gravity field models, reference equipotential surface of Mars (Areoid), and reference ellipsoid, are examples of those involvements.

Among the geodetic achievements from the Mars missions is the shape model of Mars in terms of spherical harmonics. GTM090AA developed by Smith et al. (1999a) is an example of shape models of Mars. This model, in terms of spherical harmonics to degree and order 90, is computed based on the data from the Mars Orbiter Laser Altimeter (MOLA) on the Mars Global Surveyor (MGS) spacecraft. Coordinate system of GTM090AA model is referred to IAU 1991 (Davies et al. 1992), which is an areocentric (centered at the center of mass of Mars) coordinate system, with longitudes measured positive east. This model gives the shape of Mars as absolute radius from the center of mass of Mars. The model is a numerical transform of planetary radius values from the MOLA 0.5° gridded data set into spherical harmonics using Simpson’s Rule quadrature. The gridded data set includes all MOLA nadir observations from the Orbit Insertion phase, plus Mapping Phase nadir observations from the beginning of mapping phase through the end of December 1999. Additionally, off-nadir observations of the North Pole were included from 87 N latitude northward, taken during the spring of 1998, and of both poles taken in July 1999, from 87 N and S to the poles. Figure 1 shows the topography of Mars, computed by Smith et al. (1999a) based on MGS observations.

Fig. 1
figure 1

The Global topography of Mars by Smith et al. (1999a)

Full size image

Other products of geodesists from studies related to Mars, are gravity field models which at beginning was developed based on two natural satellites of Mars, i.e. Phobos and Deimos, and later by radio tracking of successive Mars orbiters missions started by Mariner 9 in 1971. There have been two independent efforts to improve the gravitational field of Mars in terms of spherical harmonics, by Goddard Space Flight Center (GSFC) and Jet Propulsion Laboratory (JPL) with results initially reported by Smith et al. (1999b). Subsequent analysis was reported by Lemoine et al. (2001) for GSFC model GMM-2B (to degree and order 80) and by Yuan et al. (2001) for JPL model MGS75D (to degree and order 75). The newer GSFC gravitational models are GGM1025 (to degree and order 80) and GGM1041C (to degree and order 90) and the newer JPL gravitational models are MGS75E, MGS85F, MGS85F2, MGS85H2, MGS95I, and MGS95 J (where the degree and order of expansion is included in the name of models). The gravitational models from JPL and GSFC have included roughly the same MGS data sets and provide results that are mostly consistent. More details about the GSFC and JPL models can be found in Konopliv et al. (2006). Besides, the GSFC and JPL models, MGGM08A Martian gravitational field model complete to degree and order 95 developed by a European team can be mentioned (see Marty et al. 2009 for details).

Having the gravitational field and shape models of Mars available, the reference equipotential surface (i.e. Areoid) can be computed. Correlations between topography, and gravity (or Areoid), can provide important constraints on the structure and tectonic evolution of Mars, e.g. to estimate the Martian crustal and lithospheric thicknesses and elasticity. Besides, fluvial geomorphology studies and near Mars surface maneuvers like entry trajectory terminal descent, parachute deploys, and rover operations analysis, requires representation of Mars topography with respect to Areoid, i.e. planetopotential topography. The classical method for Areoid determination is as follows: First the mean radius of the reference equipotential (Areoid) at equator is determined (in case of Mars the radius 3,396,000 m is adopted). Then the radius of Areoid at a given latitude and longitude is found iteratively by evaluating the potential model at that specific location until the matching to the mean equatorial potential is obtained. Figure 2, shows the Areoid computed by Smith et al. (1999b) based on MGS75B gravitational field model.

Fig. 2
figure 2

The Areoid computed by Smith et al. (1999b) based on MGS75B gravitational field model. The heights are in meter

Full size image

In this paper, in contrast to classical methods, we compute the Areoid as an equipotential surface of the gravity field of Mars, which best fits to its shape in least square sense via solution of a constrained optimization problem. The solution of this constrained optimization problem results in both the geometry (shape) and gravity potential value of Areoid. We have borrowed this idea from Gauss-Listing (Gauss 1828; Listing 1873) definition for the geoid as an equipotential surface of the Earth gravity field which best fits to global Mean Sea Level (MSL) in least squares sense. In this way the computed Areoid represents the surface of a liquid Mars in hydrostatic equilibrium.

Another subject of interest of geodesists is approximation of the shape of planets by mathematical surfaces. Those mathematical surfaces in different approximations are usually sphere, bi-axial and tri-axial ellipsoids. Duxbury et al. (2002) describes the computation of a reference spheroid of the kind sphere, bi-axial ellipsoid and tri-axial ellipsoid, centered and not centered at the centre of mass of Mars via best fitting of those spheroids to the shape of Mars. Duxbury et al. (2002) recommends the use of bi-axial ellipsoid centered at the centre of mass of Mars as the reference ellipsoid and gives the values (3,396,190 ± 100) (m) for equatorial and (3,376,200 ± 100) (m) for polar radii of such a best fitting reference ellipsoid to the shape of Mars. In this paper in contrast to previous efforts, we will apply the concept of Somigliana-Pizzetti type reference ellipsoid to Mars. Computations of Somigliana-Pizzetti type reference ellipsoid according to Grafarend and Ardalan (1999) is based on four fundamental gravity parameters. Those parameters in case of Mars are {W 0, GM, ω, J 20}, i.e. {Areoid’s gravity potential, gravitational constant of Mars, angular velocity of Mars, second zonal spherical harmonic of gravitational field expansion of Mars}. The Somigliana-Pizzetti type reference ellipsoid unlike the so-far computed reference ellipsoids, which are merely geometric surfaces, has also physical properties, such as being an equipotential surface of the normal gravity field generated by a massive rotational ellipsoid having the mass equal to that of Mars, rotating with the angular velocity of Mars, with rotational symmetry mass distribution (Heiskanen and Mortiz 1967; Vanicek and Krakiwsky 1986; Ardalan and Grafarend 2001). Somigliana-Pizzetti type reference ellipsoid could be the hydrostatic equilibrium figure of a liquid Mars with rotational symmetry mass distribution.

In Sect. 2, the mathematical details of Areoid computations via our constrained optimization problem will be discussed. Section 3 is devoted to the error budget of Areoid computations based on our method. The numerical results of Areoid computations will be given in Sect. 4. In Sect. 5, the computations of the reference ellipsoid of Somigliana-Pizzetti type for Mars will be presented, and as the final product the computed Areoid is presented with respect to the computed reference ellipsoid. Finally, conclusions and remarks are given in Sect. 6.

2 Setup of the Constrained Optimization Problem for Areoid Computations

Having defined the Areoid as an equipotential surface of Mars gravity field, which fits to the shape of Mars in least squares sense, we setup a constrained optimization problem by minimizing the spherical radial distances d i between the shape of Mars and Areoid (Fig. 3), under the constraint that all the points on the Areoid have one gravity potential value. This can be written in mathematical notations as follows:

$$ \begin{gathered} {\text{Constraints: }}\,\left\{ \begin{gathered} f_{1} (d_{1} ,d_{2} ) = W_{{P^{\prime}_{1} }} \left( {\phi_{1} ,\lambda_{1} ,r_{1} - d_{1} } \right) - W_{{P^{\prime}_{2} }} \left( {\phi_{2} ,\lambda_{2} ,r_{2} - d_{2} } \right) = 0 \\ f_{2} (d_{1} ,d_{3} ) = W_{{P^{\prime}_{1} }} \left( {\phi_{1} ,\lambda_{1} ,r_{1} - d_{1} } \right) - W_{{P^{\prime}_{3} }} \left( {\phi_{3} ,\lambda_{3} ,r_{3} - d_{3} } \right) = 0 \\ \vdots \\ f_{k - 1} (d_{1} ,d_{n} ) = W_{{P^{\prime}_{1} }} \left( {\phi_{1} ,\lambda_{1} ,r_{1} - d_{1} } \right) - W_{{P^{\prime}_{k} }} \left( {\phi_{k} ,\lambda_{k} ,r_{k} - d_{k} } \right) = 0 \\ \end{gathered} \right. \hfill \\ {\text{Objective function:}}\, \, \sum\limits_{i = 1}^{k} {d_{i}^{2} } \to \min \hfill \\ \end{gathered} $$
(1)

where W, the gravity potential, in terms of spherical harmonics, can be computed as follows:

$$ \begin{gathered} W_{{P^{\prime}_{i} }} \left( {\phi_{i} ,\lambda_{i} ,r_{i} - d_{i} } \right) = {\frac{\text{GM}}{{r_{i} - d_{i} }}}\sum\limits_{n = 0}^{L} {\sum\limits_{m = 0}^{n} {{\frac{{R^{n} }}{{\left( {r_{i} - d_{i} } \right)^{n} }}}} } \bar{P}_{nm} \sin \left( {\phi_{i} } \right)\left( {\bar{C}_{nm} \cos \left( {m\lambda_{i} } \right) + \bar{S}_{nm} \sin \left( {m\lambda_{i} } \right)} \right) \\ + {\frac{1}{2}}\left( {r_{i} - d_{i} } \right)^{2} \omega^{2} \cos^{2} \left( {\phi_{i} } \right) \\ \end{gathered} $$
(2)

where in Eq. 2, R is the mean radius, GM is gravitational constant, ω is angular velocity of Mars, \( \bar{P}_{nm} \sin \left( {\phi_{i} } \right) \) is normalized Legendre function of the first kind, \( \bar{C}_{nm} \) and \( \bar{S}_{nm} \) are fully normalized spherical harmonics, and \( {\tfrac{1}{2}}\left( {r_{i} - d_{i} } \right)^{2} \omega^{2} \cos^{2} \left( {\phi_{i} } \right) \) is centrifugal potential. Eq. 1 can also be written in compact vectorial form as follows:

$$ \left\{ \begin{gathered} {\mathbf{f}}({\mathbf{d}}) = {\mathbf{0}} \hfill \\ {\mathbf{d}}^{T} {\mathbf{d}} \to \min \hfill \\ \end{gathered} \right. $$
(3)

where

$$ {\mathbf{d}} = \left[ {\begin{array}{*{20}c} {d_{1} } & {d_{2} } & \cdots & {d_{k} } \\ \end{array} } \right]^{T} ,\quad {\mathbf{f}} = \left[ {\begin{array}{*{20}c} {f_{1} } & {f_{2} } & \cdots & {f_{k - 1} } \\ \end{array} } \right]^{T} . $$
(4)
Fig. 3
figure 3

Shape of Mars, the Areoid and the spherical radial distances d i that must be minimized in least squares sense under the constraint that the gravity potential value on the Areoid remains constant

Full size image

The constrained optimization problem defined by Eqs. 1 or 3 must be solved for the unknown parameters d i . Here we will follow the Lagrange method after linearizing the function f in Eq. 3 as follows:

$$ {\mathbf{f}}\left( {\mathbf{d}} \right) = 0 \Rightarrow \left. {{\frac{{\partial {\mathbf{f}}}}{{\partial {\mathbf{d}}}}}} \right|_{{{\mathbf{d}} = {\mathbf{d}}_{0} }} \left( {{\mathbf{d}} - {\mathbf{d}}_{0} } \right) + {\mathbf{f}}\left( {{\mathbf{d}}_{0} } \right) \cong 0 \Rightarrow \underbrace {{\left. {{\frac{{\partial {\mathbf{f}}}}{{\partial {\mathbf{d}}}}}} \right|_{{{\mathbf{d}} = {\mathbf{d}}_{0} }} }}_{{{\mathbf{A}}_{k - 1 \times k} }}{\mathbf{d}} = \underbrace {{\left. {{\frac{{\partial {\mathbf{f}}}}{{\partial {\mathbf{d}}}}}} \right|_{{{\mathbf{d}} = {\mathbf{d}}_{0} }} {\mathbf{d}}_{0} - {\mathbf{f}}\left( {{\mathbf{d}}_{0} } \right)}}_{{{\mathbf{w}}_{k - 1 \times 1} }} \Rightarrow {\mathbf{Ad}} = {\mathbf{w}}. $$
(5)

The matrix A shown in Eq. 5 reads as:

$$ {\mathbf{A}}_{k - 1 \times k} = \left[ {\begin{array}{*{20}c} {{\frac{{\partial W_{{P^{\prime}_{1} }} }}{{\partial d_{1} }}}} & { - {\frac{{\partial W_{{P^{\prime}_{2} }} }}{{\partial d_{2} }}}} & 0 & 0 & \cdots & 0 & 0 \\ {{\frac{{\partial W_{{P^{\prime}_{1} }} }}{{\partial d_{1} }}}} & 0 & { - {\frac{{\partial W_{{P^{\prime}_{3} }} }}{{\partial d_{3} }}}} & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {{\frac{{\partial W_{{P^{\prime}_{1} }} }}{{\partial d_{1} }}}} & 0 & 0 & 0 & \cdots & 0 & { - {\frac{{\partial W_{{P^{\prime}_{k} }} }}{{\partial d_{k} }}}} \\ \end{array} } \right]. $$
(6)

Now the linearized optimization problem can be written as follows:

$$ \left\{ \begin{gathered} {\mathbf{Ad}} - {\mathbf{w}} = 0 \hfill \\ {\mathbf{d}}^{T} {\mathbf{d}} \to \min \hfill \\ \end{gathered} \right.. $$
(7)

Introducing the Lagrange function L as:

$$ L\left( {{\mathbf{d}},{\varvec{\uplambda}}} \right) = {\mathbf{d}}^{T} {\mathbf{d}} - 2{\varvec{\uplambda}}^{T} \left( {{\mathbf{Ad}} - {\mathbf{w}}} \right). $$
(8)

Necessary condition to minimize d T d under the constraint Ad – w = 0 according to Lagrange method is as follows:

$$ \left\{ \begin{gathered} {\frac{\partial L}{{\partial {\mathbf{d}}}}} = 0 \hfill \\ {\frac{\partial L}{{\partial {\varvec{\uplambda}}}}} = 0 \hfill \\ \end{gathered} \right.. $$
(9)

Therefore,

$$ \left\{ \begin{gathered} 2{\mathbf{d}}^{T} - 2{\varvec{\uplambda}}^{T} {\mathbf{A}} = 0 \hfill \\ \left( {{\mathbf{Ad}} - {\mathbf{w}}} \right)^{T} = 0 \hfill \\ \end{gathered} \right.. $$
(10)

Or

$$ \left\{ \begin{gathered} {\mathbf{d}} - {\mathbf{A}}^{T} {\varvec{\uplambda}} = 0 \hfill \\ {\mathbf{Ad}} - {\mathbf{w}} = 0 \hfill \\ \end{gathered} \right.. $$
(11)

From the first equation of Eq. 11 we have:

$$ {\mathbf{d}} = {\mathbf{A}}^{T} {\varvec{\uplambda}}. $$
(12)

Upon substitution of Eq. 12 in second equation of Eq. 11 we have:

$$ {\mathbf{AA}}^{T} {\varvec{\uplambda}} - {\mathbf{w}} = 0 $$
(13)

Since A is row full rank, therefore AA T is non-singular and invertible. Therefore, we can solve for λ in Eq. 13 to have:

$$ {\varvec{\uplambda}} = \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} {\mathbf{w}}. $$
(14)

By substituting Eq. 14 in Eq. 12, d can be obtained as follows

$$ {\mathbf{d}} = {\mathbf{A}}^{T} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} {\mathbf{w}}. $$
(15)

Next, by applying the error propagation law to Eq. 15 the covariance matrix of d, i.e. C d , can be computed, given the covariance matrix of w, i.e. C w as shown by Eq. 16.

$$ {\mathbf{C}}_{{\mathbf{d}}} = {\mathbf{A}}^{T} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} {\mathbf{C}}_{{\mathbf{w}}} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} {\mathbf{A}}. $$
(16)

The final solution of the problem can be obtained by iteration until the convergence of the problem, i.e. \( \left\| {{\mathbf{d}}_{i} - {\mathbf{d}}_{i - 1} } \right\| < \varepsilon \), is achieved. As the initial values for the unknowns (d 0) zero can also be selected.

3 Error Analysis for Areoid Computations

Now we are going to show how from Eq. 16 the accuracy of the computed Areoid can be derived. For this purpose we have to start from the computation of C w . According to Eqs. 1 and 5, w is a function of potential difference between two points. If we consider ith element of w, we have:

$$ w_{i} = f\left( {\Updelta W_{{1\left( {i + 1} \right)}} } \right) $$
(17)

where ΔW 1(i + 1) is potential difference between the 1st and the (i + 1)th points. Therefore, in general, we need the precision of potential difference between two points. This potential difference is computed based on spherical harmonic expansions, which for i and j points can be written as follows:

$$ \begin{aligned} \Updelta W_{ij} = & {\frac{\text{GM}}{{r_{i} }}}\sum\limits_{n = 0}^{L} {\sum\limits_{m = 0}^{n} {{\frac{{R^{n} }}{{r_{i}^{n} }}}} } \bar{P}_{nm} \sin \left( {\varphi_{i} } \right)\left( {\bar{C}_{nm} \cos \left( {m\lambda_{i} } \right) + \bar{S}_{nm} \sin \left( {m\lambda_{i} } \right)} \right) + {\frac{1}{2}}r_{i}^{2} \omega^{2} \cos^{2}\left( {\varphi_{i} } \right) \\ - {\frac{\text{GM}}{{r_{j} }}}\sum\limits_{n = 0}^{L} {\sum\limits_{m = 0}^{n} {{\frac{{R^{n} }}{{r_{j}^{n} }}}} } \bar{P}_{nm} \sin \left( {\varphi_{j} } \right)\left( {\bar{C}_{nm} \cos \left( {m\lambda_{j} } \right) + \bar{S}_{nm} \sin \left( {m\lambda_{j} } \right)} \right) + {\frac{1}{2}}r_{j}^{2} \omega^{2} \cos^{2} \left( {\varphi_{j} } \right). \\ \end{aligned} $$
(18)

Since in the spherical harmonic expansion of Eq. 18, the zero degree and order term is the most dominant one, we may simplify it by considering only the zero degree and order term as follows.

$$ \Updelta W_{ij} = {\frac{\text{GM}}{{r_{i} }}} + {\frac{1}{2}}r_{i}^{2} \omega^{2} \cos^{2} \left( {\varphi_{i} } \right) - {\frac{\text{GM}}{{r_{j} }}} - {\frac{1}{2}}r_{j}^{2} \omega^{2} \cos^{2} \left( {\varphi_{j} } \right). $$
(19)

Owing to the fact that among the parameters appearing in Eq. 19 the radial distances r i and r j are less accurate one as compared to other parameters φ i , φ j , GM, and ω and the fact that the effect of r i and r j on \( \sigma_{{\Updelta W_{ij} }}^{2} \) is more than the other parameters, in application of error propagation law to Eq. 19 we only consider the errors of r i and r j to have:

$$ \begin{aligned} \sigma_{{\Updelta W_{ij} }}^{2} = & \left( {{\frac{{\partial \Updelta W_{ij} }}{{\partial r_{i} }}}} \right)^{2} \sigma_{{r_{i} }}^{2} + \left( {{\frac{{\partial \Updelta W_{ij} }}{{\partial r_{j} }}}} \right)^{2} \sigma_{{r_{j} }}^{2} \\ = \left( { - {\frac{\text{GM}}{{r_{i}^{2} }}} + r_{i} \omega^{2} \cos^{2} \left( {\varphi_{i} } \right)} \right)^{2} \sigma_{{r_{i} }}^{2} + \left( { - {\frac{\text{GM}}{{r_{j}^{2} }}} + r_{j} \omega^{2} \cos^{2} \left( {\varphi_{j} } \right)} \right)^{2} \sigma_{{r_{j} }}^{2} \\ \end{aligned} $$
(20)

Or by assuming r i  = r j  = R and \( \sigma_{{r_{i} }}^{2} = \sigma_{{r_{j} }}^{2} = \sigma_{r}^{2} \), we have:

$$ \sigma_{{\Updelta W_{ij} }}^{2} = \left( {\left( { - {\frac{\text{GM}}{{R^{2} }}} + R\omega^{2} \cos^{2} \left( {\varphi_{i} } \right)} \right)^{2} + \left( { - {\frac{\text{GM}}{{R^{2} }}} + R\omega^{2} \cos^{2} \left( {\varphi_{j} } \right)} \right)^{2} } \right)\sigma_{r}^{2} . $$
(21)

According to Eq. 21 the maximum value of \( \sigma_{{\Updelta W_{ij} }}^{2} \) occurs for φ i  = φ j  = 90, and in that case we have:

$$ \sigma_{{\Updelta W_{ij} }}^{2} = 2\left( { - {\frac{\text{GM}}{{R^{2} }}}} \right)^{2} \sigma_{r}^{2} . $$
(22)

Or

$$ \sigma_{{\Updelta W_{ij} }} = \pm \sqrt 2\;{\frac{\text{GM}}{{R^{2} }}}\sigma_{r} . $$
(23)

From Eq. 23 by considering σ r  = ±13 m (Smith et al. 1999a) standard deviation of ∆W ij can be derived as \( \sigma_{{\Updelta W_{ij} }} = 68.23\,\,{\text{m}}^{ 2} / {\text{s}}^{ 2} \). Now if we assume that the potential differences are independent, then

$$ {\mathbf{C}}_{{\mathbf{w}}} = 68.23{\mathbf{I}}. $$
(24)

Having computed C w , the covariance matrix of Areoid can be computed from Eq. 16. Finally as a scalar accuracy measure of the computed Areoid we introduce a mean precision as follows:

$$ \begin{aligned} \bar{\sigma }_{\text{Areoid}}^{2} = & {\frac{{{\text{trace}}\left( {{\mathbf{C}}_{{\mathbf{d}}} } \right)}}{n}} = {\frac{{{\text{trace}}\left( {{\mathbf{A}}^{T} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} {\mathbf{C}}_{{\mathbf{w}}} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} {\mathbf{A}}} \right)}}{n}} \\ = {\frac{{{\text{trace}}\left( {{\mathbf{C}}_{{\mathbf{w}}} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} {\mathbf{AA}}^{T} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} } \right)}}{n}} \\ = {\frac{{{\text{trace}}\left( {{\mathbf{C}}_{{\mathbf{w}}} \left( {{\mathbf{AA}}^{T} } \right)^{ - 1} } \right)}}{n}}. \\ \end{aligned} $$
(25)

Or by substituting Eq. 24 in Eq. 25 we have:

$$ \bar{\sigma }_{\text{Areoid}}^{2} = 68.23{\frac{{{\text{trace}}\left( {\left( {{\mathbf{AA}}^{T} } \right)^{ - 1} } \right)}}{n}}. $$
(26)

Having computed the average accuracy of the Areoid the accuracy of the potential value of the Areoid can be computed as follows:

$$ \sigma_{{W_{0} }} = \pm {\frac{\text{GM}}{{R^{2} }}}\bar{\sigma }_{\text{Areoid}} . $$
(27)

4 Numerical Results of the Areoid Computations

In this section we are going to present the results of our numerical computations for the gravity potential value of Areoid and its shape, according to the methodology explained in previous sections. Figure 4 shows the shape of Mars with respect to its center of mass and the related minimum, maximum, mean and STD. This map has been produced from spherical harmonic expansion of Mars shape model GTM090AA (Smith et al. 1999a) to degree and order 90 from −90 to 90 latitude and 0 to 360 longitude. The resolution of sampling points in direction of latitude and longitude are selected as 2.5° and 2.5°/cosϕ, respectively, in order to generate a grid of sampling points with equal areas, i.e. to give equal weight to all regions.

Fig. 4
figure 4

The shape of Mars with respect to its center of mass based on spherical harmonic expansion to degree and order 90, in kilo meter

Full size image

Using the constrained optimization problem defined in Sect. 2 the Areoid along with its gravity potential value are computed. The gravitational field model of Mars is derived from GGM2BC80 model (lemoine et al. 2001) in terms of spherical harmonics complete to degree and order 80, with parameters which are summarized in Table 1.

Table 1 Mean radius (R), gravitational constant (GM), angular velocity (ω) and zonal spherical harmonic of second order (J 20) of Mars (Lemoine et al. 2001; Folkner et al. 1997)
Full size table

The computational steps involved in the computation of Areoid and its gravity potential value is summarized in Fig. 5.

Fig. 5
figure 5

The flowchart for the computation of geometry and gravity potential value of Areoid

Full size image

Figure 6 shows the computed Areoid with respect to the center of mass of Mars along with a summary of related statistical information. The estimated gravity potential value of the Areoid from the constrained optimization problem is W 0 = (12,654,875 ± 69) (m2/s2). The estimated mean accuracy of the Areoid according to the error analysis procedure explained in Sect. 3 is \( \bar{\sigma}_{\rm Areoid}=\pm19\,{\hbox {m}}\)σAreoid = ±19 m. We also made a plot of the topography of Mars with respect to the computed Areoid, which is shown in Fig. 7.

Fig. 6
figure 6

The computed Areoid with respect to the center of mass of Mars, in kilo meter

Full size image
Fig. 7
figure 7

The topography of Mars with respect to the computed Areoid, in kilo meter

Full size image

5 Numerical Results of the Bi-axial Somigliana-Pizzetti Type Reference Ellipsoid Computed for Mars

According to Grafarend and Ardalan (1999) for determination of a bi-axial Somigliana-Pizzetti type reference ellipsoid from the fundamental parameters {W 0, GM, ω, J 20} following equation must be solved:

$$ \left\{ \begin{gathered} {\frac{\text{GM}}{{\sqrt {a^{2} - b^{2} } }}}{\text{arc}}\cot \left( {{\frac{b}{{\sqrt {a^{2} - b^{2} } }}}} \right) + {\frac{1}{3}}\omega^{2} a^{2} - W_{0} = 0 \hfill \\ {\frac{1}{4}}{\frac{\text{GM}}{{\sqrt {a^{2} - b^{2} } }}}\left( { - 15J_{20} {\frac{{a^{2} }}{{a^{2} - b^{2} }}} + \sqrt 5 } \right)\left[ {\left( {3{\frac{{a^{2} }}{{a^{2} - b^{2} }}} + 1} \right){\text{arc}}\cot \left( {{\frac{b}{{\sqrt {a^{2} - b^{2} } }}}} \right) - 3{\frac{b}{{\sqrt {a^{2} - b^{2} } }}}} \right] \hfill \\ - {\frac{1}{3\sqrt 5 }}\omega^{2} a^{2} = 0. \hfill \\ \end{gathered} \right. $$
(28)

Having computed the potential value of the Areoid W 0 and other three parameters, i.e. {GM, ω, J 20} we used the Eq. 28 and computed a reference ellipsoid of the Somigliana-Pizzetti type for Mars. Table 2 shows the results of our computations. Table 3 offers the parameters of the best fitting bi-axial ellipsoid computed by Duxbury et al. (2002). It is important to note that the reference ellipsoid computed by Duxbury et al. (2002) is resulted from best fitting of an ellipsoid to the shape of Mars and as such is a pure geometrical figure, while our Somigliana-Pizzetti type reference ellipsoid unlike the so-far computed reference ellipsoids, has physical properties, such as being an equipotential surface of the normal gravity field generated by a massive rotational ellipsoid having the mass equal to that of Mars, rotating with the same angular velocity as that of Mars, with rotational symmetry mass distribution (Heiskanen and Mortiz 1967; Vanicek and Krakiwsky 1986; Ardalan and Grafarend 2001). Such a reference ellipsoid could be regarded as the hydrostatic equilibrium figure of a liquid Mars with rotational symmetry mass distribution. The difference between the parameters of the two reference ellipsoids is given in Table 4. As the final result of this study we have plotted of the computed Areoid with respect to the computed reference ellipsoid of Somigliana-Pizzetti type which is shown in Fig. 8.

Table 2 Estimated parameters of bi-axial Somigliana-Pizzetti reference ellipsoid of Mars
Full size table
Table 3 The best fitting bi-axial ellipsoid parameters computed by Duxbury et al. (2002)
Full size table
Table 4 Difference between the parameters of the two reference ellipsoids given in Tables 2 and 3
Full size table
Fig. 8
figure 8

The computed Areoid with respect to computed reference ellipsoid of the Somigliana-Pizzetti type, in meter

Full size image

6 Conclusions

The Areoid is defined as a reference equipotential surface, which best fits to the shape of Mars in least squares sense. Using the shape and gravitational field models of Mars, the shape and gravity potential value of Areoid are computed based on a constrained optimization problem constructed according to the above mentioned definition. Our constrains in the computation of Areoid makes it a hydrostatic equilibrium surface of a liquid Mars. From the solution of the constrained optimization problem the shape of Areoid with accuracy \( \bar{\sigma }_{\text{Areoid}} = \pm 19\,{\text{m}} \) and its gravity potential value W 0 = (12,654,875 ± 69) (m2/s2) are derived. The concept of Somigliana-Pizzetti type reference ellipsoid is applied and based on four fundamental gravity parameters of Mars {W 0, GM, ω, J 20} a bi-axial reference ellipsoid of Somigliana-Pizzetti type with semi-major axis a = (3,395,428 ± 19) (m) and semi-minor axis b = (3,377,678 ± 19) (m) is computed. The Somigliana-Pizzetti type reference ellipsoid unlike the so-far computed reference ellipsoids, which are merely geometric surfaces, has also physical properties. Somigliana-Pizzetti type reference ellipsoid could be the hydrostatic equilibrium figure of a liquid Mars with rotational symmetry mass distribution.

The theory presented in this paper can be regarded as an alternative point of view to planetopotential reference surface (Areoid) and planetary datum (reference ellipsoid) computations. Our numerical tests verified the correctness and capability of the theory and it is evident that by using gravitational and shape models of higher accuracy and resolutions more accurate results could be obtained.