A spline chaos expansion, referred to as SCE, is introduced for uncertainty quantification analysis. The expansion provides a means for representing an output random variable of interest with respect to multivariate orthonormal basis splines (B-splines) in input random variables. The multivariate B-splines are built from a whitening transformation to generate univariate orthonormal B-splines in each coordinate direction, followed by a tensor-product structure to produce the multivariate version. SCE, as it stems from compactly supported B-splines, tackles locally prominent responses more effectively than the polynomial chaos expansion (PCE). The approximation quality of the expansion is demonstrated in terms of the modulus of smoothness of the output function, leading to the mean-square convergence of SCE to the correct limit. Analytical formulae are proposed to calculate the mean and variance of an SCE approximation for a general output variable in terms of the requisite expansion coefficients. Numerical results indicate that a low-order SCE approximation with an adequate mesh is markedly more accurate than a high-order PCE approximation in estimating the output variances and probability distributions of oscillatory, nonsmooth, and nearly discontinuous functions.