# An accurate implementation of the compressibility terms in the equation of state in a low order pressure gradient scheme for sigma coordinate ocean models

Patrick Marsaleix, Francis Auclair, Claude Estournel, Cyril Nguyen, Caroline Ulses

► To cite this version:

Patrick Marsaleix, Francis Auclair, Claude Estournel, Cyril Nguyen, Caroline Ulses. An accurate implementation of the compressibility terms in the equation of state in a low order pressure gradient scheme for sigma coordinate ocean models. Ocean Modelling, 2011, 40 (1), pp.1-13. 10.1016/j.ocemod.2011.07.004. hal-02110128

HAL Id: hal-02110128

https://hal.science/hal-02110128

Submitted on 29 Sep 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 Attribution 4.0 International License
# An accurate implementation of the compressibility terms in the equation of state in a low order pressure gradient scheme for sigma coordinate ocean models

Patrick Marsaleix<sup>a,*</sup>, Francis Auclair<sup>c</sup>, Claude Estournel<sup>a</sup>, Cyril Nguyen<sup>b</sup>, Caroline Ulses<sup>c</sup>

<sup>a</sup>Laboratoire d'Aérologie, CNRS, 14 Avenue Edouard Belin, F-31400 Toulouse, France

<sup>b</sup>Observatoire Midi-Pyrénées, Université de Toulouse, 14 Avenue Edouard Belin, F-31400 Toulouse, France

<sup>c</sup>Laboratoire d'Aérologie, Université de Toulouse, 14 Avenue Edouard Belin, F-31400 Toulouse, France

In a previous study, the authors studied a low order pressure gradient force (PGF) scheme referred to as the *Primitive-Modified* scheme, that appears to be equivalent to the *Pressure-Jacobian* PGF (Lin, 1997). The scheme was successfully tested on the seamount experiment using a simplified equation of state (EOS). Yet, a complete equation of state, including compressibility effect, can raise a serious problem of accuracy. A new implementation is thus proposed in the present paper. The scheme is rewritten using a (numerically equivalent) geopotential formulation. The PGF truncation errors are removed by computing the EOS compressibility terms with potential temperature and salinity interpolated on a suitable geopotential level. The so-called *Equivalent Geopotential Formulation* (EGF) method is compared to the Finite-Volume approach proposed by Adcroft et al. (2008).

## 1. Introduction

Marsaleix et al. (2009) (hereafter M09) present a low order pressure gradient force (hereafter PGF) scheme referred to as the "*Primitive Modified*" (hereafter PM) scheme. We will see in the following that this scheme is equivalent to the *Pressure-Jacobian* PGF (Lin, 1997). At first sight the order of the PM scheme is too low to be really appropriate for ocean modelling. As a matter of fact, it does not pass the academic test of an "*ocean at rest*" unless the density field is constant, when on the other hand the wide spread *Standard-Jacobian* PGF exactly cancels for a density profile linear in *z* (Shchepetkin and McWilliams, 2003, hereafter SM03). However the authors more specifically retained two arguments in favour of the PM scheme.

First of all the discrete form used for the hydrostatic pressure is consistent with the conservation property of tracer advection and diffusion schemes, generally adopted by up-to-date ocean models. According to M09, the tracer conservation property naturally leads to consider that discrete values of potential temperature and salinity (hereafter $\theta$ and $S$) represent averaged values of "true" ($\theta$, $S$) fields over the volume of cell boxes. Assuming this formulation of the discrete ($\theta$, $S$) fields and using a simplified equation of state (neglecting non-linearity and compressibility), M09 show that the discrete pressure at the upper and lower facets of the C-grid cell

boxes is exact if the hydrostatic integral of the sea water density is computed with the “rectangular” method. Along these lines, M09 show that a substantial part of the initial PGF errors found in sigma coordinate models is cancelled if the initial ($\theta$, $S$) discrete fields are obtained by averaging the a priori true fields over the cell boxes volumes.

The second argument is the particular form of the hydrostatic correction term used to balance the pressure variation induced by the slope of grid levels. This term notably appears in sigma coordinate models or, more generally, when the grid levels are not strictly horizontal. In the following, the subscript $s$ will refer to a differentiation along a constant level in the topography-following coordinate system and the subscript $z$ to a differentiation along a constant geopotential level. Depending on models, a possible continuous formulation for the slope correction term is $g\rho'\frac{\partial z}{\partial x}|_{s}$. As far as the C-grid is concerned, the density anomaly $\rho'$ is approximated by the average of the discrete density values bounding the momentum equation grid node. Avoiding the straightforward centred average of Janjic (1977), the specificity of the PM scheme is to use a weighted average depending on the vertical size of cell boxes. M09 show that doing so enforces the accuracy of the bottom pressure torque. In case of strong interaction between the current and the bathymetry, the bottom pressure torque is actually an important mechanism of the depth-averaged flow equilibrium (Mertz and Wright, 1992). The coastal zone and the well-known seamount experiment are obviously particularly concerned by this issue. In the case of the seamount experiment, revisited by numerical

<footnote>* Corresponding author. Tel.: +33 5 61 33 27 63; fax: +33 5 61 33 27 90. E-mail address: marp@aeroobs-mip.fr (P. Marsaleix).</footnote>
studies dedicated to the PGF accuracy, the errors of the second kind
are rather expectedly responsible for the long term growth of erro-
neous currents.

In their previous study, the authors used a simplified pressure-
independent equation of state. Yet, the pressure dependent terms
of the equation of state introduce an important and non-linear
dependency of the sea water density to the depth. Thus, it is a po-
tential source of PGF errors in terrain following coordinate models.
In order to solve this problem, modellers have proposed different
methods. A classical technique consists in subtracting a back-
ground vertical profile of density (Mellor et al., 1994) from the
3D density field before computing the PGF. Another one, suggested
by SM03, consists in removing from the equation of state the terms
that depend on z only. Another option is to maintain some consis-
tency between the accuracy of the PGF and the non-linearity of the
stratification by using a high order PGF scheme. Recently, Adcroft
et al. (2008) have proposed a Finite Volume approach. In this meth-
od, the PGF involves the computation of the integral of the pressure
over the edges of the momentum cell boxes. Provided that a trac-
table equation of state has been chosen, this integral can be calcu-
lated analytically. Until now, we may consider that Adcroft et al.
(2008) proposed the most sophisticated, possibly accurate, method
to deal with the compressibility of sea water.

Along these lines, the case of the *Pressure-Jacobian* PGF is partic-
ularly relevant since (as shown in the following) its accuracy ap-
pears to be quite sensitive to the pressure-dependent terms of
the equation of state. On the other hand, the objective of the pres-
ent study is not to focus on a particular PGF scheme but to propose
an original and simple method in order to reduce the PGF errors re-
lated to the compressibility of sea water. We will see that this
method is largely independent of the chosen PGF scheme. If the

Pressure-Jacobian PGF is naturally considered, since it perfectly evidences the benefit of an accurate treatment of the compressibility of sea water, the present study will also consider other well-known PGF schemes.

Our method, which will be detailed in the following, can be summarized as follows (see also Fig. 1). A first step consists in rewriting the PGF using a (numerically equivalent) geopotential formulation inspired from SM03. A geopotential level of reference is then identified, on which the (θ, S) dependent parameters of the sea water density problem are interpolated. The equation of state is then calculated in this geopotential frame. The PGF is finally computed.

As our method also intends to be EOS independent, this study considers several equations of state, including the one recom-mended by the most recent TEOS10 report (IOC, SCOR and IAPSO, 2010). Comparisons with the modern approach of Adcroft et al. (2008) will be also presented.

