We formulate an error propagation model based on solving the Vening Meinesz-Moritz (VMM) inverse problem of isostasy. The system ofobservation equations in the VMM model defines the relation between theisostatic gravity data and the Moho depth  by means of a second-order Fredholm integralequation of the first kind. The corresponding error model (derived in aspectral domain) functionally relates the Moho depth errors with the commissionerrors of used gravity and topographic/bathymetric models. The error model alsoincorporates the non-isostatic bias which describesthe disagreement, mainly of systematic nature, between the isostatic andseismic models. The error analysis is conducted at the study area of theTibetan Plateau and Himalayas with the world largest crustal thickness. TheMoho depth uncertainties due to errors of the currently available globalgravity and topographic models are estimated to be typically up to 1-2 km,provided that the GOCE gravity gradient observables improved themedium-wavelength gravity spectra. The errors due to disregarding sedimentarybasins can locally exceed ~2 km. The largest errors (which cause a systematic bias betweenisostatic and seismic models) are attributed to unmodeled mantleheterogeneities (including thecore-mantle boundary) and other geophysical processes. These errors aremostly less than 2 km under significant orogens (Himalayas, Ural), but canreach up to ~10 km under the oceanic crust.