The increasing need of human for natural and mineral resources has attracted significant attentions towards developing new and more efficient simulation algorithms in order to enable mapping complicated geological settings. Among algorithms commonly used for simulating complex and non-linear variability, the multiple point geostatistical methods are gradually gaining more credit and applications. The aim of multiple point geostatistics is producing more reliable and geologically consistent model of natural phenomena through employing training images. In fact, training images represent spatial patterns and structures that would be expected to exist throughout ore deposit. The main objective of current study is to produce simulated realizations preserving large variability structures of third order hidden in the measured exploratory data. Since the porphyry copper deposits are inherently controlled by complex and non-linear phenomena, the central part of Sarcheshmeh copper deposit encompassing a number of 135 exploratory bore holes was considered as a good example of applying the proposed new simulation methodology. The well-known Filtersim algorithm was selected as the core engine for simulating such non-linear structures to be applied on four purposely-selected plans. In order to define training images indicative of large scale variations throughout 2D plans, a trend surface analysis of up to 3 rd degree was carried out for each of the plans separately. A second stationary training image was then produced through repeating the main 3 rd order trend in all directions.In order to maintain the governing large scale variations and in the meantime reducing the local high wavenumber noises in the course of simulation, a wavelet transform was applied on both of the training image and gridded data constructed through nearest neighbor interpolation. In this way, the simulated realizations created in wavelet domain via Filtersim algorithm are enhanced. In order to keep the image dimensions and inertia unchanged, discrete stationary wavelet transform was adopted in this study. The entire data was divided into training and testing data sets. Then, a wavelet transform was applied on both of the training images and gridded hard data up to one level and the conditional and unconditional simulations were performed for 100 realizations. Next, an e-type estimation were evaluated on simulated realizations after applying inverse discrete wavelet transform was made. Finally, through using a number of error measures such as MSE, RMSE and relative error the performance of the conditional and unconditional simulations were evaluated. The average relative error of computed for all plans was found to be in the order of 10 percent which is considered acceptable in the final stage of detailed exploration program. It was also concluded that the application of stationary discrete wavelet transform on the stationary training image has considerably reduced the undesirable edge effect created by the performance of filters in the wavenumber domain.