A triaxial reference ellipsoid for the Earth

We present a new, physically motivated triaxial reference ellipsoid for the Earth. It is an equipotential surface in the gravity field and closely approximates the geoid, akin to the conventional reference ellipsoid of revolution. According to Burša and Fialová (Studia Geophysica et Geodaetica 37(1):1–13, 1993), the triaxial reference ellipsoid is uniquely, but not exclusively, specified by the body’s total mass, the dynamic form factors of polar and equatorial flattening, the longitude of the equatorial major axis, the rotation rate, and the designated surface potential. We model the gravity field using triaxial ellipsoidal harmonics. While they are rarely considered practical for near-spherical planets, we leverage an intrinsic property that ellipsoidal harmonics yield an exact expression for the constant potential on a triaxial ellipsoid. A practical procedure is proposed to solve for the ellipsoidal parameters that converge iteratively to fulfill the exact condition of equipotentiality. We present the solution for the Earth Gravitational Model 2008.


Introduction
Earth has a mean radius of 6371 km. The polar semi-axis is shorter than the equatorial semi-axis by some 20 km. That is, to a first-order approximation in terms of polar flattening, the Earth is an oblate biaxial ellipsoid or spheroid. A flattened figure is common for large rotating planetary bodies and a manifestation of their states of (near) hydrostatic equilibrium. Sea level on Earth closely approximates a gravity equipotential surface called the geoid. The maximum deviations between these surfaces is ∼ 1 m. An oblate biaxial ellipsoid can be used in turn to approximate the geoid: The maximum vertical deviation of the geoid from this ellipsoid is ∼ 100 m. This ellipsoid is widely referred to as a reference ellipsoid (Heiskanen and Moritz 1967;Moritz 2000;Torge and Müller 2012;Jekeli 2007;Grafarend et al. 2010). In this work, we shall refer to this conventional reference ellipsoid as the biaxial reference ellipsoid. B Xuanyu Hu xuanyuhu@gmail.com 1 Gravity described in the rotating body-fixed frame consists of two distinct components, the gravitation due to body's mass distribution and the centrifugal acceleration due to its rotation (Pizzetti 1894;Heiskanen and Moritz 1967). We denote the gravity potential as where V is the gravitational potential at a field point r due to a distribution of differential body mass dM located at r, G is the gravitational constant, and is the centrifugal potential at r where ω = ωω is the rotational velocity of the body. For example, the expression in the Cartesian coordinate system whose z axis coincides withω, the unit vector of the rotation axis, is = ω 2 (x 2 + y 2 )/2.

Biaxial reference ellipsoid of Earth
An equipotential surface is a level or surface on which the gravity potential is constant. The definition of the biaxial reference ellipsoid is conveniently and most concisely expressed in terms of biaxial ellipsoidal or spheroidal harmonics (Pizzetti 1894;Heiskanen and Moritz 1967;Thong and Grafarend 1989;Jekeli 2007;Grafarend et al. 2010;Sebera et al. 2016). We denote the normal gravity potential due to this biaxial ellipsoid in the following way U (u, β) = V + = ∞ n=0 Q n (iu/E) Q n (ib/E) A n P n (sin β) where the notation conforms to that of Heiskanen and Moritz (1967). u and β are the biaxial ellipsoidal coordinates of the semiminor axis and reduced latitude, respectively. The expression does not involve the third coordinate of longitude because of rotational symmetry. On the biaxial reference ellipsoid, we have u = b. E = √ a 2 − b 2 is the focal length, with a being the semimajor axis. P n and Q n are the Legendre functions of the first and second kind, respectively, of degree n. i is the imaginary unit. A n are coefficients quantifying the degree-wise variation of the gravitational field. The second term on the right-hand side is the centrifugal potential.
We suppose that U on the biaxial reference ellipsoid reduces to a constant, namely Taking note that cos 2 β = 2 3 [1 − P 2 (sin β)], this condition yields a closed-form expression for the normal gravitational potential generated by the ellipsoid. Specifically, all A n other than those of n = 0 and n = 2 vanish, Simply put, the reduction of (3) requires that any latitudinal variation of the normal gravitational potential over the biaxial reference ellipsoid must be canceled out by the centrifugal component, thereby removing β. As a result, the resolution of the normal gravitational field does not exceed degree two.
The gravitational field is most often modeled as a series of spherical harmonics (Hobson 1931;Heiskanen and Moritz 1967;Jekeli 2007), where r , ϕ, and λ are the spherical coordinates of radius, latitude, and longitude, respectively. P nl is the associated Legendre function of degree n and order l. r 0 is a reference radius that can be chosen as that of the circumscribing sphere of the body mass. C nl and S nl are Stokes' coefficients and C 0,0 = 1 and S n0 = 0 for all n. Note that P nl and C nl , S nl are normalized in practice, e.g.,P nl = (2 − δ l0 )(2n + 1) (n−l)! (n+l)! P nl with δ l0 being the Kronecker delta. The series is truncated at some finite degree N . The Stokes' coefficients are determined from satellite tracking, ground gravity measurements, etc.
The biaxial reference ellipsoid can be specified with the aid of such a model. For example, the polar oblateness of Earth's mass distribution, indicated by J 2 = −C 2,0 , was adopted alongside a (i.e., r 0 ), G M, and ω for defining Geodetic Reference System 1980(GRS 1980) (Moritz 2000 as well as the descendent World Geodetic System 1984 (WGS 84) (Defense Mapping Agency 1991). Here, the physical parameter J 2 is equivalent to the geometric flattening, f = (a−b)/a, by virtue of Clairaut's theorem. The theory above is from an authoritative account of the topic by Heiskanen and Moritz (1967).

