  We construct a finite element method for the numerical solution of a
fractional porous medium equation on a bounded open Lipschitz polytopal domain
$\Omega \subset \mathbb{R}^{d}$, where $d = 2$ or $3$. The pressure in the
model is defined as the solution of a fractional Poisson equation, involving
the fractional Neumann Laplacian in terms of its spectral definition. We
perform a rigorous passage to the limit as the spatial and temporal
discretization parameters tend to zero and show that a subsequence of the
sequence of finite element approximations defined by the proposed numerical
method converges to a bounded and nonnegative weak solution of the
initial-boundary-value problem under consideration. This result can be
therefore viewed as a constructive proof of the existence of a nonnegative,
energy-dissipative, weak solution to the initial-boundary-value problem for the
fractional porous medium equation under consideration, based on the Neumann
Laplacian. The convergence proof relies on results concerning the finite
element approximation of the spectral fractional Laplacian and compactness
techniques for nonlinear partial differential equations, together with
properties of the equation, which are shown to be inherited by the numerical
method. We also prove that the total energy associated with the problem under
consideration exhibits exponential decay in time.
