​The pollution of groundwater, essential for supporting populations and agriculture, can have catastrophic consequences. Thus, accurate modeling of water pollution at the surface and in groundwater aquifers is vital. Here, we consider a density-driven groundwater flow problem with uncertain porosity and permeability. Addressing this problem is relevant for geothermal reservoir simulations, natural saline-disposal basins, modeling of contaminant plumes and subsurface flow predictions. This strongly nonlinear time-dependent problem describes the convection of a two-phase flow, whereby a liquid flows and propagates into groundwater reservoirs under the force of gravity to form so-called “fingers”. To achieve an accurate numerical solution, fine spatial resolution with an unstructured mesh and, therefore, high computational resources are required. Here we run a parallelized simulation toolbox ug4 with a geometric multigrid solver on a parallel cluster, and the parallelization is carried out in physical and stochastic spaces. Additionally, we demonstrate how the ug4 toolbox can be run in a black-box fashion for testing different scenarios in the density-driven flow. As a benchmark, we solve the Elder-like problem in a 3D domain. For approximations in the stochastic space, we use the generalized polynomial chaos expansion. We compute the mean, variance, and exceedance probabilities for the mass fraction. We use the solution obtained from the quasi-Monte Carlo method as a reference solution.