Accurate modeling of contamination in subsurface flow and water aquifers is crucial for agriculture and environmental protection. Here, we demonstrate a parallel algorithm to quantify the propagation of uncertainty in the dispersal of pollution in subsurface flow. Specifically, we consider the density-driven flow and estimate how uncertainty from permeability and porosity propagates to the solution. We take a two-dimensional Elder-like problem as a numerical benchmark, and we use random fields to model our limited knowledge on the porosity and permeability. We use the well-known low-cost generalized polynomial chaos (gPC) expansion surrogate model, where the gPC coefficients are computed by projection on sparse tensor grids. The numerical solver for the deterministic problem is based on the multigrid method and is run in parallel. Computation of high-dimensional integrals over the parametric space is done in parallel too.