Timestep size, timestepping growth, mesh refinement, the solver, the predictor-corrector scheme can make a difference with the error... There are lots of different ways to reduce the error! Probably the best thing to fix is the mesh refinement, especially in regions of high elevation gradients from my experience, and the type of boundary condition you use could also impact it.