  Inspired by recent works on the threshold dynamics scheme for multi-phase
mean curvature flow (by Esedo\={g}lu-Otto and Laux-Otto), we introduce a novel
framework to approximate solutions of the Muskat problem with surface tension.
Our approach is based on interpreting the Muskat problem as a gradient flow in
a product Wasserstein space. This perspective allows us to construct weak
solutions via a minimizing movements scheme. Rather than working directly with
the singular surface tension force, we instead relax the perimeter functional
with the heat content energy approximation of Esedo\={g}lu-Otto. The heat
content energy allows us to show the convergence of the associated minimizing
movement scheme in the Wasserstein space, and makes the scheme far more
tractable for numerical simulations. Under a typical energy convergence
assumption, we show that our scheme converges to weak solutions of the Muskat
problem with surface tension. We then conclude the paper with a discussion on
some numerical experiments and on equilibrium configurations.
