Previous studies have proposed the first-order reliability method (FORM) as an approach to quantitative stochastic analysis of subsurface transport. Most of these considered only simple analytical models of transport in homogeneous media. Studies that looked at more-complex, heterogeneous systems found FORM to be computationally demanding and were inconclusive as to the accuracy of the method. Here we show that FORM is poorly suited for computing point concentration cumulative distribution functions (cdfs) except in the case of a constant or monotonically increasing solute source. FORM is better equipped to predict transport in terms of the cumulative mass flux across a control surface. As a demonstration, we use FORM to estimate the cumulative mass flux cdf in two- dimensional, random porous media. Adjoint sensitivity theory is employed to minimize the computational burden. In addition, properties of the conductivity covariance and distribution are exploited to improve efficiency. FORM required eight times less CPU time than Monte Carlo simulation to generate the results presented. The accuracy of FORM is found to be minimally affected by the size of the initial solute body and the solute travel distance. However, the accuracy is significantly influenced by the degree of heterogeneity, providing an accurate estimate of the cdf when there is mild heterogeneity (σInK = 0.5) but a less accurate estimate when there is stronger heterogeneity (σInK = 1.0).