A simple time-evolution method for nuclear quantum dynamics is proposed in which the wavefunction is expanded in terms of Gaussians placed at von Neumann unit cells in phase space, with an accuracy threshold parameter epsilon. The wavepacket propagation is obtained by solving the equations of motion for expansion coefficients only. New Gaussians are added to the basis set when the change in wavefunction exceeds another threshold parameter delta (>epsilon). From the application to model tunneling systems in one to three dimensions, the present Gaussian expansion approach is able to accurately evaluate the tunneling time and tunneling direction in multidimensional systems. (C) 2018 Elsevier B.V. All rights reserved.