Many of the complex physical processes relevant for compositional multi-phase flow in porous media are well understood at the pore-scale level. In order to study CO2 storage in sub-surface formations, however, it is not feasible to perform simulations at these small scales directly and effective models for multi-phase flow description at Darcy scale are needed. Unfortunately, in many cases it is not clear how the micro-scale knowledge can rigorously be translated into consistent macroscopic equations. Here, we present a new methodology, which provides a link between Lagrangian statistics of phase particle evolution and Darcy scale dynamics. Unlike in finite-volume methods, the evolution of Lagrangian particles representing small fluid phase volumes is modeled. Each particle has a state vector consisting of its position, velocity, fluid phase information and possibly other properties like phase composition. While the particles are transported through the computational domain according to their individual velocities, the properties are modeled via stochastic processes honoring specified Lagrangian statistics. Note that the conditional expectations of the particle velocities are different for different fluid phases. The goal of this paper is to present the general framework for this alternative modeling approach. Various one and two-dimensional numerical experiments demonstrate that with appropriate stochastic rules the particle solutions are consistent with a standard two-phase Darcy flow formulation. In the end, we demonstrate how to model non-equilibrium phenomena within the stochastic particle framework, which will be the main focus of the future work.