Simulation of multiphase flow in natural subsurface formations include selection of time-step size, i.e., the discrete snapshots (steps) over which the time-dependent nonlinear process is investigated. The simulation results would naturally depend on the size of the discrete time step, according to the order of accuracy of the implemented numerical scheme. For accurate analyses the required time steps might be needed to be very small, which then leads to very costly simulations, which are often times impractical in real-world simulation applications. Taking coarse-scale time-step sizes can reduce the computational costs, however, it results in lower simulation accuracy. Note that a common simulation practice is to take constant time-step everywhere in the computational domain. To resolve this challenge, this MSc thesis report presents a novel local time-stepping (LTS) method for multiphase flow in porous media. In this method, different time-steps can be applied in one computational domain; allowing for employing small time steps in sub-domains with higher sensitivity to the time step size and bigger time steps elsewhere. The LTS method is developed for all types of sequentially-coupled and fully-coupled simulation strategies. For all of these developments, the interfacial flux continuity is the key requirement to connect the sub-domains of different time zones. The fluxes chosen for the inner stages of the coarse time-step are consistent with the chosen simulation strategy. More precisely, in the case of the implicit-pressure-explicit-saturation (IMPES) method these fluxes are lagged in time (explicit), while for the sequential implicit and fully implicit (FIM) methods they are calculated at the current (implicit) time-step. For implicit transport method, the refined zone partition matrices are developed to algebraically form the refined zone Jacobian matrix and refined zone residual vector. The numerical results show that only a fraction of the domain can be simulated with small time steps; while the rest can be simulated at a bigger one. Local time-stepping method provides accurate simulation results compared with the fine-scale time-step simulation results for all three simulation strategies. Specially, the local time-stepping for FIM simulations, which is the mostly common commercial-grade simulation approach, can preserve the accuracy of the results and also the computational efficiency. For the studied cases, the computational efficiency gain of LTS for sequential implicit method is not as significant as that of the FIM. To exploit full adaptivity in simulation, as presented in Appendix B, one can directly apply the developed method of this thesis to adaptive implicit method (AIM); where sequential and fully implicit coupling approaches are applied in different sub-domains of the reservoir. Overall, LTS casts a promising approach to optimize accuracy-efficiency tradeoff when it comes to time step selection.