Vorticity Budgets in a C-grid SW Model

In this section, we compare the analytic vorticity budget with the equivalent discretized vorticity budget for a C-grid shallow water (SW) model and explain why the two budgets do not match. We then give results for the discretized vorticity budget and discuss the implications in terms of modelling of wind driven gyres in presence of step-like coastlines.

The General Form of the Discretized Vorticity Budget

We consider the shallow water equations

  

$$\partial_{t} \mbox{\bf u}+ \mbox{\bf u}\cdot\mbox{\boldmath$\nabl... ...\nabla$ }(gh) = \frac{\mbox{\boldmath$\tau$ }}{h} + \nu \nabla^{2} \mbox{\bf u}$$ (4.1)

$$\partial_{t} h + \mbox{\boldmath$\nabla$ }\cdot(\mbox{\bf u}h) = 0$$ (4.2)

where the variables are given in Table 2.1. It is sometimes convenient to recast the non-linear terms in (4.3) in the following form:

$$\partial_{t} \mbox{\bf u}+ q \mbox{\bf k}\times (\mbox{\bf u... ...x{\boldmath$\tau$ }}{h} + \nu \nabla^{2} \mbox{\bf u}\mbox{~,}$$ (4.3)

where q and B are also given in Table 2.1. The kinematic boundary condition is no normal flow and the dynamic boundary condition is taken to be free-slip. The vorticity equation is found by taking the curl of (4.3),

$$\partial_{t} \zeta + \mbox{\boldmath$\nabla$ }\cdot (q h \mbo... ...x{\bf u}+ \frac{\mbox{\boldmath$\tau$ }}{h} \right) \mbox{~~.}$$ (4.4)

Upon integration of this equation over a closed basin, the divergence of the potential vorticity mass flux cancels out and we get

$$\partial_{t} \left( \int_{\Omega} \zeta dx dy \right) = \nu ... ...\frac{\mbox{\boldmath$\tau$ }\cdot \mbox{\bf dl}}{h} \mbox{~.}$$ (4.5)

Eqns. (4.1-4.2) or (4.2-4.3) can be discretized in different ways. To simplify the discussion, we leave the time derivative being continuous, and restrict ourselves to the C-grid. A useful general form of the SW equation is the following:

  

$$\partial_{t} u + C_u + D^-_x \Phi = \frac{\tau_{x}}{\overline{h}^{x}} + F_{x}$$ (4.6)

$$\partial_{t} v + C_v + D^-_y \Phi = \frac{\tau_{y}}{\overline{h}^{y}} + F_{y}$$ (4.7)

$$\partial_{t} h + D^+_x U + D^+_y V = 0$$ (4.8)

where $\mbox{\bf C}= (C_u,C_v)$ represents the advection-Coriolis terms, $\Phi$represents a potential function, F =(F_x,F_y) are the viscous terms and other notation is described in Section 2.2.2. The exact forms of $\mbox{\bf C}$, $\Phi$ and F depend on choices made with respect to the discretization. For example, $\Phi$ might represent the Bernoulli function or simply the pressure, depending on whether a formulation based on (4.1) or on (4.3) is employed. We first make a general point about numerical vorticity budgets and later discuss the peculiarities specific to choices for $\mbox{\bf C}$, $\Phi$ and F. From (4.6) and (4.7) we write the discretized vorticity equation

$$\partial_{t} \zeta = - D^-_x C_v + D^-_y C_u + D^-_x F_{y} ... ...-_y \left( \frac{\tau_{x}}{\overline{h}^{x}} \right) \mbox{~.}$$ (4.9)

This equation is defined at interior $\zeta$-nodes (excluding the boundary nodes), because it requires defining momentum equations at all neighboring velocity nodes (white squares in Fig.4.1). Now we want to sum over all interior $\zeta$-indices in order to get the model vorticity budget. For simplicity, we write vectors in place of x- and y- components, even though the components are not discretized at the same location (see Chapter2):

$$\partial_{t} \sum_{ij\in \Omega_{\zeta}} \zeta \Delta x \De... ...}}{\overline{h}} )\cdot\mbox{\boldmath\bf$\Delta$ l}\mbox{~,}$$ (4.10)

where $\delta\Omega_{\zeta}$ is the ensemble of indices representing the velocities nodes of the envelope of the interior vorticity node domain, $\Omega_{\zeta}$ (black nodes in Fig.4.1). We rewrite (4.10) in a more convenient form by defining

$$\mathcal{F}_{adv}= \sum_{ij\in\delta\Omega_{\zeta}} \mbox{\bf C}\cdot\mbox{\boldmath\bf$\Delta$ l}\mbox{~,}$$ (4.11)

$$\mathcal{F}_{vis}= \sum_{ij\in\delta \Omega_{\zeta}} \mbox{\bf F}\cdot\mbox{\boldmath\bf$\Delta$ l}\mbox{~,}$$ (4.12)

$$\mathcal{F}_{i}= \sum_{ij\in\delta\Omega_{\zeta}} \left(\frac... ...rline{h}}\right) \cdot \mbox{\boldmath\bf$\Delta$ l}\mbox{~~.}$$ (4.13)

Thus (4.10) becomes

$$\partial_{t} \sum_{ij\in \Omega_{\zeta}} \zeta \Delta x \Del... ...\mathcal{F}_{i}+ \mathcal{F}_{adv}+ \mathcal{F}_{vis}\mbox{~.}$$ (4.14)

$\mathcal{F}_{i}$ (flux in) is the wind input of vorticity and $\mathcal{F}_{o}$(flux out) is the sum of the viscous diffusion flux, $\mathcal{F}_{vis}$, and of the advective flux, $\mathcal{F}_{adv}$. The important point here is to note that $\mathcal{F}_{adv}$ideally should be zero since it represents an advective flux through the basin lateral boundary. It is not zero in the numerical model because the domain boundary for the model vorticity budget is located half a grid point inside the domain (see Fig. 4.1). However, as the resolution increases, the region delimiting the vorticity budget domain approaches the model boundary and, $\mathcal{F}_{adv}$ should converge to zero. How quickly this occurs will depend on the numerical formulation.

It is always possible to approximate the vorticity budget at the model boundary by using off-centered derivatives and interpolating some of the variables to the boundary. The model numerics, however, make no use of variable values found by such an interpolation and therefore, a vorticity budget calculated in this way must be considered distinct from the model vorticity budget. Such a budget might misrepresent the contribution of the different terms of the discretized equations of the model, especially if the error introduced by the coastline discretization is of lower order than are the truncation errors of the model. For this reason, we prefer to use the model vorticity budget. We note also that the truncation errors in the model vorticity budget are larger than the truncation errors in solving the shallow water equations, since vorticity is a higher order variable.

Figure 4.2 compares the rotated and non-rotated basin cases. The integrand ( $\mbox{\bf C}\cdot \mbox{\boldmath\bf$\Delta$ l}$--i.e. the local $\mathcal{F}_{adv}$) is plotted as a function of distance around the basin perimeter and the position of the steps is evident from the abrupt jumps in the integrand value. When summed along the perimeter, $\mathcal{F}_{adv}$ is non-zero and is larger for the rotated basin experiment compared to the non-rotated basin experiment. Starting from this observation, we are interested in quantifying the importance of the steps over the global vorticity budget. First, as increased resolution leads to more steps, and due to the singular behavior of $\mbox{\bf C}\cdot \mbox{\boldmath\bf$\Delta$ l}$ close to steps, it is no longer obvious that $\mathcal{F}_{adv}$ converges to zero with increasing resolution. From this point of view, $\mathcal{F}_{adv}$ is probably very sensitive to the formulation of the advective terms in (4.10), and extension in (4.1,4.3). Second, we want to investigate whether the overall circulation is sensitive to the presence of an extra term in the global vorticity budget, as $\mathcal{F}_{adv}$ can be a source or a sink term, depending on its sign.

   
Numerical Formulations

We are interested in applying different formulations for the advection-Coriolis terms since we noted that $\mathcal{F}_{adv}$ was generally non-zero for the single gyre Munk problem, with the integrand being particularly large at steps. The two advective numerical schemes that we consider are the conventional formulation (based on Eq. 4.1) and the potential enstrophy conserving formulation (based on Eq. 4.3).

In addition to testing for sensitivity to the choice of advective schemes, we also consider different formulations of the stress tensor. That the overall circulation is primarily sensitive to the formulation of the stress tensor is the main result of AM, who found that the $\delta$-$\zeta$ formulation gave better results than the conventional formulation. We refer to [Gent, 1993] and [Shchepetkin and O'Brien, 1996] for a more complete discussion on appropriate viscous stress tensor formulations for the shallow water equations and we limit ourselves to the two stress tensor formulations used by AM. Below, we review these two formulations. Thus, we are interested in testing four combinations of two advective and two diffusive formulations. Table 4.1 summarizes these four different combinations.

With respect to the advection-Coriolis terms, we compare the conventional formulation to the potential enstrophy conserving formulation of [Sadourny, 1975] . For the conventional formulation, $\Cor$ and $\Phi$ are given by

$$u D^o_x u + \overline{\overline{v}^{x}}^{y} D^o_y u - f \overline{\overline{v}^{x}}^{y}$$ 4.15

and for the potential enstrophy conserving formulation, $\Cor$ and $\Phi$ are given by

Both formulations ensure a second order accuracy to the discretized SW equations. For the conventional formulation, changes are made to incorporate the boundary conditions at second order of accuracy, by using off-centered differencings close to the boundary. No boundary condition for the vorticity is required. However, since the enstrophy conserving scheme explicitly uses the vorticity, this formulation requires that vorticity be specified at boundary points. We choose to set the relative vorticity to zero along the model boundary, which is consistent with the free-slip boundary condition along straight walls. Also, contrary to the conventional formulation, no off-centered differencing is needed at the boundary for the computation of $\Cor$.

The two numerical formulations that we consider for the viscous terms are the divergence-vorticity tensor formulation of [Madec et al, 1991] and the conventional five-point Laplacian. For the latter,

As with the conventional advection formulation, changes are made here to incorporate the boundary conditions at second order accuracy, by using off-centered differencings. Another technical remark concerns the treatment of velocity points close to tips of land. For those points (the $u$ and $v$ points of Fig.4.3, the tip of the land is half a grid cell away. Let us focus on the $u$-point. The problem is to evaluate the Laplacian of $u$ at this point. A five-point Laplacian requires knowledge of $\partial u/\partial y$ in the center of the northern and southern sides of the cell surrounding the $u$-point, and $\partial u/\partial x$ on the eastern and western sides. The problem lies with $\partial u/\partial y$ on the northern side. The usual treatment would have

$$du/dy {north} =(u_{i,j+1}-u_{ij)/Delta y$$ 4.18

which simplifies to

because of the impermeability condition which sets $u_{i,j+1}$ to zero. Alternatively, one might take impermeability to imply that the tip is a stagnation point, in which case an off-centered differencing leads to

A third logical possibility would be to apply the free slip condition at the tip to conclude that

We choose the latter (4.21), in order to let the ``fluid'' slip as much as possible along the walls since the first two conditions (4.19,4.20) tend to slow down the boundary currents. A more accurate formulation of the boundary condition close to the steps can be derived using a finite volume formulation, which treats the northern viscous flux as a mean between (4.19) and (4.21). However, this would slow down the boundary current due to the use of (4.19). In addition, more accurate treatment of the steps have limited value as the steps are artificial.

The divergence-vorticity ( $\delta$-$\zeta$ ) form of the stress tensor leads to the following form for the Laplacians

where $\delta$ is the divergence expressed at the $h$-location (center of the cell). This formulation is more general in the sense that there is no adjustment of the formulation at the boundary. Another remark concerns the case of straight walls. In that particular case, there is no difference between the $\delta$-$\zeta$ stress tensor formulation and the traditional formulation. The difference is in the treatment of steps.

To illustrate this, we consider a comparison between the two stress tensor formulations for a forward step along a north-flowing western boundary current. Choose $(i,j)$ so that, in Figure~\ref{cgrid3}, the $\zeta$-point right at the tip of the land corner would have $(i,j+1)$ indices (see Fig.~A.1 for indices arrangement). Thus, the viscous terms under the $\delta$-$\zeta$ formulation are

where the additional terms are positive. These additional terms represent a forward-acceleration. As AM noted, a serious inconvenience of the conventional formulation is that, in presence of steps, there is ``extra diffusion'' of momentum due to additional velocity points set to zero at the boundary (the impermeability condition), as compared to the straight wall case. This extra diffusion is responsible for slowing down the boundary currents. Therefore, the accelerating terms of the $\delta$-$\zeta$ formulation partly compensate the decelerating terms of the conventional formulation.

A final remark is that the divergence part of the viscous forces cancels out in the vorticity equation. Therefore, in the discretized vorticity equation, the $\delta$-$\zeta$ formulation leads to a viscous term that takes the form of the five-point Laplacian of the vorticity. This is not true of the conventional formulation.

Results

 

Table 4.1: The four combinations of advection formulations and stress tensor formulations.

A	B	C	D
enstrophy preserving advection	enstrophy preserving advection	conventional advection	conventional advection
conventional stress tensor formulation	$\delta$-$\zeta$	conventional stress tensor formulation	$\delta$-$\zeta$

By studying $\mathcal{F}_{adv}$, we want to address several issues related to the accuracy of the different combinations of the advection and diffusion formulations and their influence on the strength of the overall circulation. Firstly, a major requirement is that, whatever the geometry of the basin, $\mathcal{F}_{adv}$ should converge to zero as the resolution goes to infinity. This test allows us to rank the performances of the model for the different combinations of advective and diffusive schemes. Of particular interest will be the importance of the advective formulation. A second concern is to assess whether the size of the artificial source or sink of vorticity due to $\mathcal{F}_{adv}$ influences the overall strength of the gyres. A third concern relates to the general accuracy of model vorticity budgets.

To address these issues, we make use of the conceptual experiment proposed by AM, in which a single gyre Munk circulation is computed in rotated and non-rotated square basins. In both cases, all parameters and forcing are unchanged except for the discretized coastline. The four combinations (A, B, C, D) of numerical formulations we propose to test are detailed in Table 4.1. One remark concerns the non-rotated basin results. There, since the conventional and $\delta$-$\zeta$ stress tensor formulations are identical, the results for the B combination are identical to the results for A. The same applies for the C and D cases.

We reproduce the results of AM in Figure 4.4. This figure shows the elevation fields for the A and B cases and for no rotation and a small rotation angle of 3.4$^o$. Clearly, the A case shows circulation patterns collapsing as the number of steps along the walls increases whereas, for the B case, the circulation is quite similar to the original non-rotated circulation. The results for C are not shown but are very similar to the results for A. The results for D show a small increase in the strength of the gyre compared to A, but the original overall circulation of A-B with no rotation is not recovered (not shown).

Figure 4.5a shows the kinetic energy as a function of resolution for the various combinations and for a rotation angle of $3.4^o$. Only the B combination converges to non-rotated solutions. The A and C results are almost identical, but appear to converge to a kinetic energy that is reduced by over a factor of 2 compared to the non-rotated cases. For the D combination, kinetic energy decreases and then tends to slightly increase with increasing resolution and is generally much lower than for A-B with no rotation or B with rotation.

As mentioned, the first consistency test related to the vorticity budget is to verify that $\mathcal{F}_{adv}$ converges to zero with increasing resolution. For all rotation angles considered and for the B combination, this statement appears to be true. For the other combinations (A, C, D), such is not the case, at least for certain angles. For instance, $\mathcal{F}_{adv}$ tends to increase or stay constant for the A, C and D combinations at $3.4^o$ (Fig. 4.5b). For the D case, $\mathcal{F}_{adv}$ increases dramatically with increasing resolution ---so much that $\mathcal{F}_{adv}$ becomes larger than the wind input. Associated with this is a reverse (negative) viscous flux. This behavior may have consequences on the stability of the model. Although no obvious numerical instabilities occurred for a rotation angle of $3.4^o$, numerical instabilities cause the model to crash for other angles, for example at $-30^o$. It seems plausible that this behaviour is associated with the large (and opposing) advective and diffusive fluxes of vorticity near the model perimeter. In any event, it seems reasonable to conclude that the D combination is inappropriate. This implies that the $\delta$-$\zeta$ viscous formulation performs well only when used is conjunction with the enstrophy conserving advection. This finding complements that of AM. For the A and C combinations (Fig. 4.5b), $\mathcal{F}_{adv}$ does not converge toward zero with increasing resolution. Hence, these two combinations seem inappropriate, even if the resulting solutions are always stable.

We now address the issue of possible correlation between $\mathcal{F}_{adv}$ and the kinetic energy. Given that inertial runaway (the inability of simple models of the ocean to converge to a reasonable statistical mean solution as the eddy viscosity is decreased to the real value of the viscosity found in water) appears to be related to ``difficulties'' in balancing the global vorticity budget [Pedlosky, 1996] , it seems reasonable to ask whether the sign of $\mathcal{F}_{adv}$ is correlated with an indicator of the overall strength of the gyre, such as total kinetic energy. For example, when $\mathcal{F}_{adv}$ is negative, it adds to the wind input of vorticity and one might expect a stronger gyre to result. Some evidence that this may be the case is found by comparing B and D, which share the same formulation of the viscous terms. Figure 4.5b shows that $\mathcal{F}_{adv}$ is positive and larger for D than is the case for combination B. Thus the total wind plus advective input of vorticity to the basin is stronger in case B. As might have been anticipated, B shows a more energetic circulation (Fig.~\ref{ene4}a). It is also interesting to see whether there is any correlation between kinetic energy and the sign/strength of $\mathcal{F}_{adv}$ for a given formulation of the numerics. We restrict this discussion to the use of the B combination. From figures 4.6 and 4.7, which show the kinetic energy and advective/wind vorticity input ratio for a range of resolution and rotation angles, there does not appear to be any striking correlation. In particular, if we focus on the region of negative values of $\mathcal{F}_{adv}$ (i.e., for a case where $\mathcal{F}_{adv}$ has the same sign as the wind input), the kinetic energy for this region is not larger than the kinetic energy at the same resolution but for an opposite angle (in fact, the kinetic is slightly lower). Presumably the added advective flux in this region is locally balanced by the viscous terms, so that processes analogous to those thought to be responsible for inertial runaway do not lead to an increase in the overall strength of the gyre.

To conclude this section, we investigate the general accuracy of model vorticity budgets with respect to $\mathcal{F}_{adv}$ using the B combination, only, since this combination is the only one showing a robust convergence to zero with increasing resolution. As $\mathcal{F}_{adv}$ should ideally not be present in the vorticity budget, the viscous flux, $\fov$, can be either underestimated or overestimated (which modifies the local balance at the wall and therefore the strength of the gyre) and $\mathcal{F}_{adv}$ can be viewed as an error. From Figure 4.8 and for the range of resolution we used, $\mathcal{F}_{adv}$ varies between 5\% and 50\% of the wind input. The order of the convergence for $\mathcal{F}_{adv}$ with increasing resolution is fairly close to unity or slightly lower for all positive angles. For negative angles, we did not compute the convergence order because $\mathcal{F}_{adv}$ goes through a minimum (Figure 4.7) and had not asymptoted to an uniform convergence order at the highest resolutions we considered. A noteworthy point is that the effect of increasing the rotation angle (introducing more steps) seems to decrease the convergence order (1/2 at $20^{o}$). Paradoxically, however, the convergence order increases again to reach unity for $45^{o}$, the rotation angle at which the number of steps is maximum. In fact, at this angle $\mathcal{F}_{adv}$ even shows a negative offset compared to the $0^{o}$ angle.

Except for effects related to step-like boundaries, that the convergence order is unity follows directly from the order of discretization of the vorticity. Since the vorticity is one order higher a variable than is velocity, and since the velocity is computed at second order accuracy, it follows that the vorticity is at best accurate to first order. Therefore, $\mathcal{F}_{adv}$ can be considered an explicit first order (at best) error in the vorticity budget. For the B combination, we observe that the convergence order for $\mathcal{F}_{adv}$ varies between 1/2 and unity, depending on the rotation angle. In the 1/2 order case, errors (or discrepancies) vary between 25% (high resolution) and 50% (coarse resolution) and, in the first order case, they vary between 6% and 22%. The errors are much larger for the other combinations and can reach 100%. This implies that the accuracy of computing vorticity budgets from primitive equations models is fairly low, especially in absence of attention to the numerics. These errors may also vary a lot considerably with the discretized domain geometry.

Get the PS or PDF version here