Triaxial figure: Burša and Fialová's approach
Earth also exhibits an equatorial flattening. It is a secondorder effect but nonetheless easily perceptible given ever more precise geodetic measurements. The triaxiality of Earth is a classic topic in geodesy, whose signature in the early days was a subject of debate (Herschel 1879). A brief history can be found in Hollis (1906) and Heiskanen (1962). The equatorial flattening indicates a degree of deviation of the body from the perfect isostatic equilibrium, or mass anomalies at depths, e.g., the core-mantle boundary (Heiskanen 1962).
Earth's triaxiality has been investigated historically via astro-geodetic or gravimetric measurements (Clarke 1861; Heiskanen 1928), and since the advent of the space era, (together with) observed satellite motions (Kaula 1959;Izsak 1961;Kozai 1961). As an approximate equipotential surface, this reference figure, whether a biaxial or triaxial ellipsoid, can be determined as best fits to the geoid derived from the gravitational field model (Burša 1970;Burša and Sima 1980;Tserklevych et al. 2016;Panou et al. 2020;Soler and Han 2020). This is the same principle as determining the triaxial dimensions of other planetary bodies (Smith et al. 1999;Iz et al. 2011). On the other hand, Burša and Fialová (1993) (hereafter B&F 1993) presented a definition and solution for the Earth's triaxial figure, which is a generalization of the biaxial case as described in Sect. 1.1. The approach is based on spherical harmonics and applies the same rationale as the condition of (3). In addition to a, G M, J 2 , ω, the sectoral coefficients, C 2,2 , S 2,2 , are needed to specify the equatorial flattening. While the equatorial flattening itself can be measured by one parameter, which we denote by  B&F (1993) purpose of clarity, the physical constants of the biaxial and triaxial reference ellipsoids are compared in Table 1. Spherical harmonics in the form of Eq. (5) do not directly provide a closed-form expression for the constant potential on a triaxial ellipsoid (or a biaxial one). B&F (1993) first expanded the radius coordinate, r , of a point on the ellipsoid as another infinite spherical harmonic series with respect to r 0 , i.e., r (ϕ, λ) = r 0 + ∞ k=0 k l=0 (a kl cos lλ + b kl sin lλ)P kl (sin ϕ). The coefficients, a kl , b kl , are themselves power series of the polar and equatorial flattening as well as their product; those smaller than 10 −11 were omitted. Then, r (ϕ, λ) was incorporated into an analogue of Eq. (3) to solve for the triaxial figure on which the gravity potential reduces to a constant.

