Uncertainty quantification in complex physical models is often challenged by the computational expense of these models. One often needs to operate under the assumption of sparsely available model simulations. This issue is even more critical when models include a large number of input parameters. This “curse of dimensionality”, in particular, leads to a prohibitively large number of basis terms in spectral methods for uncertainty quantification, such as polynomial chaos (PC) methods. In this work, we implement a PC-based surrogate model construction that “learns” and retains only the most relevant basis terms of the PC expansion, using sparse Bayesian learning. This dramatically reduces the dimensionality of the problem, making it more amenable to further analysis such as sensitivity or calibration studies. The model of interest is the community land model with about 80 input parameters, which also exhibits nonsmooth input-output behavior. We enhanced the methodology by a clustering and classifying procedure that leads to a piecewise-PC surrogate thereby dealing with nonlinearity. We then obtain global sensitivity information for five outputs with respect to all input parameters using less than 10,000 model simulations−a very small number for an 80-dimensional input parameter space.