Recent developments with polynomial chaos expansions with random coefficients facilitate the accounting for subscale features, not captured in standard probabilistic models. These representations provide a geometric characterization of random variables and processes, which is quite distinct from the characterizations (in terms of probability density functions) typically adapted to Bayesian analysis. Given the importance of Bayes theorem within probability theory, it is important to synthesize the connection between these two representations. In this talk, we will describe a hierarchical Bayesian framework that introduces polynomial chaos expansions with random parameters as a consequence of Bayesian data assimilation. We will provide insight into the behavior and use of these expansions and exemplify them through a multiscale application from thermal science. Specifically, information collected from fine scale simulations is used to construct stochastic reduced order models. These coarse models are indexed in terms of specimen-to-specimen variability and also in terms of variability in their subscale features. The ability of these doubly-stochastic expansions to improve the predictive value of model-based simulations is highlighted.