The global gravity and crustal models are used in this study to determine the regional Moho model. For this purpose, we solve the Vening Meinesz-Moritz's (VMM) inverse problem of isostasy defined in terms of the isostatic gravity gradient. The functional relation between the Moho depth and the second-order radial derivative of the VMM isostatic potential is formulated by means of the (linearized) Fredholm integral equation of the first kind. Methods for a spherical harmonic analysis and synthesis of the gravity field and crustal structure models are applied to evaluate the gravity gradient corrections and the respective corrected gravity gradient, taking into consideration major known density structures within the Earth's crust (while mantle heterogeneities are disregarded). The resulting gravity gradient is compensated isostatically based on applying the VMM scheme. The VMM inverse problem for finding the Moho depths is solved iteratively. The regularization is applied to stabilize the ill-posed solution. The global geopotential model GOCO-03s, the global topographic/bathymetric model DTM2006.0 and the global crustal model CRUST1.0 are used to generate the VMM isostatic gravity gradient with a spectral resolution complete to a spherical harmonic degree of 250. The VMM inverse scheme is used to determine the regional isostatic crustal thickness beneath the Tibetan Plateau and Himalayas (compiled on a 1x1 arc-deg grid). The differences between the isostatic and seismic Moho models are modeled and subsequently corrected for by applying the non-isostatic correction. Our results show that the regional gravity gradient inversion can model realistically the relative Moho geometry, while the solution contains a systematic bias. We explain this bias by more localized information on the Earth's inner structure in the gravity gradient field compared to the potential or gravity fields.