I have been working on a solver for shallow water equations with reflective boundary conditions. I have found that it diverges very fast. As a workaround I noticed yesterday that if I smooth the momentum (water column height * velocity) field as post processing in each step, the solution no longer blows up. So it seems I have some high frequency error terms appearing somewhere without the smoothing? Where do they usually come from?
Have you plotted the solution using a visualization package or just a cut through the middle of your domain? Sometimes just looking at the solution you're generating will give you an insight into where it's arising from. Based on your description, I'd say there's a mistake in your implementation of your boundary condition.
ETA: The last sentence is probably wrong, and to recap our conversation in the comments and to amplify Gresho & Lee's pithy phrase "the Wiggles are Telling You Something": Even if you haven't used the Galerkin FEM to discretize this, the method is mathematically equivalent to several FVMs and FDMs out there, depending on the quadrature rule used, and so shows similar pathologies if you don't have enough mesh in the boundary layers.