Why is the Earth so fat?
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
I) Flatness. Here we would like to calculated analytically Lubos Motl's solution to the first order in the flatness parameter $f$,
$$0~<f~:=~1-\frac{b}{a}~\approx~\frac{a}{b}-1 ~\ll~ 1, \tag{1}$$
where $a$ and $b$ are the equatorial and polar radius of the Earth, respectively, and $a>b$. (The $\approx$ symbol will from now on mean equality up to higher-order terms in $f$.) We assume that Earth is a massive oblate ellipsoid
$$ \left(\frac{x}{a}\right)^2 + \left(\frac{y}{a}\right)^2 +\left(\frac{z}{b}\right)^2 ~\leq~ 1 \tag{2} $$
with uniform density $\rho$ and volume
$$V~=~\frac{4}{3}\pi a^2 b~\approx~\frac{4}{3}\pi b^3(1+2f).\tag{3}$$
The eccentricity is
$$e~=~\sqrt{(2-f)f}~\approx~ \sqrt{2f}~\ll ~1.\tag{4}$$
II) Quadrupole. We will assume that gravity is Newtonian. What we need to calculate is the quadrupole moment contribution to the gravitational potential
$$ U_4 ~= ~\frac{x^i Q_{ij}x^j}{r^5},\tag{5} $$
For symmetry reasons, one of the principal directions of the quadrupole moment $Q_{ij}=Q_{ji}$ must be along the polar $z$-axis, and the two other principal directions must have same eigenvalues and be in the equatorial $xy$ plane. Thus $U_4$ must be of the form
$$ U_4 ~=~ \frac{Q_a s^2+Q_b c^2}{r^3}, \tag{6}$$
where $Q_a$ and $Q_b$ are the equatorial and polar eigenvalues, respectively; where we have introduced the short-hand notation $s\equiv\sin(\theta)$ and $c\equiv\cos(\theta)$; and where $\theta\in[0,\pi]$ is the polar angle. Since the quadrupole moment $Q_{ij}$ cannot contribute to Gauss' law, we must demand that the Laplacian $\nabla^2 U_4 = 0$ vanishes, which leads to that the polar eigenvalue is minus two times the equatorial eigenvalue, $Q_b=-2Q_a$. In other words,
$$ U_4 ~\sim~ \frac{3c^2-1}{r^3}~=~\frac{2-3s^2}{r^3}.\tag{7} $$
It is therefore enough to calculate the gravitational potential at a point $P=(0,0,r)$ on the polar $z$-axis, where we have azimuthal $U(1)$ symmetry $\varphi\to\varphi+\varphi_{0}$.
III) Ring. Let us calculate the contribution to the gravitational potential energy from a ring parallel to the $xy$ plane and with center $(0,0,z)$ on the $z$-axis. Because of the $U(1)$ azimuthal symmetry, we may focus on a point on the ring with $y=0$ and $x\geq 0$, and which lie in the $xz$ plane. Let this point have 2D polar coordinates $(R,\theta)$. In other words, the point has 3D coordinates $(x,y,z)=(Rs,0,Rc)$. From the equation of an ellipse
$$ \left(\frac{x}{a}\right)^2 +\left(\frac{z}{b}\right)^2 ~=~ 1, \tag{8} $$
we obtain, after some elementary algebra,
$$ R ~\approx~ b (1 + fs^2) \quad \left\{\begin{array}{rcl}~=~ b &\mathrm{for}& \theta~=~0, \\ ~\approx~ a &\mathrm{for}& \theta~=~\frac{\pi}{2} . \end{array}\right.\tag{9}$$
Next, the distance $\ell$ from the ring to the point $P$ is given by the cosine relation
$$ \ell^2 ~=~ r^2 +R^2 - 2rRc ~\approx~ r^2 +b^2(1 + 2fs^2) - 2rb(1 + fs^2)c$$ $$~=~(r^2 +b^2 -2rbc) +2fbs^2(b-rc).\tag{10}$$
The "outer" surface area of the ring is
$$ dA ~=~ 2\pi Rs \cdot R d\theta ~\approx~ 2\pi b^2(1 + 2fs^2) (-dc), \tag{11}$$
with infinitesimal "radial" thickness
$$ dR ~\approx~ (1 + fs^2) db.\tag{12}$$
Thus the volume of the ring is
$$d^2V ~=~ dA \ dR ~\approx~ 2\pi b^2db(1 + 3fs^2) (-dc). \tag{13}$$
(Eventually we want to integrate over the polar angle $\theta$ from $0$ to $\pi$. This corresponds to integrate over $c\equiv\cos(\theta)$ from $1$ to $-1$ in the negative direction. Therefore $dc<0$ is negative.)
IV) Potential. The contribution from the ring to the gravitational potential energy $d^2 U$ at the point $P$ is
$$ d^2 U ~=~ -G \frac{d^2 M}{\ell} ~=~ - G\rho \frac{d^2 V}{\ell} ~\approx~ G\rho\frac{2\pi b^2db(1 + 3fs^2) dc}{\sqrt{(r^2 +b^2 -2rbc) +2fbs^2(b-rc)}} $$ $$~\approx~ G\rho2\pi b^2db\ dc\left(\frac{1 + 3f(1-c^2)}{\sqrt{r^2 +b^2 -2rbc}}- \frac{1}{2} \frac{2fb(1-c^2)(b-rc)}{(r^2 +b^2 -2rbc)^{\frac{3}{2}}}\right).\tag{14}$$
Integration over cosine $c$ of the polar angle $\theta\in[0,\pi]$ to form a thin shell yields (with the help of the MAPLE program)
$$ dU ~=~ \int_{c=1}^{c=-1} d^2 U ~\approx~ -G\rho\frac{4\pi b^2db}{r}\left( (1+3f) -f( 1+ \frac{2b^2}{5r^2}) - \frac{4fb^2}{15r^2}\right) $$ $$ ~=~ -G\rho\frac{4\pi b^2db}{r}\left( (1+2f) - \frac{2fb^2}{3r^2}\right).\tag{15} $$
Next we integrate over $b$ to get the polar potential of a massive ellipsoid
$$ U ~=~ \int^{b} dU ~\approx~ -G\rho\frac{4\pi b^3}{3r}\left( (1+2f) - \frac{2fb^2}{5r^2}\right) ~=~U_1+U_{4b}.\tag{16}$$
V) Monopole & quadrupole. The first term $U_1$ is, of course, precisely the monopole potential
$$U_1(r)~=~-\frac{GM}{r},\tag{17}$$
and the second term $U_{4b}$ is the polar quadrupole potential
$$U_{4b}~\approx~-\frac{2fb^2}{5r^2}U_1.\tag{18}$$
Thus we know that the gravitational quadrupole potential in an arbitrary point (not necessarily on the $z$-axis) is
$$U_{4}(r)~\approx~-\frac{(2-3s^2)fR^2}{5r^2}U_1(r),\tag{19}$$
and the full gravitational potential is
$$U(r) ~\approx~ U_1(r)+U_4(r) ~\approx~ -\frac{GM}{r}\left(1 - \frac{(2-3s^2)fR^2}{5r^2}\right). \tag{20}$$
VI) Surface. From now on, let us only consider points $P$ with
$$r~=~R~\approx~ b (1 + fs^2)\tag{21}$$
on the surface of the ellipsoid-shaped Earth. Then the monopole potential in $P$ becomes
$$ U_1(R)~ \approx~ -(1- fs^2)\frac{GM}{b}, \tag{22}$$
while the quadrupole potential in $P$ becomes
$$ U_4(R) ~\approx~ \frac{f(2-3s^2)}{5}\frac{GM}{b} \tag{23}$$
so that the full gravitational potential in $P$ becomes
$$ U(R) ~\approx~ U_1(R) + U_4(R) ~\approx~ -(1-\frac{2f(1+s^2)}{5})\frac{GM}{b}.\tag{24}$$
VII) Discussion. Let us use the North pole as a reference point, i.e., subtract the gravitational potential $U_b(b)$ at the North pole. Then the gravitational potential in $P$ becomes
$$ U(R)- U_b(b) ~\approx~ \frac{2fs^2}{5} \frac{GM}{b}~\approx~ \frac{2}{5}g hs^2,\tag{25} $$
where
$$h~:~=~a-b~=~fa~ \approx~ fb ~>~0\tag{26}$$
is the difference in the equatorial and polar radii, and
$$g~=~ |\vec{\nabla} U(R)| ~=~ \left. \frac{\partial U(r)}{\partial r}\right|_{r=R}+{\cal O}(f) ~=~ \frac{GM}{R^2}+{\cal O}(f).\tag{27}$$
We should now add the centrifugal potential
$$ V ~=~ -\frac{(\omega Rs)^2}{2}. \tag{28}$$
To the order in $f\ll 1$ that we are working, we see that the total potential $U(R)+V$ is constant (independent of the surface point $P$), if
$$\frac{2}{5}gh~=~\frac{(\omega R)^2}{2}. \tag{29}$$
Conclusion: We find a factor $\color{Red}{\frac{2}{5}}$ in difference from Mark Eichenlaub's original monopole argument.
Update: Half a year after this answer was posted on Phys.SE, on November 18th, 2011, the Wikipedia page changed its listed mathematical expression for the flatness parameter into $\frac{h}{R}\approx f\approx\frac{5}{4}\frac{\omega^2R}{g}$, and is now in full agreement with this answer.
There was some doubt about Lubos' answer (which I've accepted), so this is just a verification.
I copied the method Lubos described and found the potential difference for an ellipsoid with different eccentricities. Sure enough, for an oblate spheroid, if you make the center-equator distance a fraction $e$ larger than the center-pole distance, the potential is roughly a fraction $e/2$ smaller at the equator.
To do the whole problem, we'd have to take into account the varying density of the Earth, but as a rough estimate this seems to do the job.
For example, for the oblate spheroid
$$\left(\frac{x}{1.01}\right)^2 + \left(\frac{y}{1.01}\right)^2 + z^2 < 1$$
the average value of $1/r$ at $(0,0,1)$ is about .996, and the average at $(1.01,0,0)$ is about .991.
Python code below (please excuse the amateur-ness)
import random
import math
points = 10000000
e = .01
rad = 1+e
diam = 2*rad
pot= 0
count = 0
for i in range(1,points):
x = diam*random.random()-rad
y = diam*random.random()-rad
z = diam*random.random()-rad
r = math.sqrt((x-rad)*(x-rad)+y*y+z*z)
if x*x/(rad*rad)+y*y/(rad*rad) + z*z < 1:
pot = pot + 1/r
count = count + 1
print pot/count
pot2 = 0
count = 0
for j in range(1,points):
x = diam*random.random()-rad
y = diam*random.random()-rad
z = diam*random.random()-rad
r = math.sqrt(x*x+y*y+(z-1.0)*(z-1.0))
if x*x/(rad*rad)+y*y/(rad*rad) + z*z < 1:
pot2 = pot2 + 1/r
count = count + 1
print pot2/count