Strongly degenerate parabolic convection–diffusion equations arise as governing equations in a number of applications such as traffic flow with driver reaction and anticipation distance and sedimentation of solid–liquid suspensions in mineral processing and wastewater treatment. In these applications several parameters that define the convective flux function and the degenerating diffusion coefficient are subject to stochastic variability. A method to evaluate the variability of the solution of the governing partial differential equation in response to that of the parameters is presented. To this end, a generalized polynomial chaos (gPC) expansion of the solution is approximated by its projection onto a finite-dimensional space of piecewise polynomial functions defined on a suitable discretization of the stochastic domain, according to the basic principle of the hybrid stochastic Galerkin (HSG) approach. This approach is combined with a finite volume (FV) method, resulting in a so-called FV–HSG method, to compute the sought deterministic coefficient functions of the truncated polynomial-chaos-based expansion of the solution. Since the stochastic parameter space is now spanned by piecewise polynomial functions, one may employ the numerical result to compute the reconstruction of the numerical solution for arbitrary values of the random variables. The expectation, the variance or other stochastic quantities of the solution (as functions of time and position) can also be computed from these coefficient functions. The method is illustrated by a number of numerical examples.