Start with the Schwarzschild solution in the weak field approximation:
(cd\tau)^2=(1-\frac{2\Phi}{c^2})(cdt)^2+(1-\frac{2\Phi}{c^2})^{-1}(dr)^2+...
That should be:
{\color{red}(cd\tau)^2=(1+\frac{2\Phi}{c^2})(cdt)^2-(1+\frac{2\Phi}{c^2})^{-1}(dr)^2-...}
See post 185 of this different thread
https://www.physicsforums.com/showthread.php?t=397403&page=12 where you gave the correct equation. This error in the signs propagates all the way through the rest of your calculations.
For the case dr=d\theta=d\phi=0 you get the well known relationship:
d\tau=\sqrt{1{\color{red}+}\frac{2\Phi}{c^2}}dt
Writing the above for two different gravitational potentials \Phi_1 and \Phi_2 you obtain the well-known time dilation relationship:
\frac{d\tau_1}{d\tau_2}=\sqrt{\frac{1{\color{red}+}\frac{2\Phi_1}{c^2}}{1{\color{red}+}\frac{2\Phi_2}{c^2}}}
At the Earth surface :
\Phi_1=-\frac{GM}{R}
At the Earth center:
\Phi_2=-3/2\frac{GM}{R}
Now, due to the fact that \frac{\Phi}{c^2}<<1 you can obtain the approximation:
\frac{d\tau_1}{d\tau_2}=1{\color{red}+}\frac{\Phi_1-\Phi_2}{c^2}=1{\color{red}+}\frac{GM}{2Rc^2}{\color{red}>}1
So, f_1>f_2 where f_1 is the clock frequency on the Earth crust and f_2 is the frequency of the clock at the center of the Earth.
When corrected, your aproximation is in close agreement with the equations I gave and is on the right side of unity.
...