The inefficacy and limitations of commonly used geostatistical estimation and simulation methods based on two-point variogram models in resolving complex and nonlinear geological settings have been thoroughly investigated in recent years by both mineral and petroleum geoscientists. Due to lack of observing multiple real data simultaneously in the course of simulation, the geostatistical simulations made by two-point variograms usually produce models that are highly different and scattered, which are not representative of actual geological structures and continuity. To overcome such deficiencies, a number of approaches that incorporate the actual geological patterns beyond two-point geostatistics have been introduced. Multiple point geostatistics is one of the best methods capable of simulating nonlinear and complex structures. To alleviate the shortage of available data in real situations, the multiple point geostatistics uses training images, which are built upon our knowledge from previous geological studies and entire available exploratory data. In current study, the problem of post mineralization intrusions into the main Sarcheshmeh porphyry copper deposit is addressed. Such geological complexities may cause deviations in ore reserve estimation using conventional methods. The non-mineralized dikes, which have been injected in the Sarcheshmeh copper deposit in a number of post mineralization phases, are considered as barren occurrences. Inclusion of such barren dikes into ore reserve estimation process, will cause over estimates of ore reserve due to assigning metal content in place of barren dikes. In addition, the complex and almost irregular injection of such dikes into main ore body may result in producing nonlinear and complex structures, which cannot be modeled accurately with conventional two-point geostatistics methods. Therefore, the main objective of current study was set to evaluate the possibility of employing multiple-point geostatistics approach in ore reserve simulation and estimation followed by comparing the results with other conventional variogram based methods. To achieve such objective, it was necessary to prepare suitable training image, capable of reproducing 3D spatial distribution of dikes in estimation area. Then such training image was inputted into the multiple point simulation algorithm known as SNESIM. In order to get results that are more accurate in mapping 3D spatial distribution of dikes, the main parameters of SNESIM algorithm were to be optimized. In the next step, the optimized SNESIM algorithm along with training image was used to produce 50 realizations for dike occurrences in the entire area of reserve estimation space. Next, these 50 realizations were used to map the dike occurrences using the so-called E-type estimator. The results were satisfactory showing consistent spatial continuity and distribution of dikes in plotted plans and sections. The ore grade and reserve were estimated using ordinary and indicator kriging without including dikes in deposit blocky model. The results showed the superiority of indicator kriging in detecting dikes compared with conventional ordinary kriging. In the next step, the dikes model was embedded into the estimated grades for both methods resulting in the more accurate reserve estimations. Finally, the results were validated using the data of boreholes, which were set aside for cross validation. It was concluded that the combination of multiple point geostatistics with indicator kriging method, will result in more accurate and realistic ore reserve estimation for porphyry copper deposits that are highly heterogeneous and intruded by post mineralization barren dikes.