Summary: We present a numerical method for simulations of nonlinear surface water waves over variable bathymetry. It is applicable to either two- or three-dimensional flows, as well as to either static or moving bottom topography. The method is based on the reduction of the problem to a lower-dimensional Hamiltonian system involving boundary quantities alone. A key component of this formulation is the Dirichlet-Neumann operator which, in light of its joint analyticity properties with respect to surface and bottom deformations, is computed using its Taylor series representation. We present new, stabilized forms for the Taylor terms, each of which is efficiently computed by a pseudospectral method using the fast Fourier transform. Physically relevant applications are displayed to illustrate the performance of the method; comparisons with analytical solutions and laboratory experiments are provided.