In this study, we develop a probabilistic approach to map the parametric uncertainty to the output state uncertainty in first-order hyperbolic conservation laws. We analyze this problem for nonlinear immiscible two-phase transport in heterogeneous porous media in the presence of a stochastic velocity field. The uncertainty in the velocity field can arise from the incomplete description of either porosity field, injection flux, or both. The uncertainty in the total-velocity field leads to the spatiotemporal uncertainty in the saturation field. Given information about the spatial/temporal statistics of the correlated heterogeneity, we leverage method of distributions to derive deterministic equations that govern the evolution of single-point CDF of saturation. Unlike Buckley Leverett equation, the equation for the raw CDF function is linear in space and time. Hereby, we give routes to circumventing the computational cost of Monte Carlo scheme while obtaining the full statistical description of saturation. We conduct a set of numerical experiments and compare statistics of saturation computed with the method of distributions, against those obtained using the statistical moment equations approach and kernel density estimation post-processing of high-resolution Monte Carlo simulations. This comparison demonstrates that the CDF equations remain accurate over a wide range of statistical properties, i.e. standard deviation and correlation length of the underlying random fields, while the corresponding low-order statistical moment equations significantly deviate from Monte Carlo results, unless for very small values of standard deviation and correlation length.