  A novel implicit integration scheme for the Discrete Element Method (DEM)
based on the variational integrator approach is presented. The numerical solver
provides a fully dynamical description that, notably, reduces to an energy
minimisation scheme in the quasi-static limit. A detailed derivation of the
numerical method is presented for the Hookean contact model and tested against
an established open source DEM package that uses the velocity-Verlet
integration scheme. These tests compare results for a single collision,
long-term stability and statistical quantities of ensembles of particles.
Numerically, the proposed integration method demonstrates equivalent accuracy
to the velocity-Verlet method.
