  A numerical method for approximating weak solutions of an aggregation
equation with degenerate diffusion is introduced. The numerical method consists
of a stabilized finite element method together with a mass lumping technique
and an extra stabilizing term plus a semi--implicit Euler time integration.
Then we carry out a rigorous passage to the limit as the spatial and temporal
discretization parameters tend to zero, and show that the sequence of finite
element approximations converges toward the unique weak solution of the model
at hands. In doing so, nonnegativity is attained due to the stabilizing term
and the acuteness on partitions of the computational domain, and hence a priori
energy estimates of finite element approximations are established. As we deal
with a nonlinear problem, some form of strong convergence is required. The key
compactness result is obtained via an adaptation of a
Riesz--Fr\'echet--Kolmogorov criterion by perturbation. A numerical example is
also presented.
