Shock Capturing Methods for Shallow Water Equations

Question

I am looking for some help finding a numerical solution to the shallow water equations:

$\partial_tu+\partial_x(u^2/2+g\eta)=0$

$\partial_t \eta+\partial_x(u\eta)=0$.

where $u$ is the depth averaged horizontal velocity, $\eta$ is the depth, and $g$ is a constant.

Using mathematica's NDSolve I quickly develop a compressive shock (and consequentially the solution begins to ring and diverge) and I need a new solver to capture it accurately. It is my current understanding that I can capture the shock in at least three ways:

(a) I can capture the shock to arbitrary precision and guarantee uniqueness by using an artificial viscosity term. For this I have two questions:

1. What type of viscosity term should I use ($\partial_x^{2n} u$)? And what are the advantages and disadvantages.

2. Should I build a finite difference scheme and if so, what?

(b) I can capture the shock by having a flag that triggers when the gradient is large which then uses the Rankine–Hugoniot conditions to evolve the shock.

1. How do I decide on a trigger for the shock?

2. Should I build a finite difference scheme and if so, what?

(c) I can use a finite difference scheme that has enough numerical diffusion to not run into issues.

1. What scheme and how do I analyze the numerical diffusion contribution.

I don't know where to start. If you have any insight (maybe a solver on GitHub or a critique of that question itself) please share! Any advice on resources or methods would be helpful.