Statement of problem
In this work, we present an alternative solution for the triaxial reference ellipsoid of Earth. Hereafter, we reserve the term "ellipsoid" (and "ellipsoidal") exclusively for the triaxial figure, while the biaxial case will always be designated explicitly. The solution is based on triaxial ellipsoidal harmonics (EHs) and it has closed form (Morera 1894;Caputo 1967;Hu 2017). The closure results from the fact that the constant potential on an ellipsoid can be represented by a single term of degree zero, similar to the degree-zero spherical and biaxial EH that indicate constant potential on a sphere and a biaxial ellipsoid, respectively (Hu and Jekeli 2015). Consequently, we may rigorously impose the same condition as Eq. (3) for the triaxial case. Zagrebin (1973) formulated the Earth's normal gravity field in terms of EHs as generated by a geoid-fitting ellipsoid, i.e., of known dimensions. Hu (2017) presented an iterative approach to determine the dimensions of this reference ellipsoid for small, irregular-shaped extraterrestrial bodies, such as the Martian moon, Phobos (Hu et al. 2020). The validity of the algorithm rests on the assumption that the body's polar and equatorial flattenings are comparable. It is hence inapplicable to the Earth whose equatorial flattening ∼ 10 −6 is far smaller than the polar flattening of ∼ 10 −3 . In other words, the triaxiality of the Earth is not prominent. (By contrast, the polar and equatorial flattening of Phobos are both on the order of 0.1.) Here we present an efficient, alternative method, which finds the ellipsoidal semi-axes by searching for the equipotential level along the respective gravity vectors. Because the flattenings are never explicitly invoked, any assumption or restriction concerning their magnitude is obviated. The method presented here is complete and standalone: It yields the same results as in Hu (2017) for Phobos, but also applicable to the Earth and other planetary bodies with indistinctly triaxial figures.
The effort of developing an EH gravity field model for the Earth dates back several decades, when the triaxiality of the figure was no longer in doubt (see Walter 1970;Madden 1970, and reference therein). Since then, however, the practical interest has not quite materialized and is limited to small bodies, following the work by Garmier and Barriot (2001). Granted, a biaxial ellipsoid remains an intuitive and apposite reference for the Earth in near hydrostatic equilibrium. With the measurement precision nowadays far exceeding the (in)distinctness of the equatorial flattening, the triaxial ellipsoid has received renewed attention for being a more accurate and natural reference figure (Panou et al. 2020;Soler and Han 2020). This work exploits the innate advantages of the EHs and is among the first to offer a practical solution on the topic. It is intended to lay a piece of the groundwork for a possible, generalized approach for large and small bodies alike.
The discussion proceeds as follows. In Sect. 2, the basics of triaxial EHs are given for modeling the gravitational field and the centrifugal effect. In Sect. 3, the definition of the reference ellipsoid based on EHs is reviewed and the methodology of reference ellipsoid determination is given, with some more cumbersome derivations arranged in the appendices. We then present a benchmark solution in Sect. 4, based on the same parameters used by B&F (1993), as a means of validation. In Sect. 5, the solution compatible with the Earth Gravitational Model 2008 (EGM 2008) will be given (Pavlis et al. 2012). The geoid undulations with respect to the triaxial reference ellipsoid are presented. We also discuss our results in comparison with a least-squares solution by Panou et al. (2020).

Gravity field in terms of triaxial ellipsoidal harmonics
The gravitational potential can be expressed as an EH series as follows (Hobson 1931;Garmier and Barriot 2001;Hu 2012), where E nm (λ) and F nm (λ) are Lamé functions of the first and second kind, respectively, for degree n and order m. The infinite series is always convergent for ρ > a. ρ,μ,ν are the ellipsoidal coordinates, being the three distinct roots of the following cubic equation of t 2 , for certain real constants k > h > 0. It is imposed that Then, ρ specifies (the semimajor axis of) an ellipsoid, while μ, ν specify a onesheet and a two-sheet hyperboloids, respectively, all with focal lengths h and k. The long, intermediate, and short axes of the ellipsoid are aligned with x, y, and z axes, respectively. The Cartesian coordinates are related to the ellipsoidal via Hence, while they have the dimension of distance, μ, ν play the role of the angular arguments (e.g., latitude and longitude). There is a sign ambiguity, which, however, is not an obstacle in practice since the ellipsoidal coordinates are never used exclusively and the above expressions are not often applied for coordinate transformation. The solution of (6) refers to an ellipsoid of semi-axes a, b = √ a 2 − h 2 , and c = √ a 2 − k 2 . c nm are coefficients characterizing the gravitational field of the body associated with the ellipsoid.
There are four classes of Lamé functions, to one of which E nm belongs. The present discussion is only concerned with functions of class "K" of the following form (Hobson 1931), with d = n/2 . There are d + 1 such functions for each degree. The coefficients a (i) nm are to be determined for given h, k so that E nm satisfies the Lamé equation (Hobson 1931). The scale of a Lamé function is arbitrary, however. E nm (μ)E nm (ν) has the same role as the surface spherical harmonics, P nl (sin ϕ) cos lλ sin lλ . Note that the order m differs from the counterpart l in general. The class K is distinguished by the correspondence between the associated surface harmonic K nm (μ)K nm (ν) and the surface spherical harmonic P n,2l (sin ϕ) cos 2lλ. The Lamé function of the second kind is analogous to the radial function of r −n−1 of spherical harmonics and Q nl (iu/E) of biaxial EHs. E nm (ρ)E nm (μ)E nm (ν) of any class represents an alternative solution to Laplace's equation. It is commonly known as a solid EH, in analogy to the solid spherical harmonic of r n P nl (sin ϕ) cos lλ sin lλ (Heiskanen and Moritz 1967). A solid harmonic of degree n corresponds to a certain linear combination of mass density moments of the order r n .

