I'm trying to numerically solve a system of Shallow Water Equations in Wolfram Mathematica 14.0.0 on the rotating sphere by using the NDSolve method with a purpose of earthquake-generated Tsunami simulation. I define necessary physical constants, initial ocean floor profile, initial values of angular speeds and water surface. I also setup time-dependent ocean floor excitation law and atmospheric pressure function. My code is quite huge, so I’ve uploaded it on my Google Drive.

However, at some point of simulation I always receive an error "At t == ~18000, step size is effectively zero; singularity or stiff system suspected". So, when I try to visualize my solution, I also see that it tends to rapidly burst after some iterations. One could see my problem at the animated video below:

As it was recommended in the post, I've tried to use “Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"}}”, but I received another error “As specified the PDE cannot be discretized using the TensorProductGrid method and will be discretized using the FiniteElement method instead. Sometimes reformulating a PDE allows other methods to be used.”. Playing with the form of my PDE didn’t give any results.

Then I’ve googled some other different options in “NDSolveValue” function, such as “Method -> {PDEDiscretization}”, but none of them helped me to stabilize calculation process near the boundary of the simulation region. At the small simulation times I can clearly see my wave propagating correctly. In contrast, at the large times numerical interpolation solution explodes with enormous water surface height values without any reason.

So, give me, please, some advice how to deal with such issue. I’m out of ideas ☹