In this research, a modified form of Vening Meinesz-Moritz (VMM) theory of isostasy for the second-order radial derivative of gravitational potential, measured from the Gravity field and steady-state Ocean Circulation Explorer (GOCE), is developed for local Moho depth recovery. An integral equation is organised for inverting the GOCE data to compute a Moho model in combination with topographic/bathymetric heights of SRTM30, sediment and consolidated crystalline basement and the laterally-varying density contrast model of CRUST1.0. A Moho model from EGM2008 to degree and order 180 is also computed based on the same principle for the purpose of comparison. In addition, we compare both of them with the 3 available seismic Moho models; two global and one regional over the Indo-Pak region. Numerical results show that our GOCE-based Moho model is closer to the all seismic models than that of EGM2008. The model is closest to the regional one with a standard deviation of 5.5 km and a root mean squares error of 7.8 km, which is 2.3 km smaller than the corresponding one based on EGM2008.