I've recently finished an introductory course on the finite element method from a more mathematical perspective (following Brenner and Scott) and we were introduced to the finite element mass matrix in elliptic problems as the matrix arising from terms without a derivative. For example, a one-dimensional Helmholtz type equation with the form
$$ -u(x)'' + au(x) = f(x), \quad 0 < x < 1, \quad a>0\\ u(0) = u(1) = 0 $$
has a corresponding weak formulation that requires us to find $u$ such that
$$ \int_0^1 u' v'dx + a\int_0^1uv dx = \int_0^1fvdx \quad \forall v \in H^1_0 $$
where $H^1_0 = \{v \in H^1 : v(0) = v(1) = 0\}$.
Choosing $S \subset H^1_0$ to be a conforming finite dimensional subset with a basis $\{ \phi_i \}_{i=1}^N$, and saying $u = \sum_{j = 0}^N u_j \phi_j $, we get the linear problem
$$ (\pmb{K} + a \pmb{M})U = F $$
where $K_{ij} = \int_0^1 \phi_i' \phi_j' dx$ is the stiffness matrix and $M_{ij} = \int_0^1 \phi_i \phi_j dx$ is the mass matrix. The finite element method typically proceeds by choosing $S$ to be the space of piecewise polynomials for example. This formulation extends naturally to higher dimensions.
From this previous post: How to formulate lumped mass matrix in FEM, there are various ways to lump the mass matrix. For example, by summing the off-diagonal terms: $M_{ii} = \sum M_{ij}$.
My question is what is the justification of this? Is there mathematical reasoning why this should give a consistent method? Is there a way to quantify the error introduced by doing this? I've seen an explanation that justifies mass matrix lumping in the context of mechanics where this assumption implies that the mass of the system is concentrated at discrete points, but how does this generalize to more general elliptic PDE problems?