Sectoral and zonal EHs of the second degree
The evaluation of EHs is far more complicated than that of spherical harmonics or biaxial EHs (Garmier and Barriot 2001;Hu and Jekeli 2015;Reimond and Baur 2016). Fortunately, in this work only three EHs are of concern, namely, of degree zero and two and all belonging to class K. The degreezero EH is a constant, e.g., K 01 = 1 is usually assumed. We simplify the notation for a second-degree EH as The coefficient a m is given by where s m is a root of the following quadratic equation Hobson (1931),

Example: Earth
It is now illustrative to consider an example for the Earth.
Taking the lengths of the three semi-axes as (Panou et al. 2020), we obtain, It can be shown that the corresponding solid EHs, K 2,m = K 2,m (ρ)K 2,m (μ)K 2,m (ν), can be expressed in terms of Cartesian coordinates as (Appendix A), where p x , p y , p z , p 0 are coefficients. Those for a, b, c given above are The expressions are rescaled such that p x = ±1.

Properties
Two things are noteworthy. First, K 2,1 is formally similar to the sectoral solid spherical harmonic, r 2 P 2,2 (sin ϕ) cos 2λ = 3(x 2 − y 2 ) and K 2,2 to the zonal harmonic, r 2 P 2 (sin ϕ) = z 2 − x 2 + y 2 /2. Note that the coefficients of the respective (squared) coordinates are similar but not identical in proportion to those of spherical harmonics. The solid EHs also consist of a constant term, here factored by h 2 , which is absent in spherical harmonics. The variations of the sectoral and zonal EHs over the ellipsoid surface (i.e., for constant ρ) are illustrated in Fig. 1, where the triaxiality is exaggerated compared with Earth's. The patterns are fully analogous to those of the respective spherical harmonics on a sphere. Second, for both K 2,m we note the obvious, necessary condition (in the absence of round-off errors),  (15) is of practical interest to us in this work, as will be illustrated in the following section.

Numerical precision
The computation of EHs is susceptible to numerical instabilities from a low degree, e.g., about 15 in the case of the elongated comet (nucleus) Wirtanen (Garmier and Barriot 2001). One issue, for instance, is that the roots of the Lamé polynomials become increasingly inseparable. Our calculation is based on the 16-digit format (for double-precision floating-point numbers). While the model resolution is only up to degree two, a potential issue resides with the indistinct equatorial flattening of Earth compared with the polar flattening, e.g., h ≈ 1 20 k. To examine the numerical precision of the Lamé polynomials, we make use of the following equality (Hobson 1931), where is the angle between two vectors, r, r , and cos = r · r /(|r||r |). The expression in spherical coordinates is cos = sin ϕ sin ϕ + cos ϕ cos ϕ cos(λ − λ ). The normalization factor is defined as We set r = (a 0 0) and let r trace the equator, the equilatitudes of 30 • and 60 • on the ellipsoid, with the longitude incrementing from 0 through 180 • in all cases. Thus, varies from 0 through 180, 150, and 120 degrees, respectively. Figure 2 shows the discrepancy between the left-and right-hand sides of Eq. (16). The error is at the level of 10 −13 . Because |P 2 (cos )| ∈ [0, 1], this error indicates a loss of 3 digits in precision.

Centrifugal potential
Our goal is to solve for an equipotential triaxial ellipsoid. The centrifugal potential needs to be expressed as EHs in terms of the ellipsoidal coordinates, ρ, μ, ν. Toward this end, we seek to express a solid EH on an ellipsoid in a similar form as the centrifugal potential of x 2 + y 2 . Referring to Eq. (7), on the ellipsoid of semimajor axis a we have z 2 = c 2 1 − x 2 /a 2 − y 2 /b 2 . Consequently, Eq. (14) can be rewritten as follows, It is now clear that the centrifugal potential over the ellipsoid can be expressed by a linear combination of the surface sectoral and zonal EHs plus a constant, i.e., where φ 0 , φ 1 , φ 2 are coefficients that depend on a. Let us contract the result of Eq. (17) as K 2,m = 1 x 2 y 2 p 2,m with p 2,m being a column vector of ( p 0 + p z c 2 p x − p z c 2 /a 2 p y − p z c 2 /b 2 ) T for the corresponding harmonic.
Then, φ 0 , φ 1 , φ 2 , are the solution of the linear system: where the left-hand side corresponds to = (1 x 2 y 2 ) (0 ω 2 /2 ω 2 /2) T , and accordingly, the rows of the matrix on the right-hand side account for the constant, x 2 , and y 2 , respectively. K 2,m is evaluated via Eq. (11); it is never zero for a > k (Hobson 1931); hence, no singularity arises.

