  A fully coupled system of two second-order parabolic degenerate equations
arising as a thin film approximation to the Muskat problem is interpreted as a
gradient flow for the 2-Wasserstein distance in the space of probability
measures with finite second moment. A variational scheme is then set up and is
the starting point of the construction of weak solutions. The availability of
two Liapunov functionals turns out to be a central tool to obtain the needed
regularity to identify the Euler-Lagrange equation in the variational scheme.
