  In this paper, a least-squares finite element method for scalar nonlinear
hyperbolic balance laws is proposed and studied. The approach is based on a
formulation that utilizes an appropriate Helmholtz decomposition of the flux
vector and is related to the standard notion of a weak solution. This
relationship, together with a corresponding connection to negative-norm
least-squares, is described in detail. As a consequence, an important numerical
conservation theorem is obtained, similar to the famous Lax-Wendroff theorem.
The numerical conservation properties of the method in this paper do not fall
precisely in the framework introduced by Lax and Wendroff, but they are similar
in spirit as they guarantee that when $L^2$ convergence holds, the resulting
approximations approach a weak solution to the hyperbolic problem. The
least-squares functional is continuous and coercive in an $H^{-1}$-type norm,
but not $L^2$-coercive. Nevertheless, the $L^2$ convergence properties of the
method are discussed. Convergence can be obtained either by an explicit
regularization of the functional, that provides control of the $L^2$ norm, or
by properly choosing the finite element spaces, providing implicit control of
the $L^2$ norm. Numerical results for the inviscid Burgers equation with
discontinuous source terms are shown, demonstrating the $L^2$ convergence of
the obtained approximations to the physically admissible solution. The
numerical method utilizes a least-squares functional, minimized on finite
element spaces, and a Gauss-Newton technique with nested iteration. We believe
that the linear systems encountered with this formulation are amenable to
multigrid techniques and combining the method with adaptive mesh refinement
would make this approach an efficient tool for solving balance laws (this is
the focus of a future study).