Gravitational field: transformation from spherical to EHs
As remarked in Sect. 1.1, the gravitational field model is usually determined in the form of a spherical harmonic series. Therefore, it is required to transform a spherical harmonic model into an EH model. Here we apply an exact transformation method as discussed in Hu (2016). The method is founded on the fact that the series coefficients, spherical, biaxial or triaxial EHs, are uniquely determined by the mass distribution of the body (Jekeli 1988;Hu 2016). Therefore, the transformation can be accomplished by expressing the solid EH as a linear combination of the solid spherical harmonics. We note that the transformation from spherical to the biaxial EHs and the reversion was presented by Jekeli (1988).
Happily, we are concerned with the transformation of coefficients only up to degree two. The physical meaning of the second-degree zonal and sectoral spherical harmonic coefficients is well known, noting the trivial case of C 00 = 1 (Heiskanen and Moritz 1967), where R 2,0 , R 2,2 denote, respectively, the solid zonal and sectoral spherical harmonics. On the other hand, it can be shown that the EH coefficients are related to the body's mass density moments as follows (Hu 2016), where γ 2,m is the normalization factor. The integrand is none other than the solid EH given by Eq. (14). Thus, the physical meaning of c nm is fully analogous to that of C nm . Note that a superficial difference is that C nm is unitless, whereas c nm is defined to be of the same dimension as potential (i.e., G M/a), for the sake of convenience. We may express the second-degree EH coefficients as a linear combination of the C 2,0 , C 2,2 plus a constant. The strategy is the same as for the centrifugal potential (Sect. 2.2). We first express the solid EH as a linear combination of solid spherical harmonics. With the resulting expression incorporated into the integrand of Eq. (21), c 2,m is then related linearly to C 2,0 , C 2,2 . Let the solid EH be, The coefficients, α (Z) , α (S) , α (0) , can be determined in the same way as Eq. (19), This expression is perhaps worth a brief remark. The three equations account for the coefficients for the constant term, x 2 , and y 2 , respectively. Moreover, we note that the term of p z z 2 in K and also the integrand in Eq. (21) is not needed, since p z = −p x − p y as dictated by Laplace's equation (Eq. 15). The columns in the square matrix on the right-hand side correspond to the constant term, the zonal and the sectoral spherical harmonics, respectively; again, the terms of z 2 are redundant and thus omitted. The solutions had best be rearranged explicitly as follows Finally, upon substituting Eq. (22) with α given by (24) into the integrand of (21), we have Such is the direct, exact transformation from spherical harmonic coefficients into zonal and sectoral EH coefficients of degree two. The expressions can also be obtained via the more general numerical procedure as described in Hu (2016).

Criterion of triaxial reference ellipsoid
Once equipped with Eqs. (18) and (25), we are ready to construct the same condition to define and determine the triaxial reference ellipsoid, as by Eq. (3) for the biaxial ellipsoid.
First, let us express the gravity potential due to the ellipsoid in terms of EHs as follows, where the compacted expression of V sums up only three terms associated with c 0,1 , c 2,1 , and c 2,2 . We impose that on the reference ellipsoid, say, of semi-axes a 0 , b 0 , c 0 and corresponding focal lengths h 0 = a 2 0 − b 2 0 , k 0 = a 2 0 − c 2 0 , the gravity potential is constant, which can only hold if the arguments of μ, ν vanish. Be reminded of Eq.
(3) that on the biaxial reference ellipsoid, the latitudinal argument, β, is canceled out when the associated terms in gravitational and centrifugal potentials cancel out. Here, we expect not only latitudinal but also any longitudinal variation given by μ, ν in the gravitational field to be canceled out by the respective centrifugal terms. The criterion is, specifically, Thus, the gravitational field coefficient, c 2,m , and the centrifugal, φ m , must cancel out such that the gravity potential no longer depend on μ, ν. The above expression then reduces to three diagnostics formally equal to zero, with the other being c 0,1 + φ 0 − U 0 .

