Calculate numerical diffusivity¶
This module has functions to calculate the numerical diffusivity experienced by a tracer, associated to a specific configuration of the MITgcm. In particular, it was developed to calculate the equivalent diffusivity $kappa$, defined (here) as $kappa = kappa_{pres}+kappa_{num}$, where $kappa_{pres}$ is the prescibed or explicit tracer diffusivity one imposes on the model and $$k_{num}$$ is the additional diffusivity due to numerical truncation errors. Note that there are two $$kappa_{pres}$$ and therefore two $kappa$, one for the horizontal dimensions and one for the vertical one.
These calculations try to reproduce the method used by [1] Abernathy et al. 2010, [2] Hill et al. 2011, and [3] Leibensperger and Plumb, 2013 to determine the numerical diffusivity MITgcm Southern Ocean configurations [1,2] and a baroclinic flow simulation simulation [3].
The idea is that from the evolution equation for the variance of the tracer concentration in the model output
begin{equation} frac{1}{2}frac{partial{overline{q^{2}}}}{partial{t}}=-kappa_{h} overline{|nabla_h q|^2}-kappa_{v} overline{(frac{partial{q}}{partial {z}})^{2}} end{equation}
one can fit by a least squares regression, suitable values of $kappa_h$ and $kappa_v$ that satisfy the equation.
Calculate the volume of the domain (function: CalcDomVolume)¶
The volume of a tracer cell (remember we have an Arakawa C grid, so this changes depending on which kind of cell we are thinking about) is given by
$V(i,j,k)=depth times area = (hfacC(i,j,k)times dRf(k)) times rA(i,j) = (hfacC(i,j,k)times dRf(k)) times dXg(i,j) times dYg(i,j)$,
where hfacC is the fraction of the cell that is open (not occupied with land). So, the total volume of the domain is
$sumlimits_{i=1}^{nx}{sumlimits_{j=1}^{ny}{sumlimits_{k=1}^{nz}{(hfacC(i,j,k)times dRf(k)) times rA(i,j)}}}$
1st Term: The volume-weighted average of the squared concentration (function: CalcVariance, CalcTimeDer)¶
The first term in the variance evolution equation is $frac{1}{2}frac{partial{overline{q^{2}}}}{partial{t}}$. Note that we care about the time derivative of the variance, so that the mean concentration that usually appears in the definition of variance will not play a role here, since it is constant in time (we are not putting in or letting out any tracer).
We are going to calculate $overline{q^2}$, the volume-weighted average of the squared concentration, and then the time derivative of that using a centered difference scheme.
2nd Term: The volume-weighted average of the squared horizontal gradient (function: CalcAvgHorGrad)¶
The second term in the variance evolution equation is $-kappa_{h} overline{|\nabla_h q|^2}$. Next, we calculate the square of the horizontal gradient $|nabla_h q|^2=(frac{partial{q}}{partial{x}})^2+(frac{partial{q}}{partial{y}})^2$.
Spatial derivatives are approximated using a centered-difference scheme.
3rd Term: The volume-weighted average of the squared vertical derivative (function: CalcAvgVerGrad)¶
The third term in the variance evolution equation is $-kappa_{v} overline{(frac{partial{q}}{partial{z}})^2}$. Next, we calculate the square of the vertical gradient $(frac{partial{q}}{partial{z}})^2$.
The vertical derivative is approximated using a centered-difference scheme.