Abstract. 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. 1. Introduction. Accurate modeling of surface water wave dynamics over bottom topography is of great importance to coastal engineers and has drawn considerable attention in recent years. Here are just a few applications: linear [16,33] and nonlinear wave shoaling [25,24,28,27], Bragg reflection [35,15,36,14,30,33,4], harmonic wave generation [3,17,18], and tsunami generation [37,22]. As the publications cited above attest, the character of coastal wave dynamics can be very complex: when entering shallow water, waves are strongly affected by the bottom through shoaling, refraction, diffraction, and reflection. Nonlinear effects related to wave-wave and wave-bottom interactions can cause wave scattering and depth-induced wave breaking. In turn, the resulting nonlinear waves can have a great influence on sediment transport and the formation of sandbars in nearshore regions. The presence of bottom topography also introduces additional space and time scales to the classical perturbation problem (see, e.g., [10]).Traditionally, the water wave problem on variable depth has been modeled using long-wave approximations such as the Boussinesq or shallow water equations [26,40]. See also [39] for a weakly nonlinear formulation for water waves over topography. In recent years, progress in both mathematical techniques and computer power has led to a rapid development of numerical models solving the full Euler equations. These can be divided into two main categories: boundary integral methods [2,24,49,7,20,27] and spectral methods. In particular, efficient spectral methods based on perturbation expansions have been developed for the computation of water waves on constant or