First-order solution
Suppose we work with an initial ellipsoid, whose semimajor, semi-intermediate, and semiminor axes are given by a, b, c.
In general, this ellipsoid is not equipotential, in which case Eq. (28) does not hold. We now face one last task to determine a new set of parameters, a 0 , b 0 , c 0 , which satisfy Eq. (28) and which correspond to the reference ellipsoid. We note that the gravity field due to an ellipsoid has triaxial symmetry, namely U (x, y, z) = U (|x|, |y|, |z|). Furthermore, to determine three parameters, it suffices to deal with three field points along the respective axes, i.e., U (a, 0, 0), U (0, b, 0), and U (0, 0, c). The deviation of the potentials from U 0 is then given by, where a = a 0 − a, b = b 0 − b, c = c 0 − c measure the deviation of the current ellipsoidal dimensions from those of the target reference ellipsoid. The approximation applies as long as U is small compared with U . The three partial derivatives are associated with the signed magnitude of the "normal" gravity, which is always orthogonal to the ellipsoid. At the respective ellipsoidal extremities, only one coordinate component is nonzero, e.g., γ | x=a = (γ a 0 0) with γ a = (∂U /∂ x)| x=a . The first-order solution of the reference ellipsoid can be explicated as For the sake of completeness of our discussion, the formulas for calculating γ are given in Appendix B. However, for the Earth, this calculation can be bypassed. Instead, a zerothorder approximation suffices in practice, e.g., Note the minus sign must be included.

Practical procedure
Depending on the desired precision, the first-order solution may not always be sufficient. The above solution can be iterated until the result has converged whose variation is within a threshold of tolerance. The solution can be simplified for the Earth (and other nearly spherical bodies) as the evaluation of the normal gravity is unnecessary. We assume the following information is at hand: 1. The base or original gravitational field model is a spherical harmonic series as Eq. (5). Specifically, we will need G M, J 2 = −C 2,0 , J 2 = −C 2,0 , J 2,2 = C 2 2,2 + S 2 2,2 . 2. The rotation rate, ω, of Earth is known. 3. The target potential number, U 0 , or equivalently, the reference radius, R 0 = G M/U 0 .
Note that the six parameters are exactly as prescribed by B&F (1993). The practical procedure of the reference ellipsoid solution consists of the following steps.
Step 1. Start with an initial ellipsoid with semi-axes, a, b, c and prepare the Lamé polynomials of degree two, class K, according to Eqs. (11)-(13). Then, derive the solid EH in terms of the Cartesian coordinates as Eq. (14); the formulas for calculating coefficients p x , p y , p z , p 0 are given in Appendix A.
Step 2. Prepare the gravity field model in terms of EHs up to degree two as given by (26). This consists of two operations: (i). For the gravitational field, the zero-degree coefficient, c 0,1 is given by Eq. (21) while c 2,m are given by (25). Note that, in the latter case, J 2,2 is used in place of C 2,2 , since S 2,2 also contributes to the equatorial flattening.
Step 4. Evaluate the signed magnitude of gravity, γ a , γ b , γ c as are needed in Eq. (31). For a first-order solution, γ = −9.8 m s −2 can be used. Then, Eq. (31) is used to solve for the new parameters for the reference ellipsoid, accurate to the order of U .

Benchmark solution: comparison with Burša and Fialová (1993)
We applied the iterative method described above to solve for the Earth's reference ellipsoid. For illustrative as well as validating purposes, we performed the solution obtained by B&F (1993). The parameters are given in Table 2. Note that the target potential is U 0 = G M/R 0 . We used a = 6380 km, b = 6379 km, c = 6350 km as the starting semi-axes. The starting values are unimportant here, as long as they satisfy a > b > c and are distinct from one another. Setting the criterion to be ε = 1 mm in Eq. (33), the solution converged after 3 iterations. The solution is given in Table 3 alongside the original in B&F (1993). The full numerical output of 16 digits is shown for the present solution. The result is almost in perfect agreement with B&F (1993), whose solution retained terms down to 10 −10 . The two solu-tions are identical for the semimajor axis. The difference in (the inverse of) polar and equatorial flattening occurs only in the last digit. Compared with the corresponding semiminor axes obtained by B&F (1993), the difference is at the level of 10 −3 m, which well corresponds to 10 −10 of the semi-axes.
We performed another solution using a more stringent criterion of ε = 10 −8 m. Here, two more iterations were needed. Figure 3 shows the behavior of the residual potential, given by |c 0,1 +φ 0 −U 0 |, |c 2,1 +φ 1 |, |c 2,2 +φ 2 |, after each iteration. According to Eq. (28), the last two quantities indicate deviation from an equipotential ellipsoid while the first indicates discrepancy in scale; all three should numerically approach zero as the solution is being iteratively refined. The residuals decreased steadily and after five iterations the difference from the target U 0 was less than 10 −8 m 2 s −2 (solid black curve), which is on the order of 10 −16 of U 0 and thus the limit of the machine precision in this work.

