I'm going to work on the assumption that the disturbances aren't large enough to create shock waves, which are a whole other can of worms.
Your energy functional is close but not quite right -- the Hamiltonian for the shallow water equations is
$$H = \frac{1}{2}\int_\Omega\left(h|u|^2 + gh^2\right)dx,$$
although the Poisson brackets are a little more involved.
You're also correct to identify symplectic integrators as the tool to use if you care about the Hamiltonian structure of the problem.
But it's important to note that a symplectic integrator conserves the energy of some perturbed Hamiltonian
$$H' = H + \delta t\cdot\delta H + \mathscr{O}(\delta t^2),$$
where very often we can read off what $\delta H$ is in terms of a few Poisson brackets using the Baker-Campbell-Hausdorff formula.
In other words, some energy drift is to be expected.
Some time schemes are guaranteed to exactly conserve all linear and quadratic invariants, but the $h|u|^2$ term isn't quadratic.
That said, the property of exactly sampling from the trajectory of the modified Hamiltonian $H'$ is incredibly useful -- you get a lot of properties for free, like the structure of the system near equilibria, bounded orbits, phase space volume conservation, etc.
I haven't worked through the details of symplectic discretization of the shallow water equations specifically, so take this with a grain of salt.
There are several explicit schemes that are only symplectic for simple Hamiltonians that can be written as $H = K(p) + U(q)$ with the usual Poisson bracket.
The Hamiltonian and the bracket for SWE are neither of those.
I think that the implicit midpoint rule
$$\frac{z_{n + 1} - z_n}{\delta t} = f\left(\frac{z_{n + 1} + z_n}{2}\right)$$
for the ODE $\dot z = f(z)$ is symplectic even in unusual cases.
Now you've got an implicit equation to solve, but nothing's free.
So it may be worth trying other timestepping schemes.
Finally, upwinding interacts with the Hamiltonian structure of the problem in ways that are a little bit subtle.
Although they took more of a finite element approach, this paper worked through how to do energy-conserving discretization of SWE with upwinding.
Based on the plot you're showing, you should be asking yourself whether the energy drift that you've found indicates an error or failure to preserve the Hamiltonian structure, or merely the expected drift from integrating a slightly perturbed problem.
While the Hamiltonian that you wrote down isn't completely correct, the plot you showed indicated an energy drift on the order of 0.035 and a total of around 10, so the drift is less than 1% of the total.
That doesn't set off any red flags for me, but if you want to try a few things to debug this, you can:
- Run the simulation again at resolution $\delta x / 2$, $\delta t / 2$ and see what the energy drift looks like.
Any convergent timestepping scheme (symplectic or not) should halve the drift, so this won't tell you everything you need to know.
- Run the simulation for a very long time and see what the energy drift looks like.
To get an idea of what "long" means, you get a speed scale from the wave speed $c = \sqrt{gh}$ and a length scale from the diameter $L$ of the domain, so several multiples of $L / c$ is a good ballpark.
If the solution explodes, your scheme probably isn't symplectic.
I would agree with the suggestion in the comments to just use finite volume methods if the problem you had related to capturing shock waves, but that doesn't appear to be the case.