The paper is organized as follows: we first present the several
PGF schemes that will be used to test our method. For each of
them, we suggest an alternative geopotential formulation, numer-
ically equivalent to the usual discretization (as far as the potential
density is concerned), but offering a straightforward opportunity to
solve the problem of the PGF accuracy related to the compressibil-
ity terms of the equation of state. The latter is exposed through the
examination of the EOS proposed by Mellor (1991). Practical solu-
tions to correctly implement the compressibility effect are then
presented in the case of several EOS (Mellor, 1991; Wright, 1997;
McDougall et al., 2003). We compare their respective performance
in the case of the seamount experiment. We also compare the so-
called *Equivalent Geopotential Formulation* (EGF) method to the Fi-
nite Volume approach of Adcroft et al. (2008).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Three implementations of the Equivalent Geopotential Formulation approach for the Density-Jacobian. The density anomaly is interpolated on the geopotential level of reference as in SM03 (a). First, the temperature and the salinity are interpolated on the geopotential level of reference and then the Mellor (1991) EOS is computed (b). First, the coefficients of the Wright (1997) or the McDougall et al. (2003) EOS are computed on the tracer grid nodes and then interpolated on the geopotential level of reference, where the density anomaly is finally computed (c).</figcaption></figure>
## 2. Presentation of different PGF schemes

We recall that the PGF is not, as such, the centre of interest of this study. We consequently do not focus on a particular PGF scheme. Several schemes are in fact considered with the aim of evidencing the benefit of a new method for the treatment of the depth-dependent terms of the equation of state. We chose to deal with low order PGF schemes, easily implementable in an ocean model. Other schemes could have been chosen without changing the principle of our method.

### 2.1. Janjic (1997) scheme and its modified M09 version

The PGF scheme of models using a topography-following coordinate system can be formulated as:

$$-\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} \Big|_s \quad (1)$$

where $p$ is the pressure, $\zeta$ the free surface elevation, $\rho_0$ a constant density of reference, $\rho' = \rho - \rho_0$ with $\rho$ the density of sea water and $p'$ a hydrostatic pressure anomaly related to the density anomaly $\rho'$. The Janjic scheme, given by (2), and the M09 scheme, given by (3), are rather straightforward discrete forms of the last two terms at the right hand side of (1), namely:

$$\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{\rho'_{i,k} + \rho'_{i-1,k}}{2} \frac{z_{i,k} - z_{i-1,k}}{\Delta x} \quad (2)$$

$$\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}}{\Delta z_{i,k} + \Delta z_{i-1,k}} \frac{z_{i,k} - z_{i-1,k}}{\Delta x} \quad (3)$$

