Using *Mathematica* to differentiate the potential — after verifying that it satisfies Laplace’s equation — I found that at the surface, where

$$z^2 = (1-e^2)(a^2-r^2),$$

the exact components of the gravitational force simplify to

$$F_r\big|_\text{surface} = -C_r\,\frac{r}{a}\,\frac{GMm}{a^2}$$

$$F_z\big|_\text{surface} = -C_z\,\frac{z}{a}\,\frac{GMm}{a^2},$$

where $C_r$ and $C_z$ are dimensionless constants that depend on the eccentricity $e$ of the oblate spheroid:

$$C_r = \frac32 \frac{\sin^{-1}e-e\sqrt{1-e^2}}{e^3} = 1 + \frac{3}{10}e^2+O(e^4)$$

$$C_z = 3\,\frac{e-\sqrt{1-e^2}\sin^{-1}e}{e^3} = 1 + \frac25e^2+O(e^4).$$

Note that $r\equiv\sqrt{x^2+y^2}$ is the *cylindrical* radial coordinate.

No approximations were involved in this calculation. The point of the series expansions is to show that the expected result follows in the spherical case when $e=0$.

In the opposite limit $e\to 1$, one has $C_r \to 3\pi/4$ and $C_z \to 3$.

The following graph shows the dependence of these constants on the eccentricity. The blue curve is $C_r$ and the gold curve is $C_z$.

I have no idea whether these results are previously known.

[![enter image description here][1]][1]


  [1]: https://i.sstatic.net/rvvM5.png