Redi diffusion taper bugs: spurious dianeutral mixing

[aekiss]

ACCESS-OM2 uses Redi tracer diffusion at 1° and 0.25°, with the default dm_taper Redi tapering scheme based on Danabasoglu and McWilliams (1995) Eqn A.7a.

Fig 23.4 in sec 23.5.2.1 of Griffies (2012) shows the intended scaling of the Redi diffusivity using this tapering method in regions of large isopycnal slope (blue curve):

Griffies (2012) states that this tapering can be applied at all points, but the implementation in MOM5 has an if-statement that only applies the taper when the slope (absslope) exceeds smax (see x-flux code here in subroutine fx_flux_ndiffuse, and similar for y-flux in subroutine fy_flux_ndiffuse). Consequently as the slope increases past smax, taperA abruptly changes from 1.0 to 0.5 (because the argument of tanh is zero at that point), as shown in the red curve above.

              ! taper for steep slope regions 
              if(absslope > smax) then 
                  taperA = 0.5*(1.0 + tanh(smax_swidthr-absslope*swidthr))             
                  taperA = dm_taper_const*taperA + gkw_taper_const*(smax/absslope)**2
              else 
                  taperA = 1.0
              endif

where

swidthr       = 1.0/(swidth + epsln) !for slope taper function when dm_taper used 
smax_swidthr  = smax*swidthr         !for slope taper function when dm_taper used 
real    :: dm_taper_const   = 1.0      ! internally set to unity when dm_taper=.true.
real    :: gkw_taper_const  = 0.0      ! internally set to unity when gkw_taper=.true.

In regions with slope close to smax, this abrupt change in diffusivity might produce sharp spatial discontinuities in Redi flux, and therefore very large, transient flux divergences.

Have I misunderstood something here?

[rmholmes]

Good catch @aekiss. I've looked through the code a bit and I can't see anything wrong with your argument.

I don't recall seeing anything funny when looking at the Redi heat budget terms. It seems to me the model might be able to absorb this kind of sharp flux convergence issue without going unstable. It would be easy enough to remove the if statement and see if it makes a big difference.

[aekiss]

If I understand it correctly, the Danabasoglu and McWilliams (1995) version is also discontinuous, but in the oposite way, with the tanh function only applied to slopes smaller than a parameter $S_{max}$ that can be set independently of $S_c$ (with $S_c$ corresponding to smax in MOM5):

[aekiss]

The MOM5 implementation also calculates tapers for x and y diffusivities separately, in subroutines fx_flux_ndiffuse and fy_flux_ndiffuse, each of which only uses the slope in that direction, not the total slope as specified in Danabasoglu and McWilliams (1995) and Griffies (2012). So the taper factors (taperA) will be different in the x and y directions, producing anisotropic diffusivity in regions where the tapering is active in at least one direction.

Can this anisotropy rotate the diffusive flux vector out of the local neutral surface, causing diapycnal fluxes? Seems to me it can, but I hope I'm wrong!

[aekiss]

The MOM5 implementation of gkw_taper (quadratic tapering of Gerdes, Koberle, and Willebrand) suffers from the same anisotropy issue as dm_taper, but not the discontinuity issue, because gkw_taper does need an if-statement so that it applies only for slopes exceeding smax.

Looks like the if-statement issue with dm_taper may have been due to adding dm_taper code into existing logic for gkw_taper?

[StephenGriffies]

Here are some comments, based on incomplete memory of the work I did with this scheme.

Regions of neutral slopes larger than smax tend to be much much larger than smax. That is, in a coarse model the slopes tend to jump to a large value relatively quickly, so that there are few (if any) regions with slopes close to, but below, slopemax. This experience might be modified with finer grid spacing models, such as 0.25 degree.

Relying on tanh without the slopemax: I recall we formerly did that, but perhaps introduced slopemax for safety. Memory fails me on that question.

For the anisotropy point: I concur; we did do directions separately. We held numerical stability very central to the approach, so erred on the side of stability by checking each direction separately. This approach could be cleaned up, but it would require merging code since the triad scheme operates in each direction separately.

Are you seeing any simulation feature that led you to look into the details of the MOM5 neutral diffusion scheme?

[aekiss]

Thanks for the background @StephenGriffies, it's good to know the slope doesn't linger near smax much, at least at 1°.

I'm not aware of any model misbehaviour resulting from the way this is implemented (though I also haven't specifically looked for it). It's just something I noticed when trying to work out how to calculate Redi diffusion offline.

Re. the if-statement: the CFL linear stability condition for neutral diffusion suggests that applying tanh unconditionally should make the model a little more stable, because the diffusivity will be reduced for slopes near but below smax. This would be easy to implement and test.

Re. anisotropy and resulting dianeutral flux, having looked at the code more closely I see your point that fx_flux_ndiffuse and fy_flux_ndiffuse would need to be merged somehow, because the x and y slopes are calculated using different offsets (ip and jq, respectively), so it doesn't seem straightforward to simply calculate both slopex and slopey in each of fx_flux_ndiffuse and fy_flux_ndiffuse and use absslope = sqrt(slopex**2 + slopey**2) to calculate matching taperA in both subroutines. That said, if the spurious dianeutral flux is significant it might be worth trying. Also, using absslope = sqrt(slopex**2 + slopey**2) to calculate the taper in both directions can only reduce diffusivity and therefore increase stability relative to the current scheme.