A probabilistic approach to model macroscopic behaviour of non-wetting-phase ganglia or blobs in multi-phase flow through porous media is proposed. The key idea is to consider a set of stochastic Markov processes that can mimic the microscopic multi-phase dynamics. These processes are characterized by equilibrium probability density functions (PDFs) and correlation times, which can be obtained from microscale simulation studies or experiments. A Lagrangian viewpoint is adopted, where stochastic particles represent infinitesimal fluid elements and evolve in the physical and probability space. Ganglion mobilization and trapping are modelled by a twostate jump process with transition probabilities given as functions of ganglion size. Coalescence and breakup of ganglia influence the ganglion size distribution, which is modelled by a Langevin type equation. The joint probability density function (JPDF) of the chosen stochastic variables is governed by a high-dimensional Chapman-Kolmogorov equation. This equation can be used to derive moment (e.g. saturation, mean mobility etc.) transport equations, which in general do not form a closed system. However, in some special cases, which arise in the limit of one time scale being smaller or larger than the others, a closed set of moment transport equations can be obtained. For slowly varying and quasi-uniform flows, the saturation transport equation appears in closed form with the mean mobility fully determined, if the equilibrium PDFs are known. Furthermore, it is shown how statistical parameters such as mobilization and trapping rates and equilibrium PDFs can be obtained from the birth-death type approach, in which ganglia breakup and coalescence are explicitly considered. A two-equation transport model (one equation for the total saturation and one for the trapped saturation) is obtained in the limit of very fast coalescence and breakup processes. This model is employed to mimic hysteresis in relative permeability-saturation curves; a well known phenomenon observed in the successive processes of imbibition and drainage. For the general case, the JPDFequation is solved using the stochastic particle method, which was proposed in our previous paper (Tyagi et al. J. Comput. Phys. 227, 2008, 6696-6714). Several one-and two-dimensional numerical simulation results are presented to show the influence of correlation times on the averaged macroscopic flow behaviour.