Why is the Earth so fat? I made a naive calculation of the height of Earth's equatorial bulge and found that it should be about 10km.  The true height is about 20km.  My question is: why is there this discrepancy?
The calculation I did was to imagine placing a ball down on the spinning Earth.  Wherever I place it, it shouldn't move.
The gravitational potential per unit mass of the ball is $hg$, with $h$ the height above the pole-to-center distance of the Earth (call that $R$) and $g$ gravitational acceleration.
Gravity wants to pull the ball towards the poles, away from the bulge.  It is balanced by the centrifugal force, which has a potential $-\omega^2R^2\sin^2\theta/2$ per unit mass, with $\omega$ Earth's angular velocity and $\theta$ the angle from the north pole.  This comes taking what in an inertial frame would be the ball's kinetic energy and making it a potential in the accelerating frame.
If the ball doesn't move, this potential must be constant, so
$$U = hg - \frac{(\omega R \sin\theta)^2}{2} = \textrm{const}$$
we might as well let the constant be zero and write
$$h = \frac{(\omega R \sin\theta)^2}{2g}$$
For the Earth, 
$$R = 6.4*10^6m$$
$$\omega = \frac{2\pi}{24\ \textrm{hours}}$$
$$g = 9.8\ m/s^2$$
This gives 10.8 km when $\theta = \pi/2$, so the equatorial bulge should be roughly that large.  
According to Wikipedia, the Earth is 42.72 km wider in diameter at the equator than pole-to-pole, meaning the bulge is about twice as large as I expected.  (Wikipedia cites diameter; I estimated radius.)
Where is the extra bulge coming from?  My simple calculation uses $g$ and $R$ as constants, but neither varies more than a percent or so.  It's true the Earth does not have uniform density, but it's not clear to me how this should affect the calculation, so long as the density distribution is still spherically-symmetric (or nearly so).
(Wikipedia also includes an expression, without derivation, that agrees with mine.)
 A: 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

A: In this answer, I will present a framework to use, and then I will frame the prior answers within that framework.  Let me sum up the values we have here.  I'll use the same notation (as best as possible) as everyone else and Wikipedia for an oblate spheroid where $a$ is the large, equatorial, radius.


*

*Mark1, method in the question, $2 (a-b)  = 21.6 km$

*Mark2, method in the answer, past answer times 5/2 for $2 (a-b) = 54.0 km$


Here is my approach:
The mass of the earth can be taken a combination of two shapes, an inner sphere with radius $b$ and an edge volume, which is the oblate spheroid minus the inner sphere.  The Earth has a certain average density defined by $\rho = M/V$, but this can be divided up into two different type of materials, the core, and the crust.  The total mass requirement will then require that $M = V_{core} \rho_{core} + V_{edge} \rho_{edge}$, while the average density of the Earth requires $V = V_{core} + V_{edge}$, constraining $a$ and $b$ by one degree of freedom.  When writing a code we can say that $a$ always implies $b$ according to $b=3V/(4 \pi a^2)$, the crust density also can be taken to imply the core density.  Then the equipotential surface constraint dockets another degree of freedom, which can be used to iteratively find the value of $a$.  Illustration:

The potential from the inner sphere is easy.  I'll write it for a point on the equator and on the pole, combined with the other potentials.
$$U_a = G \frac{M_{core}}{a} + U_{a,edge} + \frac{1}{2} \omega^2 a^2$$
$$U_b = G \frac{M_{core}}{b} + U_{b,edge}$$
Obviously, the hard part is calculating the potential from the ridiculously irregular edge shape.  Before that, however, it's important to think of the physical implications of looking at the problem this way.  To begin with, what are relevant densities relevant to Earth?  Here is the average density, and then the density within about 0 to 200 km of the surface.
$$rho_{Earth} = 5520. \frac{kg}{m^3}$$
$$rho_{crust} = 3400. \frac{kg}{m^3}$$
When we actually solve the problem, we will specify the crust density, and then that will imply a density of the core sphere.  Is this method accurate?  No.  The main thing it misses is that the center sphere is not a spherically symmetric distribution of matter.  The densities will theoretically stratify according to constant potential lines.  In other words, if there were a high-density core of the Earth, it would also be an oblate spheroid.  Because of that reason, introducing the inner sphere does miss some detail, but this model could still be pretty good.
Implementing this is a little tricky, as others have pointed out, but focusing the calculations on the edge volume helps a great deal.  I also used a Monte Carlo method, but subtracted out the core volume.  That is to say, I shot points randomly into the edge volume as efficiently as reasonable.  In order to do this I used a weighted method, and it seems to have worked out okay.  With 5 million iterations I found random variations in the calculation of the potential at $a$ to have a standard deviation of about $15 m$ of gravitational potential equivalent and around $2 m$ of gravitational potential equivalent at $b$.  The reason for the higher deviation at $a$ is because it has more mass close to it, and the sampling was unbiased in the yz-plane for calculations of the $a$ potential and unbiased in the xy-plane for calculations of the $b$ potential.  Anyway, as I iterated on the $a$ value, I set a tolerance for $100 m$, because this should be considerably higher than the statistical variation and it's good enough for calculation of the bulge.  To recap, this is my method:
For calculation of potential at $a$


*

*Sample two random values for the $y$ and $z$ values between $0$ and $a$

*If these two values lie outside of the ellipse of $\frac{y^2}{a^2}+\frac{z^2}{b^2}<1$ then sample two new values and try again (gives something like a 30% loss of efficiency)

*Sample an x-value between the surface of the inner sphere and the outer spheroid.  If the (y,z) pair doesn't fall within the inner sphere, sample between $0$ and the outer spheroid surface.

*Tally the potential between the sampled (x,y,z) point and (a,0,0), using the M_edge mass

*Repeat between (-x,y,z) and (a,0,0)

*Tally the weight for this sample as the distance between the two surfaces times two.

*The potential at (a,0,0) is then the total tallied potential divided by the total tallied weight.

*Repeat a similar method to find the potential at (0,0,b)

*Numerically root find to satisfy the equipotential condition discussed above.


I did this, and for the different values of the crust density, I got the following.


*

*$\rho_{crust} = 0$, implying gravitational field is insensitive to flatness, gets $2(a-b)= 21.8 \pm 0.1 km$

*$\rho_{crust} = 3400. \frac{kg}{m^3}$, a reasonable value for the crust density, gets $2(a-b) = 34.7 \pm 0.1 km$

*$\rho_{crust} = 5520. \frac{kg}{m^3}$, a fully homogenous density of the Earth gets $2(a-b) = 54.0 \pm 0.1 km$


I thought these are good results because the first and last of them come close to the prior answers within numerical error and the reasonable value for crust density gets closer to the real value of $42 km$.
If anyone is interested, I can look into putting the code for this on github or something similar.  Otherwise, it's a little longer than the others posted here, so for now I'll avoid cluttering this space.
A: Here I would like to numerically check the theoretical prediction of a factor $\frac{2}{5}$ in difference from Mark Eichenlaub's original monopole argument. In practice, this means calculating the difference in the gravitational potential between the North pole and the Equator, and divide by the corresponding difference in monopole potentials. For numerical reasons, it is better in practice to calculate the inverse (=reciproc) fraction, which then should be compared with $\frac{5}{2}$. Since my programming skills are limited, I have just written a slow MAPLE code to do the job.
 b:=100; f:=.10; a:=b*(1+f); V1:=evalf(4*Pi*a^2*b/3); 
 xa:=a; ya:= 0; za:=0; xb:=0; yb:=0; zb:=b; 
 U1a := evalf(V1/sqrt(xa^2 + ya^2 + za^2)); 
 U1b := evalf(V1/sqrt(xb^2 + yb^2 + zb^2));

 Ua:=0;Ub:=0;V:=0;
 for x from -a-.5 by 1 to a+.5 do 
 for y from -a-.5 by 1 to a+.5 do 
 for z from -b-.5 by 1 to b+.5 do 
 if (x/a)^2 + (y/a)^2 + (z/b)^2 < 1 then 
 Ua:=Ua + 1/sqrt((x-xa)^2 + (y-ya)^2 + (z-za)^2); 
 Ub:=Ub + 1/sqrt((x-xb)^2 + (y-yb)^2 + (z-zb)^2); 
 V:=V+1;
 end if;od;od;od;

 b;f;Ua;U1a;Ub;U1b;V;V1;Ub-Ua;U1b-U1a;(U1b-U1a)/(Ub-Ua);

The result with polar radius $b:=100$ was
$$
\begin{array}{cc} 
\mathrm{Flatness}&\mathrm{Inverse \ fraction}  \\
f & \frac{U_{1b}-U_{1a}}{U_b-U_a} \\ \\
0.10 & 2.5549 \\
0.05 & 2.5353 \\
0.03 & 2.5149 \\
0.02 & 2.5085 \\
0.01 & 2.5088 \\
\end{array}
$$
The fact that the estimate does not improve from flatness $0.02$ to $0.01$ is a lattice artifact, because the lattice spacing is of same order as the difference in radius between the North pole and the Equator. 
A: 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

