# Low-order pressure gradient schemes in sigma coordinate models: The seamount test revisited

Patrick Marsaleix, Francis Auclair, Claude Estournel

► To cite this version:

Patrick Marsaleix, Francis Auclair, Claude Estournel. Low-order pressure gradient schemes in sigma coordinate models: The seamount test revisited. Ocean Modelling, 2009, 30, pp.169-177. 10.1016/j.ocemod.2009.06.011. hal-00518947

HAL Id: hal-00518947

https://hal.science/hal-00518947v1

Submitted on 5 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Low-order pressure gradient schemes in sigma coordinate models: The seamount test revisited

Patrick Marsaleix *, Francis Auclair, Claude Estournel

Laboratoire d'Aérologie, CNRS et Université de Toulouse, 14, Avenue Edouard Belin, 31400 Toulouse, France

This paper revisits the classic seamount test used in numerous previous studies to evidence the sigma errors of the pressure gradient force (PGF) and their long-term effects on circulation. Two kinds of analysis are developed. We first consider the initial PGF errors. Then, the global level of erroneous kinetic energy is computed along a 180-day simulation. The long-term circulation appears to be better correlated to the initial vorticity errors than to the initial error diagnostics.

The original feature of this study is to reconsider the currently admitted idea that *Density-Jacobian* type PGFs perform better than the primitive sigma formulation discretized in a straightforward way (hereafter *Straightforward-Primitive* PGF). Errors on the discrete hydrostatic pressure are actually closely related to the way the density field is initialized. If a mass conserving method is preferred to a straightforward initialization, the *rectangular* integral of the *Straightforward-Primitive* PGF is likely to be more accurate than the *trapezoidal* rule usually involved in *Density-Jacobian* PGFs. Errors on the vorticity field of the *Straightforward-Primitive* PGF depend on the discretization of the hydrostatic correction term. A modified version of the *Straightforward-Primitive* PGF is shown to be in better agreement with the concept of *bottom torque consistency*. The seamount tests show that this so-called *Modified-Primitive* PGF performs globally better than the current low-order *Density-Jacobian* PGFs.

## 1. Introduction

Truncation error on the pressure gradient force (PGF) is a well-known and challenging problem in ocean models using topography-following coordinates. It has given rise to a vast literature (see for instance references in Shchepetkin and McWilliams, 2003, hereafter SMW03), from which the “hydrostatic inconsistency” emerged as a key concept. The “hydrostatic inconsistency” can be viewed as a grid level slope threshold beyond which increasing vertical resolution destroys PGF accuracy (Haney, 1991). The “hydrostatic inconsistency” is a useful but incomplete indicator of the grid propensity to develop erroneous circulation, since initial PGF errors also depend on the initial density field. This led several authors to propose more complete diagnostics, both considering the grid geometry and the vertical gradient of the background density (Mellor et al., 1994). As a matter of fact, pressure gradient errors can be efficiently reduced by methods adapting the coordinate stretching function of sigma coordinate models to the initial background density field, as shown by Ciappa (2008). Alternatively, for a given grid geometry, the background density field can be optimized in order to minimize PGF errors (Auclair

et al., 2000). A purely vertical reference state, possibly deduced from a horizontal averaging of the density field, can also be removed before computing the PGF (Haney, 1991).

However, a priori estimates of the erroneous circulation, based on the assumption of a geostrophic balance with PGF errors, are somehow limited since they do not reflect the possible growth of the error level during the time of the simulation. As far as C-grid, free-surface, terrain-following type models are concerned, SMW03 and Mellor et al. (1998) (hereafter MOE98) have described the long-term effect of PGF errors on barotropic and baroclinic circulations.

Along with the development of the POM model in the eighties (Blumberg and Mellor, 1987) and its wide spreading in the coastal modelling community, the discrete PGF formulation of the POM model has generally been considered as a better alternative to the primitive sigma coordinate formulation which is given by:

$$-\frac{1}{\rho_0} \frac{\partial p}{\partial x}|_{z} = -g \frac{\partial \zeta}{\partial x} - \frac{1}{\rho_0} \left( \frac{\partial p'}{\partial x} \right)_{s} + g \rho' \frac{\partial z}{\partial x} \quad (1)$$

where *p* is the pressure, *ζ* the free-surface elevation, *ρ*₀ a reference density, *ρ'* = *ρ* − *ρ*₀ with *ρ* the density. Subscript *s* refers to a differentiation along a constant level in the topography-following coordinate system and subscript *z* simply refers to the truly horizontal differentiation. The hydrostatic pressure anomaly *p'* is given by:

<footnote>* Corresponding author. Tel.: +33 5 61 33 27 63; fax: +33 5 61 33 27 90. E-mail address: marp@aero.obs-mip.fr (P. Marsaleix).</footnote>
$$p'(x, z) = g \int_z^\zeta \rho' dz' \qquad (2)$$

In the following, the straightforward discrete form of (1) and (2), presented in Section 2.1, is referred to as the *Straightforward-Primitive* PGF.

On the other hand, the PGF used in the POM model is a discrete form derived from:

$$-\frac{1}{\rho_0} \frac{\partial p}{\partial x} \Big|_z = -g \frac{\rho(\zeta)}{\rho_0} \frac{\partial \zeta}{\partial x} - \frac{g}{\rho_0} \int_z^\zeta \left( \frac{\partial \rho}{\partial x} \Big|_s - \frac{\partial \rho}{\partial z} \frac{\partial z}{\partial x} \Big|_s \right) dz' \quad (3)$$

Because of the terms in brackets at the right-hand side (hereafter RHS) of (3) and the symmetries of its possible discretization, the POM PGF is often regarded as a *Density-Jacobian* method (SMW03). The POM PGF was first expressed using the conventional sigma coordinate system (Blumberg and Mellor, 1987). A formulation suitable for more generalized vertical coordinate systems has been proposed by Song (1998). Following Song's terminology, the standard discrete form of (3) (detailed in SMW03 and Song (1998)) is hereafter referred as to the *Standard-Jacobian* PGF.

It has been argued that the *Standard-Jacobian* PGF is more accurate than the straightforward discrete form of (1) (SMW03) since in the case of a density field linear in z and constant in x, the discrete *Standard-Jacobian* exactly vanishes (as expected) when the *Straightforward-Primitive* PGF does not. A new scheme derived from the *Standard-Jacobian* method has been suggested by Song (1998), the so-called *Weighted-Jacobian* providing an exact form of the discrete PGF when density field is both linear in z and x.

The seamount test, presented by Beckmann and Haidvogel (1993) and then revisited by several authors (Chu and Fan, 1997; MOE98, Song, 1998, SMW03, Ciappa, 2008) has enabled to test and compare different numerical forms of the PGF, for various types of grids and models. Basically, this numerical experiment consists of academical simulations of the ocean circulation around a gaussian-shaped seamount, using an initial density field which is non-linear in z and constant in x and y. Because of the vertical non-linearity of $\rho$, discrete PGFs do not exactly vanish and lead to long-term erroneous circulations. The recent seamount study of SMW03 presented a comparison of various, low- and high-order, *Density-Jacobian* type, PGF schemes. As far as low-order methods are concerned, the SMW03's study showed that the *Standard-Jacobian* performed better than the *Weighted-Jacobian*, the schemes being mainly evaluated according to the global level of erroneous kinetic energy (hereafter ERKE) during the simulation. SMW03 also showed that the *Blended-Jacobian* PGF, namely a combination of the *Standard-Jacobian* and *Weighted-Jacobian* schemes, performed better than these two schemes separately. Although a fifty-fifty blend seemed to provide the best performance, SMW03's study did not really permit to conclude whether this particular combination would still be optimal in other circumstances. Not surprisingly, higher order methods reduce the ERKE (SMW03, Chu and Fan, 1997) but it is not clear whether these sophisticated schemes are really compatible with the principle of energy conservation involving the discrete form of the density advection, notably because high-order pressure schemes are not consistent with the usual low-order definition of the potential energy (SMW03, p35-30).

The original point of the present study is to bring new insight on the *Straightforward-Primitive* PGF, generally ignored in recent studies on the grounds that this method is a priori less accurate than *Density-Jacobian* type methods. First of all we will see in Section 2 that the inherent superiority of the *Density-Jacobian* method over the *Straightforward-Primitive* method is not so obvious, notably because the answer to this question largely depends on the way the discrete density field is defined. Besides, we will see that the initialization method has possible repercussions on the writing of

the discrete hydrostatic correction term of the *Straightforward-Primitive* PGF. The so-called *Modified-Primitive* PGF possibly leads to reconsider the discretization of the density advection scheme. Secondly, SMW03 emphasised that initial PGF errors should not be regarded as being representative of errors developed during the simulations, and that any PGF scheme should in particular be evaluated based on long-term simulations. On this particular point too, the proof of the inaccuracy of the *Straightforward-Primitive* PGF has not been provided by the previous studies. We reproduced in Section 3 the long-term seamount simulations and compared the *Straightforward-Primitive* and *Modified-Primitive* methods to current low-order *Density-Jacobian* schemes, namely the *Standard-Jacobian*, the *Weighted-Jacobian* and the *Blended-Jacobian* PGFs. Following MOE98 who identified a sigma error of the second kind, the growth of the errors is examined in the light of the discrete vorticity of the vertical integral of the PGF (Section 4). We will see that the vorticity analysis provides another good reason to reconsider the writing of the discrete hydrostatic correction term of the *Straightforward-Primitive* PGF (Appendix).

## 2. Initialization issues

### 2.1. *Straightforward-Primitive* PGF: The discrete formulation

The straightforward discrete form of the primitive formulation of the PGF, corresponding to the terms in brackets at the RHS of (1), is given by:

$$\left. \frac{\partial p'}{\partial x} \right|_s + g \rho' \frac{\partial z}{\partial x} \Big|_s = \frac{p'_{i,k} - p'_{i-1,k}}{\Delta x} + g \frac{\rho'_{i,k} + \rho'_{i-1,k}}{2} \frac{z_{i,k} - z_{i-1,k}}{\Delta x} \quad (4)$$

where $(p'_{i,k}, \rho'_{i,k}, z_{i,k})$ are the discrete values of, respectively, pressure anomaly, density anomaly and depth, in the middle of the grid box corresponding to horizontal index $i$ and vertical index $k$. As in Marsaleix et al. (2008), we first defined $(z_{i,k+1/2}, z_{i,k-1/2})$, the depth at the top and bottom of grid boxes, and then computed $z_{i,k}$ according to

$$z_{i,k} = (z_{i,k-1/2} + z_{i,k+1/2})/2 \quad (5)$$

The thickness of grid box $(i,k)$ is given by:

$$\Delta z_{i,k} = z_{i,k+1/2} - z_{i,k-1/2} \quad (6)$$

The discrete form of the hydrostatic pressure is:

$$p'_{i,k} = p'_{i,k+1/2} + 0.5g \rho'_{i,k} \Delta z_{i,k} \quad (7)$$

where $p'_{i,k+1/2}$, the pressure at the top of grid box $(i,k)$, is given by a rectangular vertical integration of the density profile:

$$p'_{i,k+1/2} = g \sum_{q=k+1,k_{\max}} \rho'_{i,q} \Delta z_{i,q} \quad (8)$$

where $k_{\max}$ is the vertical index of the topmost grid box.

We can note that an equivalent expression for (7) is:

$$p'_{i,k} = \frac{p'_{i,k+1/2} + p'_{i,k-1/2}}{2} \quad (9)$$

The location of the numerical variables on the C-grid is shown on Fig. 1.

### 2.2. "Straightforward" or "Volume-Averaged" initialization

It can be shown (SMW03) that the discrete *Standard-Jacobian* would be equivalent to the *Straightforward-Primitive* scheme if the rectangular integration (8) were simply replaced by a trapezoidal scheme. Consequently, we can think that the difference in accuracy is related to the method used to compute the hydrostatic pressure and thus, it is tempting to conclude that the trapezoidal
<figure><img src="image_2.png" /><figcaption>Fig. 1. Location of variables on the vertical C-grid.</figcaption></figure>

method of the Standard-Jacobian is a priori more accurate than the rectangular integration of the Straightforward-Primitive scheme. However, this reasoning is somehow limited, notably because PGF errors also depend on how the discrete density field has been initialized. Let us thus consider the initial density field, exponential in z, constant in x and y, used by SMW03 for the seamount test. Discrete pressure can not be exact, whatever the integration method, rectangular or trapezoidal, if the non-linear analytical density has been discretized in a straightforward way, i.e.

$$ \rho'_{i,k} = \alpha e^{\frac{z_i}{\delta}} \qquad (10) $$

On the other hand, defining the discrete density in such a simple way is somehow questionable since it is generally admitted that a discrete variable should represent some average of the true field over the volume of the considered cell box (Ferziger and Peric, 2002). Practically, it means that (10) could be replaced by:

$$ \rho'_{i,k} = \frac{\int_{z_{i,k+1/2}}^{z_{i,k+1/2}} \alpha e^{\frac{z_i}{\delta}} dz'}{z_{i,k+1/2} - z_{i,k-1/2}} = \alpha \delta \frac{e^{\frac{z_{i,k+1/2}}{\delta}} - e^{\frac{z_{i,k-1/2}}{\delta}}}{z_{i,k+1/2} - z_{i,k-1/2}} \qquad (11) $$

In the following, (11) will be referred to as a Volume-Averaged initialization. If initial density is given by (11), pressure $p'_{i,k+1/2}$ computed by the rectangular method is exact, since substituting (11) into (8) (and for simplicity assuming $\zeta = 0$) leads to $p'_{i,k+1/2} = g\alpha\delta(1 - e^{\frac{z_{i,k+1/2}}{\delta}})$. We note, however, that the mid-box pressure obtained from (9), i.e. the average of the pressure at the top and bottom of grid box (i, k), is eventually not exact since the analytical pressure between intermediate levels $z_{i,k-1/2}$ and $z_{i,k+1/2}$ is not linear. In return, the fact that the pressure at intermediate levels is exact, makes the pressure at the centre of the grid box potentially more accurate than the same pressure computed with a trapezoidal integration of the straightforward discrete density (9). This becomes particularly true at deepest levels since the pressure vertical profile tends to be linear when z decreases. As a consequence, we will see in the following that the Primitive method is possibly more accurate than Density-Jacobian type PGFs, provided that a suitable discrete form of the hydrostatic correction term, $g\rho'\frac{\partial z}{\partial x}$, accounting for the Volume-Averaged initialization concept, is also adopted.

Besides, we note that the discrete integration of (11) over the whole water column, $\sum_{k=1,kmax} \rho'_{i,k} \Delta z_{i,k}$, is the same as its analytical counterpart, $\alpha\delta(1 - e^{\frac{z}{\delta}})$. In other words, the initialization method (11) has also the advantage to be more conservative (vis-à-vis to the global density) than the straightforward method (10).

We also note that most of C-grid ocean models use a conservative form of advection and diffusion processes based on discrete flux differences, ensuring that the first order volume integral $\sum_V \rho'_{i,j,k} \Delta V_{i,j,k}$ remains unchanged, provided that advection and diffusion fluxes at the boundary of V are globally zero. This implicitly supposes that $\rho'_{i,j,k}$ represents an average value over the cell box (i, j, k), and thus justifies to choose (11) rather than (10).

## 2.3. Hydrostatic cancellation

A situation of hydrostatic inconsistency is reached when (SMW03):

$$ r_x = \left| \frac{z_{i,k+1/2} + z_{i,k-1/2} - z_{i-1,k+1/2} - z_{i-1,k-1/2}}{z_{i,k+1/2} - z_{i,k-1/2} + z_{i-1,k+1/2} - z_{i-1,k-1/2}} \right| > 1 \qquad (12) $$

An example of the case $r_x = 1$ is given on Fig. 1. Pioneer studies (see for instance Haney, 1991) showed that the PGF accuracy is improved with higher vertical resolution, as long as the grid remains hydrostatically consistent. As a matter of fact, it is commonly admitted that model grids should strive to avoid hydrostatic inconsistency, or at least should not lead to unreasonably high values of $r_x$. For a given number of vertical levels, $r_x$ can be reduced by increasing horizontal resolution or smoothing the discrete bathymetry. However, within the range of variation $0 \le r_x \le 1$, the highest PGF errors do not necessarily correspond to $r_x = 1$ (Mellor et al., 1994; MOE98). Actually, PGF errors can exactly cancel in this latter case (SMW03). The so-called hydrostatic cancellation property of the PGF, can be recovered in the case of the primitive formulation of the PGF and the Volume-Averaged initialization, provided that the discrete expression (4) is rewritten as follows:

$$ \frac{\partial p'}{\partial x}|_s + g \rho' \frac{\partial z}{\partial x}|_s = \frac{p'_{i,k} - p'_{i-1,k}}{\Delta x} + g \frac{\Delta z_{i,k} \rho_{i,k} + \Delta z_{i-1,k} \rho_{i-1,k} z_{i,k} - z_{i-1,k}}{\Delta z_{i,k} + \Delta z_{i-1,k}} \frac{\Delta x}{\Delta x} \qquad (13) $$

Indeed, let us consider the example of the hydrostatic inconsistency situation shown on Fig. 1. We have $z_{i,k-1/2} = z_{i-1,k+1/2}$ and (since the Volume-Averaged initialization method leads to exact pressure at intermediate levels) $p'_{i,k-1/2} = p'_{i-1,k+1/2}$. Using 5, 6, 9, it follows that $z_{i,k} - z_{i-1,k} = (z_{i,k+1/2} - z_{i-1,k-1/2})/2$, $\Delta z_{i,k} + \Delta z_{i-1,k} = z_{i,k+1/2} - z_{i-1,k-1/2}$ and $p'_{i,k} - p'_{i-1,k} = (p'_{i,k+1/2} - p'_{i-1,k-1/2})/2$. Then the RHS of (13) becomes equivalent to $(p'_{i,k+1/2} - p'_{i-1,k-1/2} + g\Delta z_{i,k}\rho_{i,k} + g\Delta z_{i-1,k}\rho_{i-1,k})/(2\Delta x)$. Using (8), the latter is finally equal to $(p'_{i,k-1/2} - p'_{i-1,k+1/2})/(2\Delta x)$, i.e. zero. In the following, the Straightforward-Primitive terminology will be restricted to (4) and, on the other hand, scheme (13) will be referred to as the Modified-Primitive PGF.

## 2.4. Energy conservation

The principle of global energy conservation requires to discretize PGF and density advection in a consistent way. Practically, the energetic consistency has mainly two consequences. On one hand, the discrete hydrostatic pressure (7)–(9) must be consistent with vertical density advection and on the other hand, the discretization of the hydrostatic correction term, $g\rho'\frac{\partial z}{\partial x}$, appearing in the RHS of (1), has to be consistent with horizontal density advection (SMW03). As far as the Straightforward-Primitive PGF is concerned, a possible energy-conserving form of the horizontal mass flux $u\rho'$, located at the velocity grid point, is based on a simple half-half average of the density field (Marsaleix et al., 2008), when the Modified-Primitive PGF leads to compute $u\rho'$ according to

$$ u_{i-1/2,k} \frac{\Delta z_{i-1,k} \rho_{i-1,k} + \Delta z_{i,k} \rho_{i,k}}{\Delta z_{i-1,k} + \Delta z_{i,k}} \qquad (14) $$
## 3. The seamount test experiment

### 3.1. Description

The numerical domain is the same as the one described in SMW03. We used 11 vertical grid boxes. SMW03 used a generalized coordinate transformation leading to a horizontally variable vertical stretching function. Here, the vertical stretching is horizontally uniform and corresponds to the vertical stretching function of SMW03 when $h = h_{max}$. Practically, the depth of the intermediate vertical levels is computed according to $z/(h + \zeta) = (\sigma h_{min} + C(\sigma)(h_{max} - h_{min}))/h_{max}$, where, as in SMW03, $C(\sigma) = \sinh(0\sigma)/\sinh(0)$, $h_{max} = 5000$ m, $h_{min} = 500m$, $\theta = 3$, and $-1 \le \sigma \le 0$ is discretized uniformly. The depth of middle grid box levels is then computed according to (5). Bathymetry (Fig. 2) is given by:

$$h_{i,j} = h_{max} - h_{mount} \exp \left\{ - \frac{(i - i_{max}/2)^2 \Delta x^2 + (j - j_{max}/2)^2 \Delta y^2}{L^2} \right\} \quad (15)$$

with $h_{mount} = 4500m$, $L = 40 \times 10^3$ m, $\Delta x = \Delta y = 6700$ m, and $i_{max} \times j_{max}$, the number of grid boxes in the horizontal directions, is $48 \times 48$. The analytical initial density is $\rho = -3e^{2/500}$. Lateral boundaries are closed, so that boundary conditions leave the mechanical energy globally unchanged.

We used SYMPHONIE, the C-grid, free-surface, sigma coordinate model detailed in Marsaleix et al. (2008). As far as density is concerned, a linear equation of state is used and advection is computed with second-order centred schemes. Mixing is computed by the mean of a classic Laplacian type operator. In the present study, horizontal diffusivity is set to zero and vertical diffusivity is given by the turbulence closure of Gaspar et al. (1990). As far as velocities are concerned, vertical advection is centred and horizontal advection is upwind. Although horizontal viscosity is set to zero we note that the upwind horizontal advection scheme is inherently diffusive and somehow equivalent to set $O_x$ and $O_y$ viscosities to $(\frac{|u|\Delta x}{2}, \frac{|v|\Delta y}{2})$ (James, 1996). Vertical viscosity is given by the turbulence closure of Gaspar et al. (1990).

### 3.2. Initial PGF

Fig. 3 shows initial pressure gradient errors for different options (summarized by Table 1) of discrete PGF scheme and density initialization. The location of these profiles, indicated on Fig. 2, has been chosen in order to cover a wide range of topographic situations, from the deep area with gentle slope (Fig. 3a), to the seamount top area (Fig. 3d) characterized by strong hydrostatic inconsistency (SMW03). We will see in the following that initial

PGF errors have a rather small influence on the long-term circulations, which limits their interest. Keeping this observation in mind, Fig. 3 leads to the following comments:

1. The Modified-Primitive PGF combined to the Volume-Averaged initialization (case 1) leads to the smallest initial errors in most cases (Fig. 3a-c).

2. At the basis of the seamount (Fig. 3a) the bottom slope becomes small, leading to a small amplitude hydrostatic correction term, $g\rho'\frac{\partial\zeta}{\partial x}$, The Straightforward-Primitive PGF leads to small errors (comparable to those obtained in case 1) if the Volume-Averaged initialization is used (case 2B), and on the opposite, gives the largest errors if a Straightforward initialization is used (case 2).

3. Near the top of the seamount (Fig. 3d), the hydrostatic inconsistency (12) is especially high (SMW03), leading to the strongest PGF errors. The Standard-Jacobian PGF seems to be the best option. PGF errors are insensitive to the initialization method (case 1 is similar to case 1B, case 2 is similar to case 2B).

### 3.3. Long-term simulations

In this type of study, it has often been considered that the accuracy of a PGF scheme can not be judged on the initial errors only, but also on the results of long-term simulations (Mellor et al., 1994; MOE98). Indeed, it has been argued that the long-term erroneous circulation can be quite different from the geostrophic current simply deduced from the initial PGF (SMW03). As underlined by SMW03, the error level may grow along the simulation because not only the PGF scheme but all terms containing horizontal derivatives are in fact potentially altered by the terrain-following coordinate system. Thus, a long-term simulation enables to better estimate the relevance of a numerical system or another. As in SMW03, our simulations are 180-day long and the time evolution of the volume averaged global erroneous kinetic energy, ERKE = $\frac{1}{2V} \int (u^2 + v^2)dV$ (V being the total volume of the domain), is used to assess the different PGF schemes.

Results are summarized by Fig. 4. First, we note that $\gamma(0 \le \gamma \le 1)$ refers to the ratio of the two basic schemes used to build the Blended-Jacobian PGF ($\gamma = 0$ corresponds to the Standard-Jacobian and $\gamma = 1$ to the Weighted-Jacobian). Fig. 4a corresponds to the reference case described in Section 3.1. We also considered several other configurations derived from this reference case, by simply changing the vertical resolution, or the initial stratification or the shape of the seamount. Practically, Fig. 4b was obtained with 22 vertical levels, Fig. 4c with $\rho = -3e^{2/250}$ and Fig. 4d with $h_{mount} = 2500$ m. Concerning Fig. 4a, we note that curves labelled 0, 0.4, 0.5 and 1 can be compared to the time

<figure><img src="image_3.png" /><figcaption>Fig. 2. Topography and depth (m) of vertical levels. Vertical dashed lines indicate the location of vertical profiles of Fig. 3.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Vertical axes: depth (m). Horizontal axes: $10^7 \times \rho_0^{-1} \nabla_x p'$ (ms$^{-2}$). The terminology is the same as in Table 1: the discrete PGF has been computed with the *Modified-Primitive* scheme and the density has been initialized using the *Volume-Averaged* method (Case 1) or the straightforward method (Case 1B). The discrete PGF has been computed with the *Straightforward-Primitive* scheme and the density has been initialized using the straightforward method (Case 2) or the *Volume-Averaged* method (Case 2B). The discrete PGF has been computed with the *Standard-Jacobian* scheme and the density has been initialized using the straightforward method (Case 3). The location of profiles a–d is indicated by vertical dashed lines in Fig. 2.</figcaption></figure>

<table><caption>**Table 1**
Discrete PGF and density initialization options used in Fig. 3.</caption><thead><tr><th>Cases</th><th>Pressure</th><th>Initialization</th><th>Pressure gradient</th></tr></thead><tbody><tr><td>1</td><td>Rectangular</td><td>Volume-Averaged</td><td>Modified-Primitive (Eq. (13))</td></tr><tr><td>2</td><td>Rectangular</td><td>Straightforward</td><td>Straightforward-Primitive (Eq. (4))</td></tr><tr><td>3</td><td>Trapezoidal</td><td>Straightforward</td><td>Standard-Jacobian</td></tr><tr><td>1B</td><td>Rectangular</td><td>Straightforward</td><td>Modified-Primitive</td></tr><tr><td>2B</td><td>Rectangular</td><td>Volume-Averaged</td><td>Straightforward-Primitive</td></tr></tbody></table>

*Note:* SMW03 showed that the pressure gradients of Cases 2 and 3 become equivalent provided that the rectangular rule used to compute the pressure in Case 2 is replaced by a trapezoidal rule.

history of the *net kinetic energy* computed by SMW03 for the *Blended-Jacobian* PGF (see top left Fig. in SMW03, p35–8).

Fig. 4 finally leads to the following observations:

1. Despite our model is not exactly the same as the one used in SMW03, Fig. 4a leads to similar comments. Indeed, as far as *Density-Jacobian* methods are concerned, the *Weighted-Jacobian* PGF gives the worst results and for a wide range of combinations ($\gamma = 0.5$ being a priori the best option), the *Blended-Jacobian* PGF is more accurate than the *Standard-Jacobian* PGF.

2. The same comments can be done for Fig. 4b–d, except that the best configuration of the *Blended-Jacobian* PGF is not always obtained for $\gamma = 0.5$.

3. The *Modified-Primitive* PGF (13) is better than the *Straightforward-Primitive* PGF (4) and also better than the *Standard-Jacobian* PGF and the *Weighted-Jacobian* PGF. In some cases, the *Modified-Primitive* PGF (13) is even better than the best configuration of the *Blended-Jacobian* PGF (Fig. 4c and d).

4. We note that errors developed by the *Straightforward-Primitive* PGF are not, as Fig. 3 could suggest, stronger than those developed by the *Standard-Jacobian* PGF. Another unexpected result is that errors developed by *Primitive type* PGFs (13) and (4) are not very sensitive to the way density is initialized. We also observe that ERKE of Fig. 4a and b are of the same order, despite *hydrostatic inconsistency*, and thus initial PGF errors, are greater in the case corresponding to Fig. 4b (obtained with 22 vertical levels).

The latter point confirms an idea already underlined by several previous studies, namely that initial PGF errors are not systematically representative of long-term erroneous currents. As a matter of fact, analysing a vorticity balance built from the discrete equations of their model, MOE98 have identified a sigma error of the second kind, more likely responsible for the growth of erroneous circulations. A similar approach is now used to better understand the results shown on Fig. 4 and notably explained the particularly good results provided the *Modified-Primitive* PGF (13).

### 4. Bottom torque consistency

The bottom torque (or topographic torque) concept follows from the curl of the depth integrated momentum equations (Mertz and Wright, 1992). In the ideal case of the seamount test, the latter is expected to be zero. As the sea surface elevation is also supposed to remain unchanged in the ideal case (i.e. $\zeta = 0$), this means that the bottom pressure anomaly, $Pb' = g \int_{-h}^{0} \rho' dz$, should simply verify $J(Pb', h) = 0$. However, this is not likely to happen in our numerical context, because of the trancation errors associated to the discrete vorticity balance. If we consider
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** ERKE ($m^2s^{-2}$) as a function of time (days), for the reference case described in Section 3.1 (a), case with 22 vertical levels (b), case with initial density given by $\rho = -3e^{2/250}$ (c), case with bathymetry given by (15) using $h_{\text{mount}} = 2500$ m (d). Solid lines: PGF has been computed using a *Modified-Primitive* scheme with a *volume averaged* initialization (s1), a *Modified-Primitive* scheme with straightforward initialization (s2), a *Straightforward-Primitive* scheme with *volume averaged* initialization (s3), a *Straightforward-Primitive* scheme with straightforward initialization (s4). Dashed lines: PGF has been computed using the *Blended-Jacobian* method. Labels indicate the value of $\gamma$, i.e. the mixing ratio of the *Standard-Jacobian* and the *Weighted-Jacobian* PGFs. $\gamma = 0$ corresponds to the *Standard-Jacobian* PGF. $\gamma = 1$ corresponds to the *Weighted-Jacobian* PGF.</figcaption></figure>

that the variations of the coriolis parameter, the non-linear and viscosity terms are negligible, the vorticity balance of Mertz and Wright (1992, Eq. (14) p. 303) suggests that the model should

produce an adjustment of the sea surface anomaly in order to restore the equilibrium of the bottom torque, in other words $J(g\zeta, h) = -\rho_0^{-1}J(Pb', h)$.
Let us consider now the depth integrated PGF at the initial state ($\zeta = 0$) and define $\bar{p}'$ as being the depth-averaged pressure anomaly. The vertical integral for each of the last two terms at the RHS of (1) is:

$$ \int_{-h}^{0} \left( \left. \frac{\partial p'}{\partial x} \right|_{s} \right) dz = h \left. \frac{\partial}{\partial x} \int_{-1}^{0} p' ds \right. = h \left. \frac{\partial \bar{p}'}{\partial x} \right. \quad (16) $$

$$ \int_{-h}^{0} \left( g \rho' \left. \frac{\partial z}{\partial x} \right|_{s} \right) dz = \frac{g}{h} \left. \frac{\partial h}{\partial x} \right|_{-h} \int_{-h}^{0} z \rho' dz = \frac{\partial h}{\partial x} \bar{p}' - h \left. \frac{\partial Pb'}{\partial x} \right. \quad (17) $$

where, for convenience, we used the coordinate transformation $s = z/h$ in order to write $dz = h \, ds$ and $\frac{\partial z}{\partial x}|_s = h^{-1} \frac{\partial h}{\partial x} z$. We see that the hydrostatic correction term of the PGF leads to two terms at the RHS of (17). The former, $\frac{\partial h}{\partial x} \bar{p}'$, added to the RHS of (16), gives $\frac{\partial h}{\partial x} \left. \frac{\partial \bar{p}'}{\partial x} \right|_s$. The latter disappears when the curl of the depth integrated PGF is written, in order to simply obtain $J(Pb', h)$. So we clearly see that some numerical consistency between the discrete form of $\left. \frac{\partial p'}{\partial x} \right|_s$ and $g \rho' \left. \frac{\partial z}{\partial x} \right|_s$ is required in order to obtain a realistic discrete bottom torque effect. This property is referred to as the bottom torque consistency. In Appendix we show that the Modified-Primitive PGF is bottom torque consistent, in other words that the curl of the depth integrated PGF (13) has the expected form of a Jacobian of the bottom pressure anomaly and the bathymetry, namely (details of the demonstration are given in Appendix):

$$ \text{curl}_z \left( \int_{-h}^{0} \left( \frac{\nabla p'}{\rho_0} \right) dz \right)_{i-1,2,j-1/2} = \frac{(h_{i-1,j} - h_{i,j-1})(Pb'_{i,j} - Pb'_{i-1,j-1}) - (h_{i,j} - h_{i-1,j-1})(Pb'_{i-1,j} - Pb'_{i,j-1})}{2\rho_0 \Delta x \Delta y} \quad (18) $$

Although there is no particular reason for the discrete Jacobian (18) to be exactly zero in the present case (since $h$ and $Pb'$ are not linear in $x$ and $y$), we can intuitively expect that the curl of the depth integrated PGF will be smaller if the scheme is bottom torque consistent, i.e. if the RHS of (16) and the first term at the RHS of (17) disappear when the discrete vorticity balance is computed. A consequence of this should be a smaller adjustment of the sea surface pressure and thus a smaller long-term surface circulation. As a matter of fact, the bottom torque consistency has already been identified as a critical property in previous studies (Song and Wright, 1998). These authors have proposed a demonstration of the bottom torque consistency for the Standard-Jacobian PGF, but their reasoning, however, required redefining a suitable expression of the discrete bottom pressure (see Eq. (2.37) in Song and Wright, 1998). Thus, it is worthwhile to note that in the case of the Modified-Primitive PGF, our demonstration is based on the genuine discrete pressure (see Appendix). In the following, this will explain why bottom torque computations can be rather different from one scheme to another.

A mean vorticity error, obtained by averaging over the numerical domain the magnitude of the curl of the depth integrated pressure gradient, namely:

$$ \bar{\zeta}_e = \iint \left| \text{curl}_z \left( \int_{-h}^{0} \frac{\nabla p'}{\rho_0} dz \right) \right| dxdy / \iint dxdy \quad (19) $$

is presented on Table 2 for the various PGF. We show in Appendix that the bottom torque consistency of the Modified-Primitive PGF is a consequence of the particular form of its discrete hydrostatic correction term. Table 2 actually shows that the global vorticity effect of the initial PGF is twice as small in the case of the Modified-Primitive PGF than in the case of the Straightforward-Primitive PGF.

Obviously, the conclusions that can be made from Table 2 are partly limited by the fact that our vorticity analysis does not consider the non-linearity of the system. We also note that $\bar{\zeta}_e$ has been computed at the initial time so that Table 2 does not consider the

evolution of the density field either. As a matter of fact, Table 2 is not always in perfect agreement with Fig. 4. If considering the reference case, the lowest $\bar{\zeta}_e$ is obtained with the Blended-Jacobian PGF using $\gamma = 0.4$, and Fig. 4a shows that $\gamma = 0.5$ gives the lowest ERKE. The fact that the Modified-Primitive PGF performs better than the Blended-Jacobian PGF with $\gamma = 0.3$ (Fig. 4a) is also not completely consistent with Table 2.

On the other hand, Table 2 is in good agreement with Fig. 4 on several points:

1. The Primitive type PGF is more sensitive to the discrete form of the hydrostatic term, i.e. (13), (4), than to the way the density field is initialized.

2. The 22 vertical levels case and the reference case lead to similar $\bar{\zeta}_e$.

3. The Straightforward-Primitive and Standard-Jacobian PGFs give similar $\bar{\zeta}_e$.

4. $\bar{\zeta}_e$ is generally much smaller in the case of the Modified-Primitive PGF.

5. Among the considered scheme, the worst $\bar{\zeta}_e$ is obtained with the Weighted-Jacobian PGF.

6. Concerning the Blended-Jacobian PGF, the best combination (i.e. smallest $\bar{\zeta}_e$) is not the same from one topographic situation to another.

7. For $h_{\text{mount}} = 2500m$, the Modified-Primitive PGF gives particularly low $\bar{\zeta}_e$. This can be explained by the fact that both $h$ and $Pb'$ are nearly linear so that Expression (18) nearly vanishes.

Points 2 and 7 lead to the following comment: smoothing bathymetry is generally recommended to reduce hydrostatic inconsistency, but in the case of the Modified-Primitive PGF, this should be mainly regarded as a mean to reduce bottom torque inconsistency.

## 5. Conclusions

Our paper focussed on the currently admitted idea that Density-Jacobian type PGFs perform better than the primitive formulation of the PGF. We first examined the main argument that has led to this unfavourable assessment, namely the fact that the Standard-Jacobian PGF is exact when density is linear in $z$, while the primitive type PGF is not. If this argument is correct as far as the Straightforward-Primitive PGF is concerned, we showed that the discrete formulation of the latter can be improved. We first noticed that previous studies were based on a simplistic concep-

<caption>**Table 2**</caption>

$\Delta x \Delta y \bar{\zeta}_e(m^3 s^{-2})$ at the initial state as a function of the PGF scheme (column 1). (*) indicates that the density field has been initialized with the Volume-Averaged method. Column 2 corresponds to the reference case described in Section 3.1. Column 3: same as column 2 but we used 22 vertical levels. Column 4: same as column 2 but initial density is given by $\rho = -3e^{2/250}$, Column 5: same as column 2 but bathymetry is given by (15) with $h_{\text{mount}} = 2500m$.

<table><tr><td>Modified-Primitive</td><td>0.2667</td><td>0.2676</td><td>0.0583</td><td>0.0020</td></tr><tr><td>Modified-Primitive(*)</td><td>0.2678</td><td>0.2678</td><td>0.0578</td><td>0.0020</td></tr><tr><td>Straightforward-Primitive</td><td>0.4542</td><td>0.4519</td><td>0.2666</td><td>0.0272</td></tr><tr><td>Straightforward-Primitive(*)</td><td>0.4543</td><td>0.4520</td><td>0.2693</td><td>0.0276</td></tr><tr><td>Standard-Jacobian</td><td>0.4509</td><td>0.4511</td><td>0.2643</td><td>0.0270</td></tr><tr><td>Blended γ = 0.1</td><td>0.3421</td><td>0.3427</td><td>0.1609</td><td>0.0151</td></tr><tr><td>Blended γ = 0.2</td><td>0.2371</td><td>0.2373</td><td>0.1279</td><td>0.0053</td></tr><tr><td>Blended γ = 0.3</td><td>0.1704</td><td>0.1692</td><td>0.1567</td><td>0.0093</td></tr><tr><td>Blended γ = 0.4</td><td>0.1573</td><td>0.1546</td><td>0.2191</td><td>0.0204</td></tr><tr><td>Blended γ = 0.5</td><td>0.1903</td><td>0.1865</td><td>0.3022</td><td>0.0322</td></tr><tr><td>Blended γ = 0.6</td><td>0.2504</td><td>0.2464</td><td>0.4056</td><td>0.0441</td></tr><tr><td>Blended γ = 0.7</td><td>0.3354</td><td>0.3312</td><td>0.5091</td><td>0.0559</td></tr><tr><td>Blended γ = 0.8</td><td>0.4290</td><td>0.4251</td><td>0.6132</td><td>0.0677</td></tr><tr><td>Blended γ = 0.9</td><td>0.5333</td><td>0.5297</td><td>0.7229</td><td>0.0796</td></tr><tr><td>Weighted-Jacobian</td><td>0.6405</td><td>0.6366</td><td>0.8326</td><td>0.0914</td></tr></table>
tion (moreover not conservative) of the initialization of the density field. As a matter of fact, a fundamental concept in numerical modelling is that a discrete variable does not give a punctual but a mean representation of the "true" field around each grid point. Taking this basic principle into consideration ensures the consistency of the density field with the rectangular scheme used to compute the hydrostatic pressure and on the other hand leads to reconsider the discrete form of the hydrostatic correction term $g\rho' \frac{\partial z_i}{\partial x}$. The so-called *Modified-Primitive* scheme reduces the initial PGF errors.

Our paper has also examined the impact of the PGF discretization on the vorticity balance and showed that the *Modified-Primitive* PGF is consistent with the vorticity balance built from the discrete transport equations and leads to smaller vorticity errors. These errors (somehow recalling the sigma errors of the second kind reported by MOE98) are well correlated with the long-term behaviour of the different tested PGFs, and thus appear to be more relevant than the initial PGF errors. We will principally retain from these tests that long-term errors are smaller with the *Modified-Primitive* PGF than with the *Standard-Jacobian* or *Weighted-Jacobian* PGFs. In some cases, the hybrid PGF built from the combination of the *Standard-Jacobian* and *Weighted-Jacobian* PGFs, gives the best results, but the fact that the optimal value of $\gamma$ (the ratio of the two schemes), can vary from one case to the other, may lead to reconsider the relevance of the so-called *Blended-Jacobian* PGF.

However, these conclusions should be tempered by the fact that the seamount test hardly exemplifies all possible ocean situations, notably concerning the background stratification. For instance, the *Standard-Jacobian* PGF should be a better option than the *Modified-Primitive* PGF when the density is linear in $z$ and horizontally homogeneous. Some authors (Haney, 1991) have suggested to remove the purely vertical stratification before computing the PGF. This would, however, be meaningless in the case of the seamount test since all PGF schemes would exactly vanish. This is, however, an important issue for more realistic cases, since the remaining stratification (once the purely vertical stratification is removed from the density field) may significantly differ from the exponential profile used in the seamount test.

Higher order PGFs are beyond the scope of this study. Readers can refer to Chu and Fan (1997) who showed that 4th and 6th order schemes can easily be derived for Primitive type PGFs. On the other hand, the possibility to ensure the consistency of the PGF with some other important physical concepts, makes low-order schemes also attractive and, to some extent, balances their a priori lower accuracy. Apart from the aforementioned mass conservation and *bottom torque* consistencies, the energetic consistency with the density advection scheme is easily achieved with a low-order PGF. If, as noticed by SMW03, the interest of such a property is possibly doubtful in the case of the seamount test and its simplistic ideal solution (a motionless ocean), the conservation of energy is a fundamental issue of processes like gravity waves. Along these lines, our study proposed a suitable discrete form of the density advection consistent with the *Modified-Primitive* PGF.

## Acknowledgements

This study was funded by the European MFSTEP Project (EU Contract EVK3-CT-2002-00075), the European INSEA Project (SST4-CT-2005-012336), the LEFE-IDAO-2006 and Mercator OPA-Symphonique French projects. The authors thank Cyril Nguyen and the Laboratoire d'Aérologie computer team, Serge Prieur, Laurent Cabanas, Jérémy Leclercq, Didier Gazen and Juan Escobar for their support.

# Appendix A. The initial ($\zeta = 0$) bottom pressure torque

Considering (5) and (6), we have:

$$z_{i,k} = -h_i + \sum_{q=1,k-1} \Delta z_{i,q} + 0.5 \Delta z_{i,k} \quad (A1)$$

where $h$ is the bathymetry.

We now consider the numerical vertical integral $\sum_{k=1,k_{\max}} p'_{i,k} \Delta z_{i',k}$, where the index $i'$ can be arbitrarily fixed to $i' = i - 1$ or $i' = i + 1$. Using (7) and (8), we have:

$$\sum_{k=1,k_{\max}} p'_{i,k} \Delta z_{i',k} = \sum_{k=1,k_{\max}} \left[ \sum_{q=k+1,k_{\max}} \rho'_{i,q} \Delta z_{i,q} + 0.5 \rho'_{i,k} \Delta z_{i,k} \right] \Delta z_{i',k} \quad (A2)$$

The summation at the RHS of (A2) can be alternatively written:

$$\sum_{k=1,k_{\max}} p'_{i,k} \Delta z_{i',k} = \sum_{k=1,k_{\max}} \left[ g \rho'_{i,k} \Delta z_{i,k} \left( \sum_{q=1,k-1} \Delta z_{i',q} + 0.5 \Delta z_{i',k} \right) \right] \quad (A3)$$

Using (A1) and (A3) is equivalent to:

$$\sum_{k=1,k_{\max}} p'_{i,k} \Delta z_{i',k} = \sum_{k=1,k_{\max}} \left[ g \rho'_{i,k} \Delta z_{i,k} (h_{i'} + z_{i',k}) \right] \quad (A4)$$

and finally:

$$\sum_{k=1,k_{\max}} p'_{i,k} \Delta z_{i',k} = h_{i'} P b'_{i} + \sum_{k=1,k_{\max}} g z_{i',k} \rho'_{i,k} \Delta z_{i,k} \quad (A5)$$

where $Pb'_{i} = g \sum_{k=1,k_{\max}} \rho'_{i,k} \Delta z_{i,k}$ is the discrete bottom pressure anomaly.

The curl of the depth integrated pressure anomaly gradient is of the form:

$$(PX_{i-1/2,j-1} - PX_{i-1/2,j}) / \Delta y + (PY_{i-1/2} - PY_{i-1/2,j-1}) / \Delta x \quad (A6)$$

where $PX$ and $PY$ are the $O_x$ and $O_y$ components of the depth integrated PGF (the location of the related variables is given in Fig. 5). Using the *Modified-Primitive* PGF scheme (13), we have:

$$PX_{i-1/2,j} = \frac{-1}{\Delta x \rho_0} \sum_{k=1,k_{\max}} \left[ \left( p'_{i,j,k} - p'_{i-1,j,k} + g(z_{i,j,k} - z_{i-1,j,k}) \right) \times \frac{\Delta z_{i,j,k} \rho_{i,j,k} + \Delta z_{i-1,j,k} \rho_{i-1,j,k}}{\Delta z_{i,j,k} + \Delta z_{i-1,j,k}} \right] \frac{\Delta z_{i,j,k} + \Delta z_{i-1,j,k}}{2} \quad (A7)$$

which can be also written:

$$PX_{i-1/2,j} = \frac{-1}{2 \Delta x \rho_0} \sum_{k=1,k_{\max}} \left[ \begin{aligned} & p'_{i,j,k} \Delta z_{i-1,j,k} - p'_{i-1,j,k} \Delta z_{i,j,k} + g z_{i,j,k} \Delta z_{i-1,j,k} \rho_{i-1,j,k} \\ & + p'_{i,j,k} \Delta z_{i-1,j,k} - p'_{i-1,j,k} \Delta z_{i,j,k} + g z_{i,j,k} \Delta z_{i-1,j,k} \rho_{i-1,j,k} - g z_{i-1,j,k} \Delta z_{i,j,k} \rho_{i,j,k} \end{aligned} \right] \quad (A8)$$

Now, using (A5) with $i' = i - 1$ and $i' = i + 1$, the last four terms at the RHS of (A8) are rewritten:

$$\sum_{k=1,k_{\max}} \left[ \begin{aligned} & p'_{i,j,k} \Delta z_{i-1,j,k} - p'_{i-1,j,k} \Delta z_{i,j,k} + g z_{i,j,k} \Delta z_{i-1,j,k} \rho_{i-1,j,k} \\ & - g z_{i-1,j,k} \Delta z_{i,j,k} \rho_{i,j,k} \end{aligned} \right] = -h_{i,j} P b'_{i-1,j} + h_{i-1,j} P b'_{i,j} \quad (A9)$$

Using (A9) and (A8) is rewritten:

$$PX_{i-1/2,j} = \frac{-1}{2 \Delta x \rho_0} \sum_{k=1,k_{\max}} \left[ \begin{aligned} & p'_{i,j,k} \Delta z_{i,j,k} - p'_{i-1,j,k} \Delta z_{i-1,j,k} + g z_{i,j,k} \Delta z_{i,j,k} \rho_{i,j,k} \\ & - g z_{i-1,j,k} \Delta z_{i-1,j,k} \rho_{i-1,j,k} \end{aligned} \right] - \frac{-h_{i,j} P b'_{i-1,j} + h_{i-1,j} P b'_{i,j}}{2 \Delta x \rho_0} \quad (A10)$$

Similarly, $PY_{i-1/2,j}$ can be written:
<figure><img src="image_6.png" /><figcaption>Fig. 5. Location of variables on the horizontal C-grid. Dotted lines indicates the main axes of the discrete Jacobian J(Pb', h).</figcaption></figure>

$$ PY_{i,j-1/2} = \frac{-1}{2\Delta y \rho_0} \sum_{k=1, k_{\max}} [p'_{i,j,k} \Delta z_{i,j,k} - p'_{i,j-1,k} \Delta z_{i,j-1,k} \\ + g z_{i,j,k} \Delta z_{i,j,k} \rho_{i,j,k} - g z_{i,j-1,k} \Delta z_{i,j-1,k} \rho_{i,j-1,k}] \\ - \frac{-h_{i,j} Pb'_{i,j-1} + h_{i,j-1} Pb'_{i,j}}{2 \Delta y \rho_0} \quad (A11) $$

Substituting A10 and A11 into (A6) we finally obtain:

$$ \frac{(PX_{i-1/2,j-1} - PX_{i-1/2,j})}{\Delta y} + \frac{(PY_{i-1/2} - PY_{i-1/2})}{\Delta x} = \frac{(h_{i-1,j} - h_{i,j-1})(Pb'_{i,j} - Pb'_{i-1,j-1}) - (h_{i,j} - h_{i-1,j-1})(Pb'_{i-1,j} - Pb'_{i,j-1})}{2\rho_0 \Delta x \Delta y} \quad (A12) $$

which can be viewed as the numerical counterpart of the bottom torque term, $\rho_0^{-1}J(Pb', h)$, appearing in Mertz and Wright (1992, Eq (14), p303). Let us mainly retain here that:

1. The Modified-Primitive scheme (13) is bottom torque consistent since the discrete curl of the depth integrated PGF has, as expected, the form of a Jacobian of the bottom pressure and the topography.

2. The terms of the form g z Δz ρ, appearing in (A7)–(A11), are due to the density weighted-average used in the hydrostatic correction term of (13).

3. The bottom torque consistency relies on these very terms, since the cancellation of the terms appearing on the first line of (A10) and (A11), (when the curl of the depth integrated pressure is computed), is a consequence of the verification of (A5).

References

Auclair, F., Marsaleix, P., Estournel, C., 2000. Sigma coordinate pressure gradient errors: evaluation and reduction by an inverse method. Journal of Atmospheric and Oceanic Technology 17, 1348–1367.

Beckmann, A., Haidvogel, D.B., 1993. Numerical simulation of flow around a tall isolated seamount. Part I: problem formulation and model accuracy. Journal of Physical Oceanography 23, 1736–1753.

Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal circulation model. In: Heaps, N. (Ed.), Three-Dimensional Coastal Ocean Models. Coastal Estuarine Science, vol. 4. American Geophysical Union, pp. 1–16.

Chu, P.C., Fan, C., 1997. Sixth-order difference scheme for sigma coordinate ocean models. Journal of Physical Oceanography 27, 2064–2071.

Ciappa, A.C., 2008. A method for reducing pressure gradient errors improving the sigma coordinate stretching function: an idealized flow patterned after Libyan near-shore region with the POM. Ocean Modelling 23, 59–72.

Ferziger, J.H., Peric, M., 2002. Computational Methods for Fluid Dynamics, third ed. Springer. p. 423.

Gaspar, P., Gregoris, Y., Lefevre, J.M., 1990. A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: tests at station Papa and long-term upper ocean study site. Journal of Geophysical Research 95, 16,179–16,193.

Haney, R.L., 1991. On the pressure gradient force over steep topography in sigma coordinate ocean models. Journal of Physical Oceanography 21, 610–619.

James, I.D., 1996. Advection schemes for shelf sea models. Journal of Marine Systems 8, 237–254.

Marsaleix, P., Auclair, F., Floor, J.W., Herrmann, M.J., Estournel, C., Pairaud, I., Ulises, C., 2008. Energy conservation issues in sigma-coordinate free-surface ocean models. Ocean Modelling 20, 61–89.

Mellor, G.L., Ezer, T., Oey, L.-Y., 1994. The pressure gradient conundrum of sigma coordinate ocean models. Journal of atmospheric and oceanic technology 11, 1126–1134.

Mellor, G.L., Oey, L.-Y., Ezer, T., 1998. Sigma coordinate pressure gradient errors and the seamount problem. Journal of atmospheric and oceanic technology 15, 1122–1131.

Mertz, G., Wright, D.G., 1992. Interpretations of the Jebar term. Journal of Physical Oceanography 22, 301–305.

Shchepetkin, A.F., McWilliams, J.C., 2003. A method for computing horizontal pressure-gradient force in an oceanic model with non-aligned vertical coordinate. Journal of Geophysical Research 108 (3090), 35.1–35.34.

Song, Y.T., 1998. A general pressure gradient formulation for ocean models. Part I: scheme design and diagnostic analysis. Monthly Weather Review 126, 3213–3230.

Song, Y.T., Wright, D.G., 1998. A general pressure gradient formulation for ocean models. Part II: energy, momentum, and bottom torque consistency. Monthly Weather Review 126, 3231–3247.