Coarse-scale models are generally preferred in the numerical simulation of multi-phase flow due to computational constraints. However, capturing the effects of fine-scale heterogeneity on flow and isolating the impacts of numerical (artificial) dispersion, which increases with scale, are not trivial. In this paper, a particle-tracking method is devised and integrated in a scale-up workflow to estimate the conditional probability distributions of multi-phase flow functions, which can be considered as inputs in coarse-scale simulations with existing commercial packages. First, a novel particle-tracking method is developed to solve the saturation transport equation. The transport calculation is coupled with a velocity update, following the implicit pressure, explicit saturation framework, to solve the governing equations of two-phase immiscible flow. Each phase particle is advanced in a deterministic convection step according to the phase velocity, as well as in a stochastic dispersion step based on the random Brownian motion. A kernel-based formulation is proposed for computation of fluid saturation in accordance with the phase particle distribution. A novel aspect is that this method employs the kernel approach to construct saturation from phase particle distribution, which is an important improvement to the conventional box method that necessitates a large number of particles per grid cell for consistent saturation interpolation. The model is validated against various analytical solutions. Finally, the validated model is integrated in a statistical scale-up procedure to calibrate effective, or "pseudo," multi-phase flow functions (e.g., relative permeability functions). The proposed scale-up framework does not impose any length scale requirement regarding the distribution of sub-grid heterogeneities.