We introduce a one-parameter family of massive Laplacian operators $(Δ^{m(k)})_{k\in[0,1)}$ defined on isoradial graphs, involving elliptic functions. We prove an explicit formula for the inverse of $Δ^{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 $(Δ^{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.