In the context of stochastic two-phase flow in porous media, one is often interested in estimating the statistics of fluid saturations in the reservoir. In this work, we show how we can efficiently compute the probability distribution of water saturation in channelized porous systems by pursuing the methodology proposed inIbrahima et al. (Transp Porous Media 109:(1):81-107, 2015; Comput Geosci 22:(1):389-412, 2018). One of the key metrics of the developed distribution method, called the frozen streamline (FROST) method, is the logarithm of time-of-flight (log-TOF). Therefore, we dedicate asignificant portion of the paper to study the cumulative distribution function and probability density function of the log-TOF. We also demonstrate that a Gaussian mixture model can be applied for parametrization of the probability distribution of the log-TOF random field, which is expected to lead to significant computational speedup. We compare the results of the saturation statistics obtained from both the FROST method and exhaustive Monte Carlo (MC) simulations of coupled two-phase flow and transport (i.e., streamlines may change with time), and we show that the agreement between the two approaches is good and improves with grid refinement. Finally, we compare FROST with MC simulations that correspond to the exact assumptions used in FROST (i.e., the streamlines are time independent), and we demonstrate excellent agreement.