In analogy with similar effects in adiabatic compressible fluid dynamics, the effects of buoyancy gradients on incompressible stratified flows are said to be "thermal." The thermal rotating shallow water (TRSW) model equations contain three small nondimensional parameters. These are the Rossby number, the Froude number, and the buoyancy parameter. Asymptotic expansion of the TRSW model equations in these three small parameters leads to the deterministic thermal versions of the Salmon's L1 (TL1) model and the thermal quasi-geostrophic (TQG) model, upon expanding in the neighborhood of thermal quasi-geostrophic balance among the flow velocity and the gradients of free surface elevation and buoyancy. The linear instability of TQG at high wavenumber tends to create circulation at small scales. Such a high-wavenumber instability could be unresolvable in many computational simulations, but its presence at small scales may contribute significantly to fluid transport at resolvable scales. Sometimes, such effects are modeled via "stochastic backscatter of kinetic energy." Here, we try another approach. Namely, we model "stochastic transport" in the hierarchy of models TRSW/TL1/TQG. The models are derived via the approach of stochastic advection by Lie transport (SALT) as obtained from a recently introduced stochastic version of the Euler-Poincaré variational principle. We also indicate the potential next steps for applying these models in uncertainty quantification and data assimilation of the rapid, high-wavenumber effects of buoyancy fronts at these three levels of description by using the data-driven stochastic parametrization algorithms derived previously using the SALT approach.