The error is that you assume that the density distribution is "nearly spherically symmetric". It's far enough from spherical symmetry if you want to calculate first-order subleading effects such as the equatorial bulge. If your goal is to compute the deviations of the sea level away from the spherical symmetry (to the first order), it is inconsistent to neglect equally large, first-order corrections to the spherical symmetry on the other side - the source of gravity. In other words, the term $hg$ in your potential is wrong.

Just imagine that the Earth is an ellipsoid with an equatorial bulge, it's not spinning, and there's no water on the surface. What would be the potential on the surface or the potential at a fixed distance from the center of the ellipsoid? You have de facto assumed that in this case, it would be $-GM/R+h(\theta)g$ where $R$ is the fixed Earth's radius (of a spherical matter distribution) and $R+h(\theta)$ is the actual distance of the probe from the origin (center of Earth). However, by this Ansatz, you have only acknowledged the variable distance of the probe from a spherically symmetric source of gravity: you have still neglected the bulge's contribution to the non-sphericity of the gravitational field.

If you include the non-spherically-symmetric correction to the gravitational field of the Earth, $hg$ will approximately change to $hg-hg/2=hg/2$, and correspondingly, the required bulge $\Delta h$ will have to be doubled to compensate for the rotational potential. A heuristic explanation of the factor of $1/2$ is that the true potential above an ellipsoid depends on "something in between" the distance from the center of mass and the distance from the surface. In other words, a "constant potential surface" around an ellipsoidal source of matter is "exactly in between" the actual surface of the ellipsoid and the spherical $R={\rm const}$ surface.

I will try to add more accurate formulae for the gravitational field of the ellipsoid in an updated version of this answer.

**Update: gravitational field of an ellipsoid**

I have numerically verified that the gravitational field of the ellipsoid has exactly the halving effect I sketched above, using a Monte Carlo Mathematica code - to avoid double integrals which might be calculable analytically but I just found it annoying so far.

I took millions of random points inside a prolate ellipsoid with "radii" $(r_x,r_y,r_z)=(0.9,0.9,1.0)$; note that the difference between the two radii is $0.1$. The average value of $1/r$, the inverse distance between the random point of the ellipsoid and a chosen point above the ellipsoid, is $0.05=0.1/2$ smaller if the chosen point is above the equator than if it is above a pole, assuming that the distance from the origin is the same for both chosen points.

**Code:**

```
{xt, yt, zt} = {1.1, 0, 0};
runs = 200000;
totalRinverse = 0;
total = 0;
For[i = 1, i < runs, i++,
x = RandomReal[]*2 - 1;
y = RandomReal[]*2 - 1;
z = RandomReal[]*2 - 1;
inside = x^2/0.81 + y^2/0.81 + z^2 < 1;
total = If[inside, total + 1, total];
totalRinverse =
totalRinverse +
If[inside, 1/Sqrt[(x - xt)^2 + (y - yt)^2 + (z - zt)^2], 0];
]
res1 = N[total/runs / (4 Pi/3/8)]
res2 = N[totalRinverse/runs / (4 Pi/3/8)]
res2/res1
```

**Description**

Use the Mathematica code above: its goal is to calculate a single purely numerical constant because the proportionality of the non-sphericity of the gravitational field to the bulge; mass; Newton's constant is self-evident. The final number that is printed by the code is the average value of $1/r$. If {1.1, 0, 0} is chosen instead of {0, 0, 1.1} at the beginning, the program generates 0.89 instead of 0.94. That proves that the gravitational potential of the ellipsoid behaves as $-GM/R - hg/2$ at distance $R$ from the origin where $h$ is the local height of the surface relatively to the idealized spherical surface.

In the code above, I chose the ellipsoid with radii (0.9, 0.9, 1) which is a prolate spheroid (long, stick-like), unlike the Earth which is close to an oblate spheroid (flat, disk-like). So don't be confused by some signs - they work out OK.

**Bonus from Isaac**

Mariano C. has pointed out the following solution by a rather well-known author:

http://books.google.com/books?id=ySYULc7VEwsC&lpg=PP1&dq=principia%20mathematica&pg=PA424#v=onepage&q&f=false

There are three misconceptions that I can see in your reasoning.

The poles are the only places on Earth where you are **not** accelerating due to the Earth's rotation, so you have that backwards.

You seem to think that the normal force has to be the same at the poles and the equator, which isn't true. Finding the normal force at the poles does not give you the normal force at the equator.

The forces involved are vectors, not scalars, and they are not collinear (except for the special case of the equator). The gravity and normal force are approximately collinear with the Earth's radius everywhere, but the centripetal (or centrifugal) force is not; it points towards (away from) the axis of rotation. So you need to do some vector math/trigonometry to get the actual values.

You seem to be struggling with the distinction between the centripetal force and the centrifugal force. It seems like you have the right idea, but its hard to tell due to the other issues. Let me try to explain what those are.

The centrifugal force is a "fictitious force" (meaning there is no object that exerts this force) that appears to arise in a rotating coordinate system; considering the centrifugal force in a rotating coordinate system maintains the usefulness of Newton's 2nd law.

The centripetal force is the required amount & direction of force that the net force on an object must satisfy in order for the object to move in a circle of a certain radius at a certain speed.

## Best Answer

Equatorial bulging of a planet is caused by the combination of gravity and centrifugal force. To show this I will first make a few assumptions:

A small volume, $dV$, experiences two volumetric accelerations, namely gravitational and centrifugal, and normal forces by the neighboring liquid in therms pressure. The sum of all accelerations on $dV$ should add up to zero to comply with the second assumption (the centrifugal acceleration already accounts for the fact that the reference frame is rotating). At any point on the surface there is a constant pressure, because above it there would be the vacuum of space. This means that the neighboring liquid, also at the surface, has the same pressure and therefore can not exert any force on each other in the plane op the surface. The only direction that liquid can exert force on each other is in the normal direction to the surface. However the sum of all accelerations still should add up to zero and therefore the sum of the gravitational and centrifugal acceleration should also point in the normal direction of the surface.

The magnitude of the gravitational acceleration, $a_g$, is defined by assumption three and its direction is always radially inwards. The magnitude of the centrifugal acceleration, $a_c$, is equal to:

$$ a_c = \omega^2 \sin\phi\ \|\vec{x}\|, $$

where $\phi$ is equal to $\pi/2$ minus the latitude; its direction is always parallel to the plane of the equator and its line of action always goes through the axis of rotation. These accelerations are illustrated in the figure below.

For the next part I will define local unit vectors $\vec{e}_r$ and $\vec{e}_t$, where $\vec{e}_r$ points into the local radial outwards direction and $\vec{e}_t$ is perpendicular to it, lies in the plane spanned by the axis of rotation and $\vec{x}$ and faces the direction closest to the equator. The direction of vectors also correspond with the grey vectors in the figure above. Using these unit vectors, the vector sum of the gravitational and centrifugal acceleration can be written as

$$ \vec{a}_g + \vec{a}_c = \vec{e}_r \left(\omega^2 \sin^2\!\phi\ \|\vec{x}\| - G\frac{M}{\|\vec{x}\|^2}\right) + \vec{e}_t\ \omega^2 \sin\phi\ \cos\phi\ \|\vec{x}\|. $$

If there would be no bulging then the normal vector should always point radial outwards. However the normal vector has to point in the opposite direction as the equation above, which means that for $\omega>0$ it will not point in the same direction as $-\vec{e}_r$ for all values of $\phi$. This means that the surface will have a small slope, $\alpha$, relative to $\vec{e}_t$

$$ \alpha = \tan^{-1}\left(\frac{\omega^2 \sin\phi\ \cos\phi\ \|\vec{x}\|}{G\frac{M}{\|\vec{x}\|^2} - \omega^2 \sin^2\!\phi\ \|\vec{x}\|}\right). $$

A slope means a change of height, and thus radius, when displacing horizontally. To simplify the expression, $r$ will substitute $\|\vec{x}\|$. For a slope $\alpha$ the change of the radius, $dr$, for a small change in $\phi$, $d\phi$, will be equal to:

$$ dr = \tan\alpha\ r\ d\phi. $$

By substituting in the equation for $\alpha$ the following differential equation can be obtained

$$ \frac{dr}{d\phi} = \frac{\omega^2 \sin\phi\ \cos\phi\ r^2}{G\frac{M}{r^2} - \omega^2 \sin^2\!\phi\ r} . $$

When $\phi$ is equal to $0$ or $\frac{\pi}{2}$, the poles and the equator respectively, this equation will be zero, however for any value in between, it will be positive, since when denominator would become negative it would mean that the centrifugal force will be bigger than gravity and the liquid would be flung into space. So this planet would have the smallest radius near the poles after which the radius will increase with $\phi$ until you reach the equator.