Large-scale tsunami propagation simulations from the fault region to the coast are conducted using a three-dimensional (3-D) parallel unstructured mesh finite element code (Fluidity-ICOM). Unlike conventional 2-D approximation models, our tsunami model solves the full 3-D incompressible Navier-Stokes (NS) equations. The model is tested against analytical solutions to simple dispersive wave propagation problems. Comparisons of our 3-D NS model results with those from linear shallow water and linear dispersive wave models demonstrate that the 3-D NS model simulates the dispersion of very short wavelength components more accurately than the 2-D models. This improved accuracy is achieved using only a small number (three to five) of vertical layers in the mesh. The numerical error in the wave velocity compared with the linear wave theory is less than 3% up to kH = 40, where k is the wave number and H is the sea depth. The same 2-D and 3-D models are also used to simulate two earthquake-generated tsunamis off the coast of Japan: the 2004 off Kii peninsula and the 2011 off Tohoku tsunamis. The linear dispersive and NS models showed good agreement in the leading waves but differed especially in their near-source, short wavelength dispersive wave components. This is consistent with the results from earlier tests, suggesting that the 3-D NS simulations are more accurate. The computational performance on a parallel computer showed good scalability up to 512 cores. By using a combination of unstructured meshes and high-performance computers, highly accurate 3-D tsunami simulations can be conducted in a practical timescale.