We have developed a new Monte Carlo sampling method for simulating flow and transport in unsaturated porous media, characterized by van Genuchten–Mualem constitutive relations. Instead of sampling each individual soil parameter from its probability space and then running Monte Carlo simulations using realizations of rock parameters directly, we calculate retention curves from realizations of rock parameters, take subsamples from these retention curves, and run simulations using parameter realizations corresponding to these selected retention curves. The retention curve subsampling methodology was applied to three-dimensional simulations of conservative tracer transport beneath Material Disposal Area G at the Los Alamos National Laboratory. Convergence of the proposed sampling method was assessed by comparing statistics of breakthrough curves observed at a compliance boundary with those obtained using between 25 and 1000 Latin hypercube sampling (LHS) Monte Carlo simulations. Our example shows that 25 model runs based on selected retention curves could adequately approximate the results from our assumed truth (1000 LHS Monte Carlo simulations), while LHS alone required in excess of 50 realizations to achieve the same quality result. Another finding from this work is that the median of the breakthrough curves was more meaningful than the arithmetic mean of curves, and the former was nearly identical to the breakthrough curve derived from mean rock properties.