Triaxial reference ellipsoid for EGM 2008
As a further demonstration, we performed a solution for the reference ellipsoid compatible with EGM 2008 (Pavlis et al. 2012). A milestone accomplishment, EGM 2008 has assimilated both gravity and altimetric data and since its release served as a benchmark solution for various further model syntheses. 1 The six parameters required for the reference ellipsoid solution are summarized in Table 4. Tidal effects were neglected, i.e., the solution was for the tide-free case. The potential was set to be U 0 = 62636851.7146 m 2 s −2 , which corresponds to the number for WGS 84 (biaxial) reference ellipsoid. The resulting parameters of the reference ellipsoid are given in Table 5.
It is informative to compare the result with another solution by Panou et al. (2020). They derived the reference ellipsoids for various gravitational field models via a geometric method, i.e., as the best-fitting shapes to the corresponding geoid models. In comparison, the solution here is exact for the designated parameters; it is blind to the geoid model  (beyond degree two) and involves no minimization of distances between the geoid and the reference ellipsoid.
Our results are in reasonable agreement with those obtained by Panou et al. (2020). The difference of the semimajor axis is 2 cm, within 1 σ uncertainty of their leastsquares solution. The discrepancy of other two axes is larger, at about 7 and 9 cm, respectively, but still within 2 σ . Hence, even for a preliminary treatment, the physical solution here is consistent with the geometric solution. Figure 4 shows the geoid undulations of EGM 2008 with respect to the WGS84 biaxial reference ellipsoid (top), the derived triaxial ellipsoid in Table 5 (middle), and the difference between them. The gridded undulations with respect to WGS 84 were obtained from Office of Geomatics, National Geospatial-Intelligence Agency (https://earth-info.nga.mil/ index.php?action=home) and interpolated at an interval of 1 • in both longitude and geodetic latitude. The geoid undulations vary between − 106.52 and 84.61 m. We evaluated the area-weighted root mean square (WRMS) as follows, where N i is the geoid undulation at the ith grid point, A i is the surface area of the corresponding grid cell and A ≈ i A i . We used spherical approximation A ≈ cos β β λ; the discrepancy from the ellipsoidal expression is on the order of f 2 (with f being the polar flattening) hence immaterial. The WRMS of the undulations with respect to the WGS84 biaxial ellipsoid is 30.60 m.
We re-measured the undulations of the geoid grid points with respect to the triaxial reference ellipsoid along the ellipsoidal normal via the method by Lin and Wang (1995) adapted for the triaxial case (see also Ligas 2012). Note that this ellipsoid is rotated by λ 0 = −14.93 • with respect to the prime meridian, which must be accounted for in the measurements. In this case, the maximum and minimum were moderated to 69.67 and − 72.31 m, respectively. The WRMS was also reduced to 24.71 m. Thus, the RMS misfit is reduced by more than 19%. The behavior of these statistics is similar to that reported by Panou et al. (2020). Table 6 provides a summary of the comparison. The difference of a given statistic therein is calculated as which always indicates an improvement of fitting. The difference between the two sets of geoid undulations displays a sectoral pattern of the second degree, more prominent at lower latitudes (bottom panel of Fig. 4). The pattern is conceivable, as it reflects the presence of Earth's equatorial flattening, a feature unaccounted for in the biaxial case (thus present in the geoid undulations). On the other hand, it is worth noting that the semiminor axis of the triaxial reference ellipsoid, c 0 = 6356752.33 m (right column in Table 5), differs from that of the WGS 84 biaxial ellipsoid, 6356752.31 m, by only 2 cm, which explains the vanishing difference near the poles.
The Supplementary Information contains the data file of geoid undulations with respect to the biaxial and triaxial reference ellipsoids as displayed in Fig. 4. The results are also presented in pole-centered, orthographic map projections, which incur less area distortion in the polar regions.

Conclusion and discussion
We have presented a new method to determine the Earth's triaxial reference ellipsoid from its gravity field model. The condition of equipotentiality is formulated by a finite EH series up to degree two, in analogy with the well-known theory for the biaxial case. The formulation was previously proposed for irregular-shaped Phobos as described in Hu (2017). The algorithms differ from the original in two aspects. First, the present solution approaches the equipotential reference ellipsoid via iterative corrections of the semi-axes along the respective normal gravity vectors, whereas the method in Hu (2017) does not involve evaluation of the normal gravity. The use of a constant gravity as an approximation is fully justified in the case of Earth. Second, the transformation of the gravitational field from spherical harmonics into triaxial EHs is founded on the uniqueness of mass distribution and therefore the exactness of the respective expressions (Hu 2016). This is more efficient than the approach in Hu (2017) based on the solution of a triaxial ellipsoidal boundary value problem, with gridded boundary values given by spherical harmonics.
In the case of Phobos, the two methods are numerically identical. But, the method in Hu (2017) cannot be applied (i.e., won't converge) for the Earth. Compared with the method by B&F (1993) involving an infinite spherical harmonic series, an obvious advantage here is that the solution is exact, i.e., incurs no truncation errors. The exactness is a consequence or expression of the fact that the EHs up to degree two are sufficient to capture the ellipsoidicity of the gravity field (equipotentials), which exhibits an infinite bandwidth in the spherical, or any non-ellipsoidal, spectrum. A potential deterrence in practice is that the EHs are computationally more complicated than spherical or biaxial ellipsoidal harmonics. However, the procedure only requires evaluation of three terms up to degree two and of the same class, for which the Lamé polynomials are but quadratic functions. The computation is thus deceptively straightforward. In a benchmark case, we obtained results in remarkable agreement with those by B&F (1993) for the given truncation errors at the level of millimeters.
Our solution also differs in nature from least-squares solutions, such as those by Tserklevych et al. (2016) and Panou et al. (2020),which adjust not only dimensions but also origin and orientation of the ellipsoid so as to best fit the geoid. For instance, the longitude of the equatorial elongation, λ 0 , is directly estimated by the least-squares approach (Panou et al. 2020). In comparison, the same parameter is fully specified by the sectoral Stokes' coefficients, C 2,2 , S 2,2 , in the present solution. The ellipsoidal origin is also fixed at the center-ofmass of the Earth, as is (always) stipulated in the gravitational field model in order to suppress the degree-one coefficients. Despite the essential difference, our result for the EGM2008 differs from that of Panou et al. (2020) by no more than a decimeter and within the 2σ uncertainty. These comparisons indicate that the result here is a valid alternative and can be directly used as a triaxial reference ellipsoid for the Earth. The procedure should also be practical for solutions compatible with other gravitational field models or for other planetary bodies.
The present solution is for the tide-free case. At least the permanent tidal effect can be addressed without extra effort. All that is required is to use C 2,0 for the "zero-tide" case, namely one that has been corrected for permanent tide (see IERS Conventions, Petit and Luzum 2010). Then, the procedure outlined in Sect. 3.3 is applied without any alteration.
respect to the ellipsoidal coordinates are The issue here is only with the derivative of the Lamé function of the second kind. We restrict ourselves to n = 0, 2, class K, whereas the derivative of the first kind is simply, The calculation of higher-, odd-degree harmonics, and those of other classes, is more complicated, e.g., the sign ambiguity of the Lamé functions must be carefully addressed (Garmier and Barriot 2001). For the present discussion concerning only normal gravity, the above formulas suffice. At the ellipsoidal extremities the evaluation is much simpler. At (x, y, z) = (a, 0, 0), we have ρ = a, μ = k, ν = h; thus, only ∂ρ/∂ x in Eq. (B3) is nonzero. Similarly, we have ρ = b, μ = k, ν = 0 at (x, y, z) = (0, b, 0) and ρ = c, μ = h, ν = 0 at (x, y, z) = (0, 0, c). Hence, only ∂ρ/∂ y and ∂ρ/∂z remain in the respective cases. Thus, ∂ V /∂μ, ∂ V /∂ν need not be calculated. It confirms that only one coordinate component is nonzero along the corresponding axis, e.g., γ | x=a = (γ a 0 0) as in Eq. (31).
The above formulas are not needed for deriving Earth's reference ellipsoid, and an approximate value of − 9.8 m s −2 can be used for the gravity without any perceptible difference in performance. However, in the case of a distinctly triaxial object (e.g., Phobos), the formulas are more accurate and ensure better convergence of the solution.