In 2013 the Grozny State Oil Technical University together with members of the "Geothermal resources" consortium started a pilot project to build a geothermal plant on the basis of the most promising XIII layer of the Khankala thermal waters deposit. Planned capacity of the projected facility is 22.8 GJ/h, which will heat a greenhouse complex. Geostatistical analysis and numerical simulations are required for achieving sustainability in exploitation of the resource and in order to define new drilling sites. In this paper intrinsic random functions of k order kriging (IRF-k) is used for the XIII layer map creation and estimation of the temperature distribution within it. The results obtained by geostatistical approach are used for numerical simulation of the hydro-thermal process relied to reinjection of cold thermal water back into the productive layer. The work carried out allows highlighting the most promising zones for exploitation and shows that the geothermal water temperature in the production well starts to decrease after 6 years. (C) 2015 Elsevier Ltd. All rights reserved.