  We introduce a time discretization for Wasserstein gradient flows based on
the classical Backward Differentiation Formula of order two. The main building
block of the scheme is the notion of geodesic extrapolation in the Wasserstein
space, which in general is not uniquely defined. We propose several possible
definitions for such an operation, and we prove convergence of the resulting
scheme to the limit PDE, in the case of the Fokker-Planck equation. For a
specific choice of extrapolation we also prove a more general result, that is
convergence towards EVI flows. Finally, we propose a variational finite volume
discretization of the scheme which numerically achieves second order accuracy
in both space and time.