where $(p'_{i,k}, \rho'_{i,k}, z_{i,k})$ are the discrete values of, respectively, pressure anomaly, density anomaly (i.e. $\rho'_{i,k} = \rho_{i,k} - \rho_0$) and depth, in the middle of the grid box corresponding to horizontal index $i$ and vertical index $k$. Vertical resolution and middle grid box depth are deduced from the depth of the upper and lower facets, i.e. $z_{i,k} = (z_{i,k+1/2} + z_{i,k-1/2})/2$ and $\Delta z_{i,k} = z_{i,k+1/2} - z_{i,k-1/2}$ (see also Fig. 2). Similarly, the middle grid box pressure is given by:

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

The pressure at grid box facets is given by the hydrostatic integral of the sea water density anomaly using the rectangular method:

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

where $k_{\max}$ corresponds to the vertical index of the surface level. In comparison to its primitive form, proposed by Janjic (1977), where the sea water density at the momentum equation grid node is approximated by the half-half average $(\rho'_{i,k} + \rho'_{i-1,k})/2$, the modified M09 scheme (3), is based on a weighted average depending on the vertical size of respective grid boxes, i.e. $(\Delta z_{i,k} \rho'_{i,k} + \Delta z_{i-1,k} \rho'_{i-1,k})/ (\Delta z_{i,k} + \Delta z_{i-1,k})$. As shown by M09, this leads to a more accurate representation of the bottom torque and thus limits the growth of long term errors in the case of the seamount experiment.

### 2.2. Pressure-Jacobian PGF

Considering (4) and noting that (5) leads to $p'_{i,k-1/2} - p'_{i,k+1/2} = g \rho'_{i,k} \Delta z_{i,k}$ it is possible to rewrite (3) in the following manner:

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

<figure><img src="image_3.png" /><figcaption>Fig. 2. Location of variables of the Equivalent Geopotential Formulation of the Pressure-Jacobian on the vertical C-grid.</figcaption></figure>
This expression is characterized by a symmetry that is also found in the Pressure-Jacobian PGF of Lin (1997) (hereafter L97). Despite its context restricted to the atmospheric case, L97 points up some arguments corroborating the M09 approach. L97 notably defends the idea that a physically motivated approach, even if it is based on a simple finite-volume method, can be competitive compared to some analytic, yet sophisticated, methods. L97 also insists on the consistency that should be maintained with the discrete density. Practically, and as in M09, the latter should represent a mean density between the vertices. As a consequence, the tracer equations should preferably use a finite-volume advection diffusion scheme predicting the volume mean variables. Since (3) and (6) are numerically equivalent, both will be referred to as Pressure-Jacobian in the following.

## 2.3. Standard Density-Jacobian PGF

The Standard Density-Jacobian scheme is based on the following continuous expression for the PGF:

$$-\frac{1}{\rho_0} \left. \frac{\partial p}{\partial x} \right|_z = -\frac{1}{\rho_0} \left( \rho_s g \frac{\partial \zeta}{\partial x} - g \int_z^\zeta J(\rho, z) ds \right) \quad (7)$$

where $\rho_s$ is the density at the surface and where

$$J(\rho, z) = \left. \frac{\partial \rho}{\partial x} \right|_s \frac{\partial z}{\partial s} - \left. \frac{\partial \rho}{\partial s} \frac{\partial z}{\partial x} \right|_s \quad (8)$$

A usual discretization of $J(\rho, z)$ is given by Song (1998), namely:

$$-\Delta x \Delta s J(\rho, z) = \frac{1}{4} \left\{ (\rho'_{i,k+1} + \rho'_{i,k} - \rho'_{i-1,k+1} - \rho'_{i-1,k}) \times (Z_{i,k+1} - Z_{i,k} + Z_{i-1,k+1} - Z_{i-1,k}) - (Z_{i,k+1} + Z_{i,k} - Z_{i-1,k+1} - Z_{i-1,k}) \times (\rho'_{i,k+1} - \rho'_{i,k} + \rho'_{i-1,k+1} - \rho'_{i-1,k}) \right\} \quad (9)$$

where the position of the variables on the C grid is shown in Fig. 1.

## 2.4. Equivalent geopotential formulations

### 2.4.1. The Pressure-Jacobian case

Analogously to SM03, who proposed a reformulation of the Standard Density-Jacobian scheme in a geopotential frame (see Eq. (2.6) in SM03), an alternative (but numerically equivalent) discretization for the right hand side of (6) is obtained by differencing the hydrostatic pressure along a geopotential level. Practically, the new formulation of the Pressure-Jacobian PGF at the $(i - 1/2, k)$ momentum grid node is given by:

$$\frac{p_{i,k}^{*R} - p_{i-1,k}^{*L}}{\Delta x} \quad (10)$$

where, as shown in Fig. 2, $(p_{i,k}^{*R}, p_{i-1,k}^{*L})$ are two pressures computed at a same depth $z_{i-1/2,k}^*$ equal to:

$$z_{i-1/2,k}^* = \frac{Z_{i,k+1/2} Z_{i-1,k+1/2} - Z_{i,k-1/2} Z_{i-1,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}} \quad (11)$$

The height $z_{i-1/2,k}^*$ actually corresponds to the height of intersection of diagonals of the trapezoidal element bounded by vertices $z_{i-1,k-1/2}$, $Z_{i,k-1/2}$, $Z_{i,k+1/2}$, and $Z_{i-1,k+1/2}$ (see Fig. 2). The hydrostatic pressures at depth $z_{i-1/2,k}^*$ are given by:

$$\begin{aligned} p_{i-1,k}^{*L} &= g \sum_{k'=k,k_{\max}-1} \rho_{i-1,k'+1/2}^{*L} \Delta z_{i-1/2,k'+1/2}^{*} + p_{i-1,k_{\max}}^{*L} \\ p_{i,k}^{*R} &= g \sum_{k'=k,k_{\max}-1} \rho_{i,k'+1/2}^{*R} \Delta z_{i-1/2,k+1/2}^{*} + p_{i,k_{\max}}^{*R} \end{aligned} \quad (12)$$

with

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

The sea water density is a weighted average of the surrounding discrete values, namely:

$$\begin{aligned} \rho_{i,k+1/2}^{*R} &= \alpha_{i,k+1/2}^{*R} \rho'_{i,k+1} + (1 - \alpha_{i,k+1/2}^{*R}) \rho'_{i,k} \\ \rho_{i-1,k+1/2}^{*L} &= \alpha_{i-1,k+1/2}^{*L} \rho'_{i-1,k+1} + (1 - \alpha_{i-1,k+1/2}^{*L}) \rho'_{i-1,k} \end{aligned} \quad (14)$$

where the weights of the combination are given by:

$$\begin{aligned} \alpha_{i,k+1/2}^{*R} &= \frac{(z_{i-1/2,k+1}^* - Z_{i,k+1/2}^*)}{(z_{i-1/2,k+1}^* - Z_{i-1/2,k}^*)} \\ \alpha_{i-1,k+1/2}^{*L} &= \frac{(z_{i-1/2,k+1}^* - Z_{i-1,k+1/2}^*)}{(z_{i-1/2,k+1}^* - Z_{i-1,k}^*)} \end{aligned} \quad (15)$$

We note that $0 < \alpha_{i-1,k+1/2}^{*L} < 1$ and $0 < \alpha_{i,k+1/2}^{*R} < 1$ provided that $z_{i-1/2,k}^* < Z_{i-1,k+1/2}^* < z_{i-1/2,k+1}^*$ and $z_{i-1/2,k}^* < Z_{i,k+1/2}^* < z_{i-1/2,k+1}^*$ respectively, the interpolation of the density turning into extrapolation otherwise. The first level under the sea surface, $k = k_{\max}$, is a particular case where we have $p_{i,k_{\max}}^{*R} = g \rho_{i,k_{\max}} (\zeta_i - z_{k_{\max}}^*)$ and $p_{i-1,k_{\max}}^{*L} = g \rho_{i-1,k_{\max}} (\zeta_i - z_{k_{\max}}^*)$.

As shown in Fig. 2, each grid node $(i, k)$ leads to the computation of two pressures, $p_{i,k}^{*R}$ and $p_{i,k}^{*L}$, respectively involved in the computation of the pressure gradient at grid nodes $(i - 1/2, k)$ and $(i + 1/2, k)$. We note that $p_{i,k}^{*R}$ and $p_{i,k}^{*L}$ have no particular reason to be equal.

### 2.4.2. The Janjic (1977) case

The equivalent geopotential formulation of (2) is very similar to that of the Pressure-Jacobian scheme except for the depth of the geopotential level $z^*$ (given by (11) in the previous case) which is now given by:

$$z_{i-1/2,k}^* = \frac{Z_{i-1,k} + Z_{i,k}}{2} \quad (16)$$

that is the average of the level depth corresponding to the discrete density values $\rho_{i-1,k}$ and $\rho_{i,k}$. The other relations given for the Pressure-Jacobian case remain the same.

### 2.4.3. The Standard Density Jacobian case

An Equivalent Geopotential Formulation of the Standard Density-Jacobian scheme is given by SM03. We briefly recall it (details and useful comments can be found in SM03) but using the notations specific to the present study. The discrete form (9) is numerically equivalent to

$$-\Delta x \Delta s J(\rho, z) = -A(\rho^{*R} - \rho^{*L}) \quad (17)$$

where

$$A = \Delta x \frac{Z_{i,k+1} - Z_{i,k} + Z_{i-1,k+1} - Z_{i-1,k}}{2} \quad (18)$$

The discrete along-geopotential differencing involves $\rho^{*L}$ and $\rho^{*R}$ obtained from a linear vertical interpolation (or extrapolation) of the surrounding discrete densities to a common level $z^*$, namely

$$\rho_{i-1,k+1/2}^{*L} = \frac{\rho'_{i-1,k}(Z_{i-1,k+1} - Z^*) + \rho'_{i-1,k+1}(Z^* - Z_{i-1,k})}{Z_{i-1,k+1} - Z_{i-1,k}} \quad (19)$$

$$\rho_{i,k+1/2}^{*R} = \frac{\rho'_{i,k}(Z_{i,k+1} - Z^*) + \rho'_{i,k+1}(Z^* - Z_{i,k})}{Z_{i,k+1} - Z_{i,k}} \quad (20)$$

The common level $z^*$ is given by:

$$z_{i-1/2,k+1/2}^* = \frac{Z_{i,k+1} Z_{i-1,k+1} - Z_{i,k} Z_{i-1,k}}{Z_{i,k+1} - Z_{i,k} + Z_{i-1,k+1} - Z_{i-1,k}} \quad (21)$$

As noticed by SM03, the latter corresponds to the level of intersection of diagonals of the trapezoidal element depicted by Fig. 1a. We note that this expression is analogous to that of the Pressure-Jacobian given by (11), the depth of the density grid nodes, e.g. $(Z_{i,k+1}, Z_{i-1,k+1}, Z_{i,k}, Z_{i-1,k})$, now replacing that of the cellphone facets, e.g. $(Z_{i,k+1/2}, Z_{i-1,k+1/2}, Z_{i,k-1/2}, Z_{i-1,k-1/2})$.
### 3. Accuracy issues related to the type of equation of state

M09 used a equation of state ignoring pressure terms. This approximation may be defendable in the case of a coastal ocean model because the depth of continental shelves is usually small enough to neglect the compressibility effect on sea water density. Over a large domain, including an abyssal area, a complete equation of state can be essential.

As far as the Pressure-Jacobian scheme is concerned, the pressure computed at the upper and lower facets of cell boxes is exact, what ever the complexity of the density profile, because the rectangular method used to compute the hydrostatic integral (5) is consistent with the discrete density field provided that the latter is considered "true on average over the cell boxes volumes" (M09). However the pressure at the middle of grid boxes, given by (4), a half-half average of discrete pressure at facets, is not exact unless the vertical pressure profile is linear (i.e. if the density is constant). Yet, in the case of a complete equation of state, a strongly non-linear pressure profile is normally expected, due to the pressure effect on the density profile, causing a major source of PGF errors when the Pressure-Jacobian scheme is used in a straightforward way.

On the other hand, the density $ρ(θ, S, p)$ can be expressed as the sum of $ρ^0 = ρ(θ, S, p_r)$, a potential density computed with the equation of state at a pressure of reference $p_r$ (for convenience we arbitrarily chose $p_r = 0$), and $δρ^C(θ, S, p)$ a density anomaly related to the compression terms of the equation of state, namely:

$$ρ(θ, S, p) = ρ^0(θ, S) + δρ^C(θ, S, p) \quad (22)$$

The equation of state of Mellor (1991) was, to our knowledge, the first equation relying on the potential temperature (rather than the temperature) and thus particularly well suited to ocean models. The contribution of the compression terms to the Mellor (1991) EOS is given by:

$$δρ^C = 10^4 \frac{p}{c^2} (1 - 0.2 \frac{p}{c^2}) \\ c = 1449.2 + 1.34(S - 35) + 4.55θ - 0.045θ^2 + 0.00821p + 15 \times 10^{-9}p^2 \quad (23)$$

where $p$ is the pressure in decibars. Practically, the latter is not precisely known when the equation of state is computed by the ocean model, since the pressure itself requires the knowledge of the density. Following several other studies (SM03), $p$ is approximated to $-z$ (i.e. assuming $p \approx -10^{-4}gρ_0z$ with $10^{-4}gρ_0 \approx 1$). Besides we observed that using $p$ rather than $-z$ in the EOS does not produce noticeable changes in the results of the realistic experiments presented in Section 6.2. Thus a main characteristic of $δρ^C$ is its strong dependency on depth. According to (23) and using constant temperature and salinity, the variation of $δρ^C$ between the sea surface and the bottom of the abyssal plane ($z \approx -5000$ m) is of the order of $20 \text{ kg m}^{-3}$. Considering the sea water properties presented by Gill (1982, Fig. 3.2), this is sensibly stronger than the variations of the potential density, which tends to become homogeneous at abyssal depth. Conversely, variations of $δρ^C$ at a constant depth are often weaker than those of the potential density. This is notably the case of the two realistic stratifications considered in Section 6 (North-Western Mediterranean and North-Eastern Atlantic). In a sigma coordinate model, this means that taking compressibility into account generates PGF errors that may be greater than the physical signal itself. Indeed the part of $δρ^C$ that depends on $z$ only (which should not create any horizontal pressure gradient) possibly leads to an erroneous PGF (i.e. if the last terms at the right hand side of (2) and (3) fail to balance the other terms) higher than the PGF theoretically expected because of $(θ, S)$ gradients. Besides, we note that several authors have proposed to isolate the part of $δρ^C$ that

depends on $z$ only and remove it from the rest of the equation of state (Dukowicz, 2001, SM03).

On the other hand, and provided that the PGF scheme is linear relative to $ρ$, it may split into two parts (as suggested by SM03 but only as theoretical possibility), namely:

$$PGF(ρ^0 + δρ^C) = PGF(ρ^0) + PGF(δρ^C) \quad (24)$$

which can be specifically treated according to the nature of their respective errors. In other words $PGF(ρ^0)$ and $PGF(δρ^C)$ can be computed using different schemes. Despite a complex dependency of $δρ^C$ to $z$, apparently expressed by (23), the vertical profile of $δρ^C$ reveals a certain linearity relative to depth (Gill, 1982, Fig. 3.2). The Standard Density-Jacobian PGF (SM03, Song, 1998), which is exact when the density profile is linear in $z$, may thus be considered to compute $PGF(δρ^C)$, when keeping the Pressure-Jacobian scheme to compute $PGF(ρ^0)$. We will however show that the performances of this possible hybrid PGF are lower than those of a revisited form of the Pressure-Jacobian scheme dealing with $δρ^C$ in a more suitable way. The latter is presented in the following section.

### 4. A discrete approach adapted to a complete equation of state

In order to reduce the PGF errors associated with $δρ^C(θ, S, z)$, we developed a method based on the Equivalent Geopotential Formulaion of the PGF described in Section 2.4. Depending on the EOS, different implementations of this method will be considered in order to rationalize the computational cost. In the following, we distinguish the case of the Mellor (1991) EOS from the EOS proposed by the Wright (1997) or McDougall et al. (2003).

#### 4.1. Mellor (1991) EOS

As far as $δρ^C(θ, S, z)$ is concerned, the accuracy problem of the PGF scheme mainly comes from the fact that a small variation in $z$ can generate a substantial change in $δρ^C$. Even if $θ$ and $S$ are constant, the variations of $δρ^C$ related to $z$ cannot be exactly balanced by (2), (3), (9), when the PGF should theoretically vanish. The idea is thus to build a scheme that would not deal with discrete values of $δρ^C$ computed with different values of $z$. This is easily achievable if we now consider the alternative geopotential formulation of the PGF. From now, we distinguish $PGF(ρ^0)$ and $PGF(δρ^C)$. At this stage the choice of the scheme used for $PGF(ρ^0)$ does not really matter. It can be one of the three schemes described in Section 2, or any other scheme provided that the separation in a $δρ^C$ dependent term and a $ρ^0$ dependent term is easily tractable. We actually consider $PGF(δρ^C)$ on which we apply the Geopotential approach. The method closely follows Section 2.4 except that we change the order of the calculus. In Section 2.4, the EOS is first computed from $(θ, S, z)$ given at tracer grid nodes and then interpolated on the common $z^*$ level. We actually propose to interpolate $(θ, S)$ on the common $z^*$ level first, and then to compute $δρ^C$ from the interpolated temperature and salinity using the common $z^*$ level. This significantly reduces the PGF errors. It notably ensures that the PGF exactly vanishes if $(θ, S)$ are constant. As far as the Equivalent Geopotential Formulation of the Density-Jacobian is concerned, this new order is illustrated by Fig. 1b (to be compared to Fig. 1a) and the interpolation of $(θ, S)$ is derived from (19) and (20), namely:

$$ \begin{align} S_{i-1,k+1/2}^{L} &= \frac{S_{i-1,k}(Z_{i-1,k+1} - Z^*) + S_{i-1,k+1}(Z^* - Z_{i-1,k})}{(Z_{i-1,k+1} - Z_{i-1,k})} \nonumber \\ S_{i,k+1/2}^{R} &= \frac{S_{i,k}(Z_{i,k+1} - Z^*) + S_{i,k+1}(Z^* - Z_{i,k})}{Z_{i,k+1} - Z_{i,k}} \tag{25} \\ \theta_{i-1,k+1/2}^{L} &= \frac{\theta_{i-1,k}(Z_{i-1,k+1} - Z^*) + \theta_{i-1,k+1}(Z^* - Z_{i-1,k})}{Z_{i-1,k+1} - Z_{i-1,k}} \nonumber \\ \theta_{i,k+1/2}^{R} &= \frac{\theta_{i,k}(Z_{i,k+1} - Z^*) + \theta_{i,k+1}(Z^* - Z_{i,k})}{Z_{i,k+1} - Z_{i,k}} \nonumber \end{align} $$
As far as the Janjic (1997) PGF or the Pressure-Jacobian PGF are concerned, the interpolation of $(\theta, S)$ is derived from (14), that is:

$$ 
\begin{aligned}
\theta_{i,k+1/2}^R &= \alpha_{i,k+1/2}^R \theta_{i,k+1} + (1 - \alpha_{i,k+1/2}^R) \theta_{i,k} \\
S_{i,k+1/2}^R &= \alpha_{i,k+1/2}^R S_{i,k+1} + (1 - \alpha_{i,k+1/2}^R) S_{i,k} \\
\theta_{i-1,k+1/2}^L &= \alpha_{i-1,k+1/2}^L \theta_{i-1,k+1} + (1 - \alpha_{i-1,k+1/2}^L) \theta_{i-1,k} \\
S_{i-1,k+1/2}^L &= \alpha_{i-1,k+1/2}^L S_{i-1,k+1} + (1 - \alpha_{i-1,k+1/2}^L) S_{i-1,k}
\end{aligned} 
\quad (26) 
$$

where $(\alpha_{i-1,k+1/2}^R, \alpha_{i,k+1/2}^R)$ are given by (15). We may note this subtle difference with the Density-Jacobian case: the common geopotential level $z_{i-1/2,k}^L$ defined by (16) or (11) corresponds to the pressure $(p^R, p^L)$ level. As the discrete integral (12) used to compute the hydrostatic pressure assumes that the interpolated temperature and salinity are located between two consecutive pressure levels, a common geopotential level for $(\theta, S)$ is deduced from two consecutive pressure levels according to:

$$ z_{i-1/2,k+1/2}^L = \frac{z_{i-1/2,k}^L + z_{i-1/2,k+1}^L}{2} \quad (27) $$

Using (27) and (26), we then obtain the two discrete densities, $\delta\rho^C(\theta_{i,k+1/2}^R, S_{i,k+1/2}^R, z_{i-1/2,k+1/2}^L)$ and $\delta\rho^C(\theta_{i-1,k+1/2}^L, S_{i-1,k+1/2}^L, z_{i-1/2,k+1/2}^R)$, involved in the computation of the hydrostatic pressure (12).

## 4.2. Computational cost

We note from Fig. 2 that $p_{i,k+1}^R$ and $p_{i,k+1}^L$ do not share the same geopotential level, leading $\delta\rho^C(\theta_{i,k+1/2}^R, S_{i,k+1/2}^R, z_{i-1/2,k+1/2}^L)$ and $\delta\rho^C(\theta_{i,k+1/2}^L, S_{i,k+1/2}^L, z_{i+1/2,k+1/2}^R)$ (yet sharing the same horizontal position) to be subject to different interpolation rules. If we consider the two horizontal directions, it finally appears that $\delta\rho^C$ needs to be computed four times per each tracer grid point. As far as the Mellor (1991) EOS is concerned, the numerical expression of $\delta\rho^C$ is rather simple and thus the computational cost of the method remains affordable. Moreover, as the numerical expression of $\delta\rho^C$ is much simpler than that of $\rho^0$, the choice of separating $\delta\rho^C$ and $\rho^0$ ($\rho^0$ being computed once per grid point at tracer locations) is particularly judicious from a computational cost point of view. The cases of the EOS of Wright (1997) and McDougall et al. (2003) are quite different and require a particular implementation described in the next section.

## 4.3. Wright (1997) EOS

Wright (1997) EOS is given by:

$$ \rho(\theta, S, p) = \frac{p + p_0}{\lambda + \alpha_0(p + p_0)} \quad (28) $$

where $p$ is the pressure and $(p_0, \lambda, \alpha_0)$ are polynomial functions of $\theta$ and $S$. Its very simple pressure dependency makes Wright (1997) EOS attractive, notably for analytic approaches (Adcroft et al., 2008). In the present study, as for the previous case, the density $\rho(\theta, S, p)$ is expressed as the sum of a potential density $\rho^0 = \rho(\theta, S, 0)$ and a density anomaly related to the compression terms $\delta\rho^C(\theta, S, p) = \rho(\theta, S, p) - \rho(\theta, S, 0)$. As far as (28) is concerned, and using the approximation $p = -z \times 10^4$ (note that the pressure in Wright 1997 EOS is given in Pascal), this leads to:

$$ \delta\rho^C = \frac{a_1 z}{b_0 + b_1 z} \quad (29) $$

with $(a_1, b_0, b_1) = (-10^4\lambda, (\lambda + \alpha_0 p_0)^2, -10^4\alpha_0(\lambda + \alpha_0 p_0))$. The comparison of expressions (23) and (29) shows that the pressure dependency of Wright (1997) EOS is simpler than that of Mellor (1991) EOS. On the other hand, the dependency to $(\theta, S)$ of (29) is significantly more complex than (23) since $(p_0, \lambda, \alpha_0)$ are given by polynomial functions of $(\theta, S)$ based on 15 coefficients (see Table 1 in

Wright 1997). The related numerical cost is moreover aggravated by the fact that our method leads to compute $\delta\rho^C$ four times per each tracer grid point. However it is possible to reduce the computational cost if the coefficients $(a_1, b_0, b_1)$ of (29) are first computed at tracer grid points and then interpolated at the common height $z^*$. Using (29), $\delta\rho^C$ is finally given by

$$ 
\begin{aligned}
\delta\rho_{i,k+1/2}^C &= \frac{a_{i1,k+1/2}^R z_{i-1/2,k+1/2}^*}{b_{0i,k+1/2}^R + b_{1i,k+1/2}^R z_{i-1/2,k+1/2}^*} \\
\delta\rho_{i-1,k+1/2}^C &= \frac{a_{i-1,k+1/2}^L z_{i-1/2,k+1/2}^*}{b_{0i-1,k+1/2}^L + b_{1i-1,k+1/2}^L z_{i-1/2,k+1/2}^*}
\end{aligned} 
\quad (30) 
$$

where the coefficients $(a_{i1,k+1/2}^R, b_{0i,k+1/2}^R, b_{1i,k+1/2}^R, a_{i-1,k+1/2}^L, b_{0i-1,k+1/2}^L, b_{1i-1,k+1/2}^L)$ are obtained as follows: first $(a_1, b_0, b_1)$ are computed at tracer grid points and are then interpolated on the common height $z_{i-1/2,k+1/2}^*$ according to:

$$ 
\begin{aligned}
a_{i1,k+1/2}^R &= \alpha_{i,k+1/2}^R a_1(\theta_{i,k+1}, S_{i,k+1}) + (1 - \alpha_{i,k+1/2}^R) a_1(\theta_{i,k}, S_{i,k}) \\
b_{0i,k+1/2}^R &= \alpha_{i,k+1/2}^R b_0(\theta_{i,k+1}, S_{i,k+1}) + (1 - \alpha_{i,k+1/2}^R) b_0(\theta_{i,k}, S_{i,k}) \\
b_{1i,k+1/2}^R &= \alpha_{i,k+1/2}^R b_1(\theta_{i,k+1}, S_{i,k+1}) + (1 - \alpha_{i,k+1/2}^R) b_1(\theta_{i,k}, S_{i,k}) \\
a_{i-1,k+1/2}^L &= \alpha_{i-1,k+1/2}^L a_1(\theta_{i-1,k+1}, S_{i-1,k+1}) + (1 - \alpha_{i-1,k+1/2}^L) a_1(\theta_{i-1,k}, S_{i-1,k}) \\
b_{0i-1,k+1/2}^L &= \alpha_{i-1,k+1/2}^L b_0(\theta_{i-1,k+1}, S_{i-1,k+1}) + (1 - \alpha_{i-1,k+1/2}^L) b_0(\theta_{i-1,k}, S_{i-1,k}) \\
b_{1i-1,k+1/2}^L &= \alpha_{i-1,k+1/2}^L b_1(\theta_{i-1,k+1}, S_{i-1,k+1}) + (1 - \alpha_{i-1,k+1/2}^L) b_1(\theta_{i-1,k}, S_{i-1,k})
\end{aligned} 
\quad (31) 
$$

The weights of the interpolation, $\alpha_{i,k+1/2}^R$ and $\alpha_{i-1,k+1/2}^L$, are given by (15). The implementation of the method is illustrated by Fig. 1c (to be compared with Fig. 1b).

## 4.4. McDougall et al. (2003) EOS

The McDougall et al. (2003) EOS is of the form:

$$ \rho = \frac{a_0 + a_1 p + a_2 p^2}{b_0 + b_1 p + b_2 p^2 + b_3 p^3} \quad (32) $$

where $p$ is the pressure in decibars and $(a_0, a_1, a_2, b_0, b_1, b_2, b_3)$ are given by polynomial functions of $(\theta, S)$ defined by 25 coefficients updated by Jackett et al. (2006). We also note that a similar formula is available for the Conservative Temperature and the Absolute Salinity (IOC SCOR and IAPSO, 2010). The pressure dependant density anomaly is given by $\delta\rho^C = \rho - a_0/b_0$. As the structure of (32) is similar to that of (28), the order of the polynomials being simply higher in the case of the McDougall et al. (2003) EOS, the method used for Wright (1997) EOS still applies here. Practically $(a_0, a_1, a_2, b_0, b_1, b_2, b_3)$ are first computed on the tracer grid points and then interpolated on the common height $z^*$, as in (31). We can note that in the special case where both PGF($\delta\rho^C$) and PGF($\rho^0$) are based on the same scheme (for instance both Pressure-Jacobian), there is no particular need to separate $\delta\rho^C$ and $\rho^0$. As far as the Wright EOS or the McDougall et al. EOS are concerned, dealing with a reunified density $\rho = \rho^0 + \delta\rho^C$ contributes to lower the computational cost of our method.

# 5. A high-order Finite-Volume method derived from Adcroft et al. (2008)

Adcroft et al. (2008) (hereafter A08) proposed a method to deal with the compressibility of sea water based on an analytical, Finite Volume (hereafter FV), approach of the PGF. Although we use a substantially different model (the A08 model is based on isopycnal layers with momentum equations making use of the Montgomery potential), this FV method can be applied here. We may note that the general principle of the FV approach of these authors follows
that of the Lin, 1997) Pressure-Jacobian PGF. Indeed the pressure force is obtained from the integral of the pressure over the edges of momentum cell boxes, namely:

$$A^{-1} \times \left( \underbrace{\int_{z_{i,k-1/2}}^{z_{i,k+1/2}} p \, dz}_{I} + \underbrace{\int_{z_{i,k+1/2}}^{z_{i-1,k+1/2}} p \, dz}_{II} + \underbrace{\int_{z_{i-1,k+1/2}}^{z_{i-1,k-1/2}} p \, dz}_{III} + \underbrace{\int_{z_{i-1,k-1/2}}^{z_{i,k-1/2}} p \, dz}_{IV} \right) \quad (33)$$

In the discrete approach of Lin (1997), the four terms in brackets and the area A are approximated by:

$$\begin{align*}
I &= \frac{p_{i,k+1/2} + p_{i,k-1/2}}{2} (z_{i,k+1/2} - z_{i,k-1/2}) \\
II &= \frac{p_{i-1,k+1/2} + p_{i,k+1/2}}{2} (z_{i-1,k+1/2} - z_{i,k+1/2}) \\
III &= \frac{p_{i-1,k-1/2} + p_{i-1,k+1/2}}{2} (z_{i-1,k-1/2} - z_{i-1,k+1/2}) \\
IV &= \frac{p_{i,k-1/2} + p_{i,k-1/2}}{2} (z_{i,k-1/2} - z_{i-1,k-1/2}) \\
A &= \Delta x \frac{\Delta z_{i,k} + \Delta z_{i-1,k}}{2}
\end{align*} \quad (34)$$

where the location of variables is shown in Fig. 5. Actually, using (34) in (33) leads to the right hand side of (6).

Alternatively, the pressure integral can be written:

$$\int_{z_1}^{z_2} p(z) dz = \underbrace{\int_{z_1}^{z_2} \left( \int_{z}^{\zeta} g \rho^0 dz' \right) dz}_{I} + \underbrace{\int_{z_1}^{z_2} \left( \int_{z}^{z_2} g \delta \rho^C dz' \right) dz}_{II} + \underbrace{\int_{z_1}^{z_2} \left( \int_{z_2}^{\zeta} g \delta \rho^C dz' \right) dz}_{III} \quad (35)$$

where $(z_1, z_2)$ represents any of the four pairs of levels in (33) and where, as previously, we distinguish a potential density and a depth-dependent density anomaly, that is $\rho = \rho^0 + \delta\rho^C$. In the A08 approach, the depth dependent part of the sea water density (e.g. terms II and III in (35)) is subjected to an analytical treatment. As the latter possibly leads to complicated expressions, the simplicity of the Wright (1997) EOS is, in these circumstances, rather attractive. Even though, the problem remains complex because of the possible variations, over the range of the vertical integral, of the $(\theta, S)$ dependent coefficients of the EOS. This point is discussed by A08 who proposed practical solutions.

As far as the pressure integral is calculated over the two vertical segments (see terms I and III in 33 and Fig. 5), the temperature and the salinity (and consequently the EOS coefficients) can notably be assumed constant within the cell box. In this case, and using (29), the term II of (35) leads to:

$$g \left| \frac{a_1 b_0}{b_1^3} (b_0 + b_1 z) (\ln(b_0 + b_1 z) - 1) - \frac{a_1}{b_1} \left( \frac{b_0}{b_1} z + \frac{z^2}{2} \right) \right|_{z_1}^{z_2} \quad (36)$$

where the symbol $|z_1^2|$ stands for the following operation: $\int f(z)|z_1^2| = f(z_2) - f(z_1)$, and where the EOS coefficient $(a_1, b_0, b_1)$ are computed using $(\theta_{i,k}, S_{i,k})$ (case $(z_1, z_2) = (z_{i,k-1/2}, z_{i,k+1/2})$) or using $(\theta_{i-1,k}, S_{i-1,k})$ (case $(z_1, z_2) = (z_{i-1,k-1/2}, z_{i-1,k+1/2})$). On the other hand, the density integral of the term III of (35) can be calculated using a method combining a hybrid discrete-analytical approach, namely:

$$\underbrace{\int_{z_1}^{z_2} \left( \int_{z_2}^{\zeta} g \delta \rho^C dz \right) dz}_{III} = \int_{z_1}^{z_2} p^C(z_2) dz = p^C(z_2) |z|_{z_1}^{z_2} \quad (37)$$

where, considering the case $z_2 = z_{i,k+1/2}$, the pressure anomaly $p^C(z_2)$ is given by:

$$p^C(z_{i,k+1/2}) = p^C(z_{i,k+3/2}) + \int_{z_{i,k+1/2}}^{z_{i,k+3/2}} g \delta \rho^C dz \quad (38)$$

The $z_{i,k_{max+1/2}}$ level corresponds to the sea surface where we have $p^C(z_{i,k_{max+1/2}}) = 0$, while the second term at the right hand side of (38) is given by an analytical calculus expression based on (29), that is:

$$\int_{z_{i,k+1/2}}^{z_{i,k+3/2}} g \delta \rho^C dz = \left| \frac{a_1}{b_1} \left\{ z + \frac{b_0}{b_1} \left[ 1 - \ln(b_0 + b_1 z) \right] \right\} \right|_{z_{i,k+1/2}}^{z_{i,k+3/2}} \quad (39)$$

where the EOS coefficients $(a_1, b_0, b_1)$, assumed constant, are computed with $(\theta_{i,k+1}, S_{i,k+1})$.

As far as the pressure integral concerns the lower and upper facets of the cell box (see terms II and IV in 33 and Fig. 5), the analytical approach is hampered by the horizontal variations of temperature and salinity. This was reported by Adcroft et al. (2008) who proposed to overcome the problem through a high-order numerical approach and some simple hypothesis on subgrid variations of the fields. Here, we assume that the $(\theta, S)$ field at a given vertical grid level varies linearly between two neighbouring grid nodes. According to A08, a sixth-order quadrature apparently provides enough accuracy. As the computational cost is not really an issue for us at this stage, we use an even higher (e.g. 10th-order) scheme. Besides, we checked that this alternative numerical approach could be extended to the four terms of (33) with a very satisfying accuracy (the PGF can notably be considered cancelled for constant $(\theta, S)$). The interest is obviously to render the computation of the pressure integral more “EOS independent” and, insofar as the Mellor (1991) and McDougall et al. (2003) EOSs are significantly more complex than the Wright (1997) EOS, to render the FV approach easily implementable whatever the complexity of the considered EOS.

# 6. Tests and discussion

## 6.1. The Seamount experiment

We now repeat the seamount experiment presented in M09. The initial profile of salinity and potential temperature is:

$$\begin{align}
S &= 35 \nonumber \\
\theta &= 20e^{z/500} \tag{40}
\end{align}$$

which roughly corresponds to the potential density profile used in SM03. As in M09, the initial discrete fields are provided by the average of the true fields over the volume of the cell boxes, i.e. $\Theta_{i,k} = 20 \times 500 \times (\exp(z_{i,k+1/2}/500) - \exp(z_{i-1,k-1/2}/500))/(z_{i,k+1/2} - z_{i,k-1/2})$ (see also expression (11) in M09). The other fields are given by the ocean at rest. The numerical model (SYMPHONIE) is derived from

**Table 1**

PGF schemes used for the Seamount Test. Left column: PGF($\rho^0$), Right column: PGF($\delta\rho^C$). We recall that FV stands for “Finite Volume” and EGF for “Equivalent Geopotential Formulation”.

<table><thead><tr><th>PGF(ρ<sup>0</sup>)</th><th>PGF(δρ<sup>C</sup>)</th></tr></thead><tbody><tr><td>Pressure-Jacobian (Lin, 1997)</td><td>Pressure-Jacobian (Lin, 1997)</td></tr><tr><td>Pressure-Jacobian (Lin, 1997)</td><td>Density-Jacobian (Song, 1998) (green in Fig. 3)</td></tr><tr><td>Pressure-Jacobian (Lin, 1997)</td><td>High Order FV scheme (blue in Fig. 3)</td></tr><tr><td>Pressure-Jacobian (Lin, 1997)</td><td>EGF Janjic PGF</td></tr><tr><td>Pressure-Jacobian (Lin, 1997)</td><td>EGF Density-Jacobian</td></tr><tr><td>Pressure-Jacobian (Lin, 1997)</td><td>EGF Pressure-Jacobian (red in Fig. 3)</td></tr><tr><td>EGF Pressure-Jacobian</td><td>EGF Pressure-Jacobian</td></tr><tr><td>Density-Jacobian (Song, 1998)</td><td>Density-Jacobian (Song, 1998)</td></tr><tr><td>Density-Jacobian (Song, 1998)</td><td>High Order FV scheme</td></tr><tr><td>Density-Jacobian (Song, 1998)</td><td>EGF Janjic PGF</td></tr><tr><td>Density-Jacobian (Song, 1998)</td><td>EGF Density-Jacobian</td></tr><tr><td>Density-Jacobian (Song, 1998)</td><td>EGF Pressure-Jacobian</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Vertical axes: depth (m), Horizontal axes: $10^7 \times \rho_0^{-1} \nabla_s p^\prime (\text{m s}^{-2})$. At the location of the profile $h = 5000 \text{ m (a)}, h = 3500 \text{ m (b)}, h = 1800 \text{ m, (c)}, h = 700 \text{ m (d)$. See also Fig. 2 in M09. Usual Density Jacobian discretization (green), Equivalent Geopotential Formulation of the Pressure-Jacobian (red), Finite Volume approach (blue).</figcaption></figure>

Marsaleix et al. (2008) and M09 but we now use the complete equation of state proposed by McDougall et al. (2003).

The ocean should remain at rest but the PGF errors are responsible for erroneous currents. We follow their long term evolution with a 180-days run similar to that presented in M09 and SM03.

We compare several low order schemes listed in Table 1. As far as the potential density is concerned, $\text{PGF}(\rho^\prime)$ is computed using the genuine formulation of the Pressure-Jacobian or the Density-Jacobian schemes (left column in Table 1). In each case, $\text{PGF}(\delta\rho^C)$ is computed with the usual genuine scheme, or using the various Equivalent Geopotential Formulations (hereafter EGF) proposed in Section 4, or the high order FV approach described in Section 5 (right column in Table 1). We also consider a case where both $\text{PGF}(\rho^\prime)$ and $\text{PGF}(\delta\rho^C)$ are computed using the Equivalent Geopotential Formulation of the Pressure-Jacobian.

Fig. 3, which can be directly compared to Fig. 3 of M09, presents initial PGF errors at different locations of the seamount (depicted in Fig. 2 of M09). In order to evidence the PGF errors caused by the compression terms only, Fig. 3 shows the part of the PGF solely due to $\delta\rho^C$ (the other part being shown in Fig. 3 of M09). Note that

since the three EGF schemes give similar results, only one of them (the EGF Pressure-Jacobian) is presented in Fig. 3 for the sake of clarity. As expected, PGF errors are weak, whatever the scheme, for small depth since compression terms become negligible. Errors are also small at the basis of the seamount because the slope of the $s$ levels is small (Fig. 3a). The largest errors are finally found when the slope of the bathymetry and the depth are both significant (see Fig. 3b and c). In this latter case, errors are rather of the same order as PGF errors related to $\rho^\prime$ shown in M09. Globally the errors of the EGF method are smaller than those of the Standard-Jacobian scheme or the high order FV approach. These three methods are anyway much more accurate than the usual Pressure-Jacobian discretization (errors not shown). In some areas the FV approach is the most accurate option (Fig. 3a), in others, it can be rather disappointing (Fig. 3d).

Fig. 4 shows the time variation of the total kinetic energy. The Pressure-Jacobian based on the full usual discretization (Fig. 4a, red curve) is particularly subject to the PGF errors related to the compressibility of the sea water. The error level is notably much higher than that of the full Density-Jacobian PGF (Fig. 4b, magenta
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Total kinetic energy (m² s⁻²) as a function of time (days). Fig. 4a: $PGF(ρ^0)$ is computed on the basis of the usual Pressure-Jacobian discretization. $PGF(δρ^C)$ is computed with the usual Pressure-Jacobian (red), the usual Density-Jacobian (magenta), the EGF formulation of the Janjic (1997) scheme (green), the high order Finite Volume method (blue). Fig. 4b: $PGF(ρ^0)$ is computed on the basis of the usual Density-Jacobian discretization. $PGF(δρ^C)$ is computed with the usual Density-Jacobian (magenta), the EGF formulation of the Janjic (1997) scheme (green), the high order Finite Volume method (blue).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Contour of the pressure integral used to compute the *u*-component of momentum.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Location of the vertical sections corresponding to Fig. 7 (red line) and Fig. 8 (blue line).</figcaption></figure>

curve). As far as the computation of $PGF(ρ^0)$ relies on the usual Pressure-Jacobian discretization, all the methods proposed to reduce the error level are efficient (Fig. 4a). The EGF approach clearly leads to the best results (green curve). The three EGF schemes, including the full EGF approach (e.g. both $PGF(ρ^0)$ and $PGF(δρ^C)$ using the EGF discretization), lead to very similar results (for clarity, only one of them is thus shown in Fig. 4a). Among these three schemes, the EGF of the Janjic PGF seems slightly better than the

two others. The PGF errors of the FV approach (Fig. 4, blue<fnref n="1" /> curve) are comparable to those obtained when the computation of $PGF(δρ^C)$ is based on the Density-Jacobian scheme (Fig. 4a, magenta curve).

<footnote n="1">¹ For interpretation of color in Figs. 1–4 and 6–8, the reader is referred to the web version of this article.</footnote>
As far as $PGF(\rho^0)$ is computed on the basis of the usual Density-Jacobian discretization, the reduction of the compressibility errors by the FV and EGF methods is less significant. Indeed, the error level of the full Density-Jacobian PGF (e.g. both $PGF(\rho^0)$ and $PGF(\delta\rho^C)$ based on the usual density-jacobian discretization) is of the same

order of the errors obtained when $PGF(\delta\rho^C)$ is computed using the FV or EGF methods. This is not really surprising insofar as the previous test (e.g. Fig. 4a) shows that computing $PGF(\delta\rho^C)$ with the usual density-Jacobian discretization significantly reduces the PGF errors of the Pressure-Jacobian. Anyway, the benefit of the

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** South-North pressure gradient divided by the Coriolis frequency (m s⁻¹) for the Mediterranean case (location given by the red line in Fig. 6). Colour resolution: 5 × 10⁻⁴ m s⁻¹. Horizontal axis: latitude (°). Vertical axis: depth (m). The PGF related to the potential density has been arbitrarily cancelled in order to evidence the effect of the compressibility terms of the equation of state. Reference PGF computed on a geopotential grid using the *Pressure-Jacobian* scheme (a), PGF computed on a sigma coordinate grid (b, c, d, e) using the usual *Pressure-Jacobian* discretization (b), the *Equivalent Geopotential Formulation* of the *Pressure-Jacobian* scheme (c), the high-order Finite-Volume scheme (d), the usual *Density-Jacobian* discretization (e). Note the particular colour resolution used in Fig. 7b.</figcaption></figure>
EGF approach (green curve) is small but indubitable. On the other hand, the interest of the high order FV approach (blue curve) is possibly doubtful.

Several reasons for the partial success of the FV method depicted by Figs. 3d and 4b could be invoked. We may for instance note that the benefit of the high order treatment of the depth dependency of the problem is possibly lost because of the lower order approach

used for several other crucial parameters. The fact that the pressure integral is computed under the approximation of a constant temperature and salinity, is notably not consistent with the exponential profile (40). On the other hand, the pressure integral is computed in our case over a cell box contour that is linearly interpolated from four vertices, when a more complex approach of the cell box geometry is possibly required by high order methods (SM03).

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Same as Fig. 7 but for the Atlantic case (location given by the blue line in Fig. 6). Colour resolution: $10^{-2}$ m s$^{-1}$. Note the particular colour resolution used in Fig. 8b.</figcaption></figure>
## 6.2. Realistic stratifications

We now present the PGF errors in the cases of realistic stratifications. The interest of these additional tests is to assess whether or not the error amplitude is smaller than the amplitude of the expected physical signal. Indeed, the ideal case of an ocean at rest (i.e. the Seamount experiment) does not permit to conclude on this question. Two regions are considered: the North Occidental Mediterranean Sea and the North Atlantic Ocean (Fig. 6). The stratification is provided by the NEMO OPA general circulation model on a geopotential grid using 72 levels in the Mediterranean case (Tonani et al., 2009) and 50 levels in the Atlantic case (Dombrowsky et al., 2009; Hurlburt et al., 2009). The horizontal resolution is respectively 1/16° and 1/12°. The (θ, S) fields are interpolated on the grid of our model (horizontal resolution: 2.5 km). Two vertical grids are used. First of all, a geopotential grid, with a vertical distribution identical to that of the NEMO OPA model, enables to compute the PGF without truncation errors related to the sigma coordinate, thus providing a reference solution (Figs. 7a and 8a). The second vertical grid is based on 40 sigma levels regularly distributed. The PGF is thus subject to the truncation errors of the sigma system. In order to evidence the errors related to the compressibility terms of the EOS, PGF(ρ°) is excluded, so that the computation of the PGF relies on δρ° only. The PGF has been divided by the Coriolis Frequency so that Figs. 7 and 8 may be interpreted in terms of geostrophic current balance. The latter is particularly weak, notably in Fig. 7a, since it hardly reaches a few cm/s. Clearly, even a small error on the PGF computed in the sigma coordinate may be sufficient to hide the physical signal. Along these lines, it is clear that the usual Pressure-Jacobian discretization should be avoided since the induced truncation errors completely hide the physical signal (Figs. 7b and 8b). On the other hand the PGF is in good agreement with the reference solution when it relies on the other schemes (e.g. Equivalent Geopotential Formulation, high order FV, density Jacobian). In these three cases, the error level is thus lower than the physical signal. Consequently, even though the previous section shows that the EGF method performs better than the two others, it does not seem unreasonable to use any of these three methods on a realistic case.

We note that the agreement with the reference solution is rather better in the Mediterranean case. Actually, the North-Western Mediterranean stratification tends to become rather homogeneous with decreasing z, placing the PGF scheme in a favourable situation since the geopotential formulation exactly cancels for (θ, S) constant. In the North-East Atlantic case, (θ, S) vertical variations are not negligible at deep levels where the compressibility effect becomes effective. As a consequence, even if the comparison of Fig. 8a, c-e is globally satisfying, we note that truncation errors tend to emerge from the physical signal in regions of strong hydrostatic inconsistency (Fig. 8, latitude = 44.2).

## 7. Conclusion

Even if the PGF scheme studied in M09 (a numerical equivalent of the L97 Pressure-Jacobian) was the starting point of the present study, we may retain that the method proposed to remove the truncation error associated to the compressibility terms of the EOS is not limited to a particular PGF scheme. Actually, this method can be applied to any scheme that can be rewritten using an Equivalent Geopotential Formulation. We actually showed that it can be easily applied to the Standard-Jacobian or the Janjic (1977) schemes.

By placing the computation of the compressibility terms on a common geopotential level z*, this method removes the errors related to the pressure dependency of the density. A major result is

indeed the exact cancellation of the PGF when (θ, S) are constant. Because our discrete approach is not really limited by the complexity of the EOS, it appears competitive compared to some analytic approaches. As a matter of fact, the Finite Volume method developed by Adcroft et al. (2008) seems rather difficult to apply to any EOS significantly more complex than the Wright (1997) EOS. A high order numerical scheme can be substituted to the analytical treatment proposed by Adcroft et al. (2008) with a reasonable level of accuracy, but the efficiency of the method remains anyway lower than that of the EGF method. This latter comment should be however tempered by the fact that for a given numerical experiment, different models can lead to different conclusions, because models are not only sensitive to the PGF scheme but also to the way the other processes are computed (Coriolis, advection, turbulence and so on). We consequently can not exclude that the FV method would better perform in other models.

On the other hand, the computation of (θ, S) at the common height z* remains subject to the errors of the interpolation scheme (26). The same remark can be done concerning the interpolation of the (aᵢ, bᵢ) coefficients of Wright or McDougall EOSs using (31). As a matter of fact, these linear interpolations become particularly inaccurate in situations of hydrostatic inconsistency (the interpolation turning into extrapolation) and thus could certainly be improved by a higher order method.

Last we note that we did not use the well-known method that consists in removing a background density profile (horizontally homogeneous) from the PGF computation. The reason is that this method can be applied to any of the schemes considered in this paper and thus does not help to differentiate the different approaches. This does not mean that this good sigma-modelling practice should be avoided. As a matter of fact, it can be combined to the EGF method. For instance, the PGF term retained in the momentum equations can be:

$$ \text{PGF} = \text{PGF}(\rho) - \text{PGF}(\rho(\theta_0(z), S_0(z), z)) \quad (41) $$

This ensures the cancellation of the PGF when (θ, S) = (θ₀(z), S₀(z)), a result that is normally not verified with a trivial approach because of truncation errors. As far as the pressure-dependent density δρ° is concerned, we recall that PGF(δρ°(θ₀, S₀, z)) vanishes anyway if (θ₀, S₀) are simple constants. Applying (41) is thus interesting provided that (θ₀, S₀) are function of z and, of course, provided that these reference profiles are reasonably representative of the temperature and salinity fields. We also thank several anonymous reviewers for their helpful comments.

## Acknowledgements

This study was supported by the French project ANR COMODO. We used the 3D ocean circulation SYMPHONIE developed by the SIROCCO oceanography group of the Toulouse University and CNRS. We also thank the Laboratoire d'Aérologie computer team, Serge Prieur, Laurent Cabanas, Jérémy Leclercq, Didier Gazen and Juan Escobar for their support.

## References

Adcroft, A., Hallberg, R., Harrison, M., 2008. A finite volume discretization of the pressure gradient force using analytic integration. Ocean Modell. 22, 106–113.

Dombrowsky, E., Bertino, L., Brassington, G.B., Chassignet, E.P., Davidson, F., Hurlburt, H.E., Kamachi, M., Lee, T., Martin, M.J., Mei, S., Tonani, M., 2009. GODAE systems in operation. Oceanography 22, 80–95.

Dukowicz, J.K., 2001. Reduction of density and pressure-gradient errors in ocean simulations. J. Phys. Oceanogr. 31, 1915–1921.

Gill, A.E., 1982. Atmosphere-Ocean Dynamics. Academic Press, 662 pp.

Hurlburt, H.E., Brassington, G.B., Drillet, Y., Kamachi, M., Benkiran, M., Bourdallé-Badie, R., Chassignet, E.P., Jacobs, G.A., Le Galloudec, O., Lellouche, J.-M., Metzger, E.J., Oke, P.R., Pugh, T.F., Schiller, A., Smedstad, O.M., Tranchant, B., Tsujino, H., Usui, N., Wallcraft, A.J., 2009. High-resolution global and basin-scale ocean analyses and forecasts. Oceanography 22, 110–127.
IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of seawater
– 2010: Calculation and use of thermodynamic properties. Intergovernmental
Oceanographic Commission, Manuals and Guides N 56, UNESCO (English), 196
pp.

Jackett, D.R., McDougall, T.J., Feistel, R., Wright, D.G., Griffies, S.M., 2006. Algorithms for density, potential temperature, conservative temperature, and the freezing temperature of seawater. J. Atmos. Ocean. Technol. 23, 1709–1728.

Janjic, Z.I., 1977. Pressure gradient force and advection scheme used for forecasting with steep and small scale topography. Contrib. Atmos. Phys. 50, 186–199.

Lin, Shian-Jiann, 1997. A finite-volume integration method for computing pressure gradient force in general vertical coordinates. Quart. J. Roy. Meteorol. Soc. 123, 1749–1762.

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

Marsaleix, P., Auclair, F., Estournel, C., 2009. Low-order pressure gradient schemes in sigma coordinate models: the seamount test revisited. Ocean Modell. 30, 169–177.

McDougall, T.J., Jackett, D.R., Wright, D.G., Feistel, R., 2003. Accurate and computationally efficient algorithms for potential temperature and density of sea water. J. Atmos. Ocean. Technol. 20, 730–741.

Mellor, G.L., 1991. An equation of state for numerical models and estuaries. J. Atmos. Ocean. Technol. 1991, 609–611.

Mellor, G.L., Ezer, T., Oey, L.-Y., 1994. The pressure gradient conundrum of sigma coordinate ocean models. J. Atmos. Ocean. Technol. 11, 1126–1134.

Mertz, G., Wright, D.G., 1992. Interpretations of the Jebar term. J. Phys. Oceanogr. 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. J. Geophys. Res. 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. Mon. Weather Rev. 126, 3213–3230.

Tonani, M., Pinardi, N., Fratianni, C., Pistoria, J., Dobricic, S., Pensieri, S., de Alfunso, M., Nittis, K., 2009. Mediterranean forecasting system: forecast and analysis assessment through skill scores. Ocean Sci. 5, 649–660.

Wright, D.G., 1997. An equation of state for use in ocean models: Eckart's formula revisited. J. Atmos. Ocean. Technol. 14, 735–740.