Improve accuracy and stability of tridiagonal solve for vertical mixing

[cbegeman]

There is a potential accuracy and stability issue with our current implementation when diffusivities are high for one layer and low for a neighboring layer but their thicknesses are of similar magnitude. This issue is documented in Schopf and Loughe, 1995, Appendix E.

Both MOM6 and POP implement the fix noted in this publication, but MPAS-Ocean uses the older algorithmic approach to the tridiagonal solve. We have not yet, as far as I am aware, noticed issues with MPAS-Ocean simulations that are clearly linked with this issue, but some vmix schemes may generate large enough diffusivities to cause this issue.

Thanks go to Alistair Adcroft and Bob Hallberg for prompting us to check this.

[cbegeman]

@vanroekel Your input on this would be appreciated.

[cbegeman]

@mark-petersen We discussed this issue elsewhere. You're invited to elaborate here.

[mark-petersen]

@cbegeman, thanks for bringing this to our attention. When I implemented the tridiagonal solve, it was just a straight version that I derived myself. From my 2012 design document:

I did not use the Formulation in Appendix E of the Schopf and Loughe 1995 paper. But I see that POP does.
You can grab the POP ref manual here https://www2.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf
When I wrote the MPAS version, I was looking at the POP manual here:

you can see the similarities with my write-up. Just dz -> h

But there is a section further down with this topic exactly, and I did not use that part. I simply didn’t look for an additional section in the POP manual, and the topic didn't come up in the design or code review. Here it is in the POP manual:

It looks like the variables are simply redefined with some scaling coefficients to keep them closer to order 1, so you don’t lose digits when adding numbers of widely varying size. The 1995 paper says this is important for vertical convection. And KPP does turn on high values of vertical diffusivity in unstable regions.

It sounds like it is important, since it was in POP and MOM (presumably) uses it. The consequence of not using it seems to be that high values of diffusivity (kappa) produce inaccurate solutions, like it doesn’t actually mix correctly. I have not seen that symptom in MPAS-Ocean specifically, but it could be causing some of our other problems like our strangely high heat uptake, and there might be follow-on consequences with AMOC and sea ice. (Easy to just list all your problems, of course). I think this is important enough that we need to discuss it and consider adding this change.

[vanroekel]

Thanks for noting this @cbegeman and for the detailed comment @mark-petersen. Given this is a lack of accuracy and stability at high diffusivity I can certainly see how this could be potentially important. But it's not super clear what it might do. I wonder if it would be somewhat straightforward to code the current tri-diagonal solve and the revised algorithm then run a temperature profile with different diffusivities to see what the sensitivity is. This might be a low cost way to see if this is important. The comment at the end of Appendix E saying 64 bit arithmetic is sufficiently accurate over 10 orders of magnitude for diffusivity suggests this may not be a critical fix. But given Bob and Alistair mentioned it it seems very worthy to look into. @cbegeman did Bob or Alistair happen to mention any impacts of this change?

[cbegeman]

@vanroekel Bob and Alistair mentioned that the original algorithm led to instability at ice-shelf fronts. I believe this had arisen not only from differences in diffusivity but also in layer thickness between vertical layers. We're not currently running with large differences in layer thickness vertically, so the main question as I see it is how often we generate large differences in diffusivity between vertical levels in realistic simulations. Do we have the output to evaluate that?

[vanroekel]

Ah thanks for the explanation. In my experience, we have large diffusivity differences extremely frequently, almost every simulation, from 0 to 1 for interior convection. Any monthly average file with timeMonthly_avg_vertDiffTopOfCell (should be standard). It likely won't go zero to 1 as it gets averaged, but you should see very big, noisy swings. Which as I write this and think about this, I wonder if that near grid scale diffusivity shifting I see could be due (in part) to the issue raised here.

I also want to say that I'm not at all against this change. I'd just like to get a sense of how important the change might be for some of our simulations. E.g. does it help any of the issues @mark-petersen notes in his longer comment.

[cbegeman]

In my experience, we have large diffusivity differences extremely frequently, almost every simulation, from 0 to 1 for interior convection.

@vanroekel What do you mean "0 to 1"? I thought we were talking about orders of magnitude differences in diffusivity?

[cbegeman]

@mark-petersen and @vanroekel Would you like me to look into this or is this something one of you have a particular interest in?

[vanroekel]

@vanroekel What do you mean "0 to 1"? I thought we were talking about orders of magnitude differences in diffusivity?

@cbegeman for vertical diffusivity I haven't seen much beyond 10-20 m^2/s for diffusivity. Those are very large values for vertical mixing

[vanroekel]

On your other question if you have the bandwidth and are interested in digging into this @cbegeman that would be great in my view.