  This paper is concerned with the recovery of (approximate) solutions to
parabolic problems from incomplete and possibly inconsistent observational
data, given on a time-space cylinder that is a strict subset of the
computational domain under consideration. Unlike previous approaches to this
and related problems our starting point is a regularized least squares
formulation in a continuous infinite-dimensional setting that is based on
stable variational time-space formulations of the parabolic PDE. This allows us
to derive a priori as well as a posteriori error bounds for the recovered
states with respect to a certain reference solution. In these bounds the
regularization parameter is disentangled from the underlying discretization. An
important ingredient for the derivation of a posteriori bounds is the
construction of suitable Fortin operators which allow us to control oscillation
errors stemming from the discretization of dual norms. Moreover, the
variational framework allows us to contrive preconditioners for the discrete
problems whose application can be performed in linear time, and for which the
condition numbers of the preconditioned systems are uniformly proportional to
that of the regularized continuous problem.
  In particular, we provide suitable stopping criteria for the iterative
solvers based on the a posteriori error bounds. The presented numerical
experiments quantify the theoretical findings and demonstrate the performance
of the numerical scheme in relation with the underlying discretization and
regularization.
