We present an implementation of a Bayesian mixture model using Hamiltonian Monte Carlo techniques to search for the spatial separation of Galactic dust populations. Utilizing intensity measurements from the Planck High Frequency Instrument, we apply this model to high-latitude Galactic dust emission. Our analysis reveals a strong preference for a spatially varying two-population dust model over a one-population dust model, when the latter must capture the total variance in the sky. Each dust population is well characterized by a single-component spectral energy distribution (SED) and accommodates small variations. These populations could signify two distinct components or may originate from a one-component model with different temperatures resulting in different SED scalings. While no spatial information is built into the likelihood, our investigation unveils large-scale spatially coherent structures with high significance, pointing to a physical origin for the observed spatial variation. These results are robust to our choice of likelihood and input data. Furthermore, this spatially varying two-population model is the most favored from Bayesian evidence calculations. Incorporating IRAS 100 μm to constrain the Wein side of the blackbody function, we find the dust populations differ at the 2.5σ level in the spectral index (β d ) versus temperature (T d ) plane. The presence of multiple dust populations has implications for component separation techniques frequently employed in the recovery of the cosmic microwave background.