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.
