We derive analytical expressions for the one-point cumulative distribution function and probability density function of water saturation in the one-dimensional stochastic immiscible two-phase (Buckley-Leverett) problem. The sources of uncertainty are the spatial distributions of porosity and total velocity. The derived distribution functions involve integrals of the input random parameters. Comparisons with standard Monte Carlo simulation demonstrate that the method is applicable for input parameters with large variance and arbitrary correlation lengths. We also show that the proposed method is superior to the low-order statistical moment equations approach. We also outline a streamline-based strategy to extend our distribution-based method to multiple spatial dimensions.