Using Mathematica to differentiate the exterior potential — after first 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 for the surface gravity are previously known.