- To control computational instabilities associated with
poor numerical design. In order to do this, the coefficient of
viscosity is usually made unphysically large.
*I make it a point not to do this for my modelling. In the context of DieCAST, it must be noted that an artificial viscosity is used to stabilize the advection scheme --- which has the effect of making the advection scheme formally less than fourth-order accurate, although it is still accurate. I prefer to use the QUICK advection scheme which is 3rd-order accurate and can be generalized to higher order if required.* - Some people think that horizontal eddy-viscosity is a physically
meaningful way to parameterize sub-grid-scale phenomena.
*I know that this is just not so. For some reason oceanographers have not been able to figure out that the sub-grid-scale is a function of truncation error. Instead they have gotten all confused by the work of Osborn Reynolds, which has little to do with computational methods. Or perhaps they think that that the Munk western boundary-layer is physically meaningful. It isn't. Munk was a smart guy but he was wrong on this occasion.*

A control-volume model (like DieCAST) might
define u_{i} as a component of velocity averaged
over the volume of the i'th cell. (Many oceanographers
crudely think of the
index i indicating the center of the cell --- but this is only second-order
accurate thinking.)
In order to accurately calculate the viscosity terms one must
compute viscous fluxes at the cell faces at each side of the
cell

K u'_{i+0.5}

and

K u'_{i-0.5}

Here I have used the notation +/-0.5 to indicate locations to the right/left
of the cell center. The prime ' indicates the first derivative.
A double prime '' wil be used to indicate a second derivative.
Molecular viscosity is denoted K.

So the problem reduces to computing the derivative of u at the cell face
from the cell-averaged values of u.
This may seem a little messy because we are dealing with both face and
cell-averaged quantities.
Another difficulty arises when there are boundaries.
I will use integral reconstruction to make things simple by replacing
algebra with simple geometry. Consider u_{i} as the cell-averaged
quantity for i=0,1,2,3,4, ... I cells. If the cell i=0
is on land it has no velocity u_{0}=0.
** Denote the cell-with as D**

- First one cumulatively sums to exactly obtain the integral of u at the
cell faces. This involves no truncation error. Denote the integral
of u with the symbol S (S for summation). Thus we obtain:

S_{-0.5}= 0

S_{0.5}= u_{0}D

S_{1.5}=(u_{0}+u_{1})D

S_{2.5}=(u_{1}+ u_{2})D

S_{3.5}=(u_{1}+ u_{2}+ u_{3})D

S_{4.5}=(u_{1}+ u_{2}+ u_{3}+ u_{4})D - Next we recognise that the second derivative of f
_{i+/-0.5}is merely the first derivative of u_{i+/-0.5}, ie u'_{i+/-0.5}= S''_{i+/-0.5}. Consider a specific instance where we seek u'_{i+/-0.5}with i=1, corresponding to the boundary face of the cell.

u'_{0.5}= S''_{0.5}= a S_{-0.5}+ b S_{0.5}+ c S_{1.5}+d S_{2.5}+ e S_{3.5}

We can substitute for S to obtain

u'_{0.5}= b u_{0}D

+ c u_{0}D + c u_{1}D

+ d u_{0}D + d u_{1}D + d u_{2}D

+ e u_{0}D + e u_{1}D + e u_{2}D + e u_{3}D

which can be rewritten

u'_{0.5}= (b+c+d+e) u_{0}+ (c+d+e) u_{1}D + (d+e) u_{2}D + e u_{3}D

The Matlab script coef.m can be used to calculate values for a,b,c,d,e to represent a second derivative at whatever order of accuracy is required at whichever point the second derivative is required.- Considering the second derivative S''
_{0.5}at second-order accuracy we use: x=[-0.5:1:1.5]; xt=0.5; C=coef(x,xt); C(3,:) to obtain the stencil:

(a,b,c,d,e)=(1,-2,1,0,0)/D^{2}

so that u'_{0.5}= (- u_{0}+ u_{1})/D which is the usual second-order accurate differencing.The boundary condition u

_{0}=0 gives the viscous flux at the boundary:

K u'_{i+0.5}= K u_{1}/D - Considering the second derivative S''
_{0.5}at fourth-order accuracy we use: x=[-0.5:1:3.5]; xt=0.5; C=coef(x,xt); 12*C(3,:) to obtain:

(a,b,c,d,e)=(11, -20, 6, 4, -1 )/(12D^{2})

so that u'_{0.5}=(-11u_{0}+9u_{1}+3u_{2}-u_{3}) /(12 D)The boundary condition u

_{0}=0 gives the viscous flux at the boundary:

K u'_{i+0.5}= K (9u_{1}+ 3u_{2}-u_{3})/DSimilarly, if we wish to find u'

_{1.5}then we use: x=[-0.5:1:3.5]; xt=1.5; C=coef(x,xt); 12*C(3,:)

(a,b,c,d,e)=(-1, 16, -30, 16, -1)/(12D^{2})

so that u'_{0.5}=(u_{0}-15u_{1}+15 u_{2}-u_{3}) /(12 D)The boundary condition u

_{0}=0 gives the viscous flux one cell in from the boundary:

K u'_{i+0.5}= K (-15u_{1}+ 15u_{2}-u_{3})/DFurther away from the boundary, similar analysis gives

so that u'_{1.5}=(u_{1}-15u_{2}+15u_{3}-u_{4}) /(12 D)If u

_{5}=0 because it is on a land point then x=[1.5:1:5.5]; xt=4.5; C=coef(x,xt); 12*C(3,:)

(a,b,c,d,e)=(-1, 4, 6, -20, 11 )/(12D^{2})

so that u'_{4.5}=(u_{2}-3u_{3}-9u_{4}+11u_{5}) /(12 D)Similarly, x=[1.5:1:5.5]; xt=3.5; C=coef(x,xt); 12*C(3,:)

(a,b,c,d,e)=(-1, 16, -30, 16, -1)/(12D^{2})

so that u'_{3.5}=(u_{2}-15u_{3}+15u_{4}-u_{5}) /(12 D)

- Considering the second derivative S''