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. 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? 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). 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.
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 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 an ellipsoid with "radii" $(r_x,r_y,r_z)=(0.9,0.9,1.0)$. 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$ 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:**
<blockquote>
{xt, yt, zt} = {1.1, 0, 0};<br><br>
runs = 200000;<br>
totalRinverse = 0;<br>
total = 0;<br>
For[i = 1, i <= runs, i++,<br>
x = RandomReal[]*2 - 1;<br>
y = RandomReal[]*2 - 1;<br>
z = RandomReal[]*2 - 1;<br>
inside = x^2/0.81 + y^2/0.81 + z^2 < 1;<br>
total = If[inside, total + 1, total];<br>
totalRinverse =
totalRinverse +
If[inside, 1/Sqrt[(x - xt)^2 + (y - yt)^2 + (z - zt)^2], 0];<br>
]<br><br>
res1 = N[total/runs / (4 Pi/3/8)]<br>
res2 = N[totalRinverse/runs / (4 Pi/3/8)]<br>
res2/res1<br>
</blockquote>
**Description**
Use the Mathematica code above. The final number that is printed by the code is the average $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.