Treated traditionally by the Ehrenfest approximation, the dynamics of a one-dimensional molecular crystal model with off-diagonal exciton-phonon coupling is investigated in this work using the Dirac-Frenkel time-dependent variational principle with the multi-D2 Ansatz. It is shown that the Ehrenfest method is equivalent to our variational method with the single D2 Ansatz, and with the multi-D2 Ansatz, the accuracy of our simulated dynamics is significantly enhanced in comparison with the semi-classical Ehrenfest dynamics. The multi-D2 Ansatz is able to capture numerically accurate exciton momentum probability and help clarify the relation between the exciton momentum redistribution and the exciton energy relaxation. The results demonstrate that the exciton momentum distributions in the steady state are determined by a combination of the transfer integral and the off-diagonal coupling strength, independent of the excitonic initial conditions. We also probe the effect of the transfer integral and the off-diagonal coupling on exciton transport in both real and reciprocal space representations. Finally, the variational method with importance sampling is employed to investigate temperature effects on exciton transport using the multi-D2 Ansatz, and it is demonstrated that the variational approach is valid in both low and high temperature regimes.