Driven by growing momentum in two-dimensional geophysical flow modeling, this paper introduces a general family of "thermal" rotating shallow-water models. The models are capable of accommodating thermodynamic processes, such as those acting in the ocean mixed layer, by allowing buoyancy to vary in horizontal position and time as well as with depth, in a polynomial fashion up to an arbitrary degree. Moreover, the models admit Euler-Poincare variational formulations and possess Lie-Poisson Hamiltonian structure. Such a geometric property provides solid fundamental support to the theories described with consequences for numerical implementation and the construction of unresolved motion parametrizations. In particular, it is found that stratification halts the development of small-scale filament rollups recently observed in a popular model, which, having vertically homogeneous density, represents a special case of the models presented here.