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
- Are there any other rivers in the world for which this is true?
The Mekong, at least after it has left the Tibetan plateau.
Location Latitude Elevation (m) Radial distance (km)
Source 33° 42' 30" 5224 6376.8
Manwan Reservoir 24° 45' 15.5" 997 6375.4
Ruak River mouth 20° 21' 16" 336 6375.9
Mekong delta 9° 27' 30" 0 6377.6
While the Mekong does flow "downhill" (toward the center of the Earth) at the start thanks to that four kilometer drop in altitude across the Tibetan plateau, it's pretty much all uphill from the Manwan Reservoir on.
- If one uses the actual gravitational potential for a non-rotating oblate spheroid, is the Gulf of Mexico at a higher gravitational potential than Minnesota?
There's a nice what-if scenario at esri, the Environmental Systems Research Institute, that asks and answers this question: If the Earth Stood Still: Modeling the absence of centrifugal force. The scenario investigates what would happen to the Earth's waters if the Earth somehow stopped rotating over the course of a few decades. That's a long enough span of time that the waters would have a chance to adjust, but far too short a span of time for isostasy to readjust the shape of the Earth.
Not only would the Mississippi stop flowing southward, the waters of the oceans would flow away from the equator. This would leave a globe-spanning equatorial bulge of land surrounded by two polar oceans. The end result:
Source: https://www.esri.com/news/arcuser/0610/nospin.html
Best Answer
TL;DR version: There are three reasons. In order of magnitude,
The poles are closer to the center of the Earth due to the equatorial bulge. This strengthens gravitation at the poles and weakens it at the equator.
The equatorial bulge modifies how the Earth the gravitates. This weakens gravitation at the poles and strengthens it at the equator.
The Earth is rotating, so an Earth-bound observer sees a centrifugal force. This has no effect at the poles and weakens gravitation at the equator.
Let's see how the two explanations in the question compare to observation. The following table compares what a spherical gravity model less centrifugal acceleration predicts for gravitational acceleration at sea level at the equator ($g_{\text{eq}}$) and the north pole ($g_{\text{p}}$) versus the values computed using the well-established Somigliana gravity formula $g = g_{\text{eq}}(1+\kappa \sin^2\lambda)/\sqrt{1-e^2\sin^2 \lambda}$.
$\begin{matrix} \text{Quantity} & GM/r^2 & r\omega^2 & \text{Total} & \text{Somigliana} & \text{Error} \\ g_\text{eq} & 9.79828 & -0.03392 & 9.76436 & 9.78033 & -0.01596 \\ g_\text{p} & 9.86431 & 0 & 9.86431 & 9.83219 & \phantom{-}0.03213 \\ g_\text{p} - g_\text{eq} & 0.06604 & \phantom{-}0.03392 & 0.09995 & 0.05186 & \phantom{-}0.04809 \end{matrix}$
This simple model works in a qualitative sense. It shows that gravitation at the north pole is higher than at the equator. Quantitatively, this simple model is not very good. It considerably overstates the difference between gravitation at the north pole versus the equator, almost by a factor of two.
The problem is that this simple model does not account for the gravitational influence of the equatorial bulge. A simple way to think of that bulge is that it adds positive mass at the equator but adds negative mass at the poles, for a zero net change in mass. The negative mass at the pole will reduce gravitation in the vicinity of the pole, while the positive mass at the equator will increase equatorial gravitation. That's exactly what the doctor ordered.
Mathematically, what that moving around of masses does is to create a quadrupole moment in the Earth's gravity field. Without going into the details of spherical harmonics, this adds a term equal to $3 J_2 \frac {GMa^2}{r^4}\left(\frac 3 2 \cos^2 \lambda - 1\right)$ to the gravitational force, where $\lambda$ is the geocentric latitude and $J_2$ is the Earth's second dynamic form. Adding this quadrupole term to the above table yields the following:
$\begin{matrix} \text{Quantity} & GM/r^2 & r\omega^2 & J_2\,\text{term} & \text{Total} & \text{Somigliana} & \text{Error} \\ g_\text{eq} & 9.79828 & -0.03392 & \phantom{-}0.01591 & 9.78027 & 9.78033 & -0.00005 \\ g_\text{p} & 9.86431 & 0 & -0.03225 & 9.83206 & 9.83219 & -0.00013 \\ g_\text{p} - g_\text{eq} & 0.06604 & \phantom{-}0.03392 & -0.04817 & 0.05179 & 0.05186 & -0.00007 \end{matrix}$
This simple addition of the quadrupole now makes for a very nice match.
The numbers I used in the above:
$\mu_E = 398600.0982\,\text{km}^3/\text{s}^2$, the Earth's gravitational parameter less the atmospheric contribution.
$R_\text{eq} = 6378.13672\,\text{km}$, the Earth's equatorial radius (mean tide value).
$1/f = 298.25231$, the Earth's flattening (mean tide value).
$\omega = 7.292115855 \times 10^{-5}\,\text{rad}/\text{s}$, the Earth's rotation rate.
$J_2 = 0.0010826359$, the Earth's second dynamic form factor.
$g_{\text{eq}} = 9.7803267714\,\text{m}/\text{s}^2$, gravitation at sea level at the equator.
$\kappa = 0.00193185138639$, which reflects the observed difference between gravitation at the equator versus the poles.
$e^2 = 0.00669437999013$, the square of the eccentricity of the figure of the Earth.
These values are mostly from Groten, "Fundamental parameters and current (2004) best estimates of the parameters of common relevance to astronomy, geodesy, and geodynamics." Journal of Geodesy, 77:10-11 724-797 (2004), with the standard gravitational parameter modified to exclude the mass of the atmosphere. The Earth's atmosphere has a gravitational effect on the Moon and on satellites, but not so much on people standing on the surface of the Earth.