Abstract The popularity of parallel reservoir simulation has increased in recent years due to availability of multi-core CPUs and the opportunity for run-time reduction, especially in complex and fractured reservoir models that include large amounts of data. However, the accuracy of classical domain decomposition techniques significantly decreases as the number of processing units increases. The scope of this work is mainly concentrated on development of multi-threaded algorithms to overcome the mentioned problem in parallel reservoir simulations. In this work, at first, a newly developed, parallel, three dimensional, fully implicit, three-phase black-oil reservoir simulator is introduced. For solving large scale linear systems produced in giant black-oil models, which are probable cases in real world reservoir models, BiCG-Stabilized linear solver, enhanced with various preconditioners including, Blocked Incomplete LU factorization, Constrained Pressure Residual using Algebraic Multi-Grid as its pressure smoother is constructed with a variety of parallelization techniques applied in each part. The presented multi-threaded algorithms can utilize any user-defined number of CPU cores. The algorithm has been evaluated by SPE10 with 1122000 grid blocks and a real data file containing 277875 grid blocks. The result indicates that increasing the number of CPU threads does not have any negative effect on the simulation performance in case of high memory bandwidth limits. For SPE10 data file, 1.70, 2.4, and 2.3 speedup ratios (in comparison with a single thread) has been observed for 2, 4, and 6 parallel threads, respectively. In the real case, the simulation speed is increased by a factor of 1.49, 1.74, and 1.75 for 2, 4, and 6 parallel threads, respectively. Using multi-threaded algorithms as a novel approach in parallel simulations, the solution configuration does not change from the single core case and follows the same exact path. This however is not the case for other simulators that use domain-decomposition techniques. These techniques insert numerical errors in the solution of linear systems, which lead to convergence problems in the whole system.