Constructing polynomial chaos expansions in the context of uncertainty quantification requires the computation of inner products that are typically high dimensional integrals involving an expensive to evaluate model and a set of orthogonal polynomials. This thesis explores the approach of first interpolating the model and then either performing pseudospectral projection (PSP) on this surrogate or, by exploiting the simpler structure of the interpolant, deriving analytical formulas for the inner products. To keep the number of model evaluations low, several variants of the sparse grid combination technique are proposed for interpolation, in particular the single dimension spatially adaptive refinement strategy. Comparisons for several test functions and the hydrological model HBV, with the number of model evaluations as reference, show that this approach often produces better results than direct PSP.