The objective of this research is to develop a model that will adequately simulate the dynamics of tsunami propagating across the continental shelf. In practical terms, a large spatial domain with high resolution is required so that source areas and runup areas are adequately resolved. Hence efficiency of the model is a major issue. The three-dimensional Reynolds averaged Navier-Stokes equations are depth-averaged to yield a set of equations that are similar to the shallow water equations but retain the non-hydrostatic pressure terms. This approach differs from the development of the Boussinesq equations where pressure is eliminated in favour of high-order velocity and geometry terms. The model gives good results for several test problems including an oscillating basin, propagation of a solitary wave, and a wave transformation over a bar. The hydrostatic and non-hydrostatic versions of the model are compared for a large-scale problem where a fault rupture generates a tsunami on the New Zealand continental shelf. The model efficiency is also very good and execution times are about a factor of 1.8 to 5 slower than the standard shallow water model, depending on problem size. Moreover, there are at least two methods to increase model accuracy when warranted: choosing a more optimal vertical interpolation function, and dividing the problem into layers. Copyright (c) 2005 John Wiley & Sons, Ltd.