We study the consistency of the estimator in spatial regression with partial differential equa-tion (PDE) regularization. This new smoothing technique allows to accurately estimate spatial fields over complex two-dimensional domains, starting from noisy observations; the regularizing term involves a PDE that formalizes problem specific information about the phenomenon at hand. Differently from classical smoothing methods, the solution of the infinite-dimensional estimation problem cannot be computed analytically. An approximation is obtained via the finite element method, considering asuitable triangulation of the spatial domain. We first consider the consistency of the estimator in the infinite-dimensional setting. We then study the consistency of the finite element estimator, resulting from the approximated PDE. We study the bias and variance of the estimators, with respect to the sample size and to the value of the smoothing parameter. Some final simulation studies provide numerical evidence of the rates derived for the bias, variance and mean square error.