  We introduce a one-parameter family of massive Laplacian operators
$(\Delta^{m(k)})_{k\in[0,1)}$ defined on isoradial graphs, involving elliptic
functions. We prove an explicit formula for the inverse of $\Delta^{m(k)}$, the
massive Green function, which has the remarkable property of only depending on
the local geometry of the graph, and compute its asymptotics. We study the
corresponding statistical mechanics model of random rooted spanning forests. We
prove an explicit local formula for an infinite volume Boltzmann measure, and
for the free energy of the model. We show that the model undergoes a second
order phase transition at $k=0$, thus proving that spanning trees corresponding
to the Laplacian introduced by Kenyon are critical. We prove that the massive
Laplacian operators $(\Delta^{m(k)})_{k\in(0,1)}$ provide a one-parameter
family of $Z$-invariant rooted spanning forest models. When the isoradial graph
is moreover $\mathbb{Z}^2$-periodic, we consider the spectral curve of the
characteristic polynomial of the massive Laplacian. We provide an explicit
parametrization of the curve and prove that it is Harnack and has genus $1$. We
further show that every Harnack curve of genus $1$ with
$(z,w)\leftrightarrow(z^{-1},w^{-1})$ symmetry arises from such a massive
Laplacian.
