# Defining a Simplified Yet “Realistic” Equation of State for Seawater

Fabien Roquet, Gurvan Madec, Laurent Brodeau, J. Nycander

► To cite this version:

Fabien Roquet, Gurvan Madec, Laurent Brodeau, J. Nycander. Defining a Simplified Yet “Realistic” Equation of State for Seawater. Journal of Physical Oceanography, 2015, 45 (10), pp.2564-2579. 10.1175/JPO-D-15-0080.1 . hal-01234084

HAL Id: hal-01234084

https://hal.science/hal-01234084

Submitted on 12 Nov 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
# Defining a Simplified Yet “Realistic” Equation of State for Seawater

FABIEN ROQUET

Department of Meteorology (MISU), Stockholm University, Stockholm, Sweden

GURVAN MADEC

LOCEAN Sorbonne Universités (UPMC, University of Paris 6)/CNRS/IRD/MNHN, Paris, France

LAURENT BRODEAU AND J. NYCANDER

Department of Meteorology (MISU), Stockholm University, Stockholm, Sweden

(Manuscript received 27 April 2015, in final form 24 June 2015)

## ABSTRACT

There is a growing realization that the nonlinear nature of the equation of state has a deep impact on the global ocean circulation; however, the understanding of the global effects of these nonlinearities remains elusive. This is partly because of the complicated formulation of the seawater equation of state making it difficult to handle in theoretical studies. In this paper, a hierarchy of polynomial equations of state of increasing complexity, optimal in a least squares sense, is presented. These different simplified equations of state are then used to simulate the ocean circulation in a global 2°-resolution configuration. Comparisons between simulated ocean circulations confirm that nonlinear effects are of major importance, in particular influencing the circulation through determination of the static stability below the mixed layer, thus controlling rates of exchange between the atmosphere and the ocean interior. It is found that a simple polynomial equation of state, with a quadratic term in temperature (for cabbeling), a temperature-pressure product term (for thermobaricity), and a linear term in salinity, that is, only four tuning parameters, is enough to simulate a reasonably realistic global circulation. The best simulation is obtained when the simplified equation of state is forced to have an accurate thermal expansion coefficient near the freezing point, highlighting the importance of polar regions for the global stratification. It is argued that this simplified equation of state will be of great value for theoretical studies and pedagogical purposes.

## 1. Introduction

The equation of state (EOS) of seawater is a thermodynamic relationship whose definition derives from the first law of thermodynamics applied to a parcel of seawater, although it is determined empirically in practice. In this paper, we consider the International Thermodynamic Equation of Seawater—2010 (TEOS-10) standard (IOC et al. 2010), which defines seawater density ρ as a function of Absolute Salinity S_A (see McDougall et al.

2012), in situ temperature, and gauge pressure. The TEOS-10 variant, which is a function of Conservative Temperature Θ (see McDougall 2003) instead of in situ temperature and geopotential depth instead of gauge pressure (defined as Z ≈ p × 1 m dbar⁻¹), will be preferred here (defined in Roquet et al. 2015) because it is more suitable for ocean models, making the Boussinesq approximation (Young 2010):

$$ \rho = \rho(\Theta, S_A, Z). \qquad (1) $$

For simplicity, we will refer to Absolute Salinity as “salinity” and to Conservative Temperature as “temperature” in the following.

The EOS is a nonlinear function of salinity, temperature, and pressure, which introduces considerable difficulty in the analysis of the ocean circulation, as it then becomes impossible to define a truly conservative

Denotes Open Access content.

Corresponding author address: Fabien Roquet, Department of Meteorology (MISU), Stockholm University, S-106 91 Stockholm, Sweden.

E-mail: fabien.roquet@gmail.com

DOI: 10.1175/JPO-D-15-0080.1
<figure><img src="image_2.png" /><figcaption>FIG. 1. Climatological distribution of the thermal expansion coefficient (kg m⁻³ K⁻<fnref n="1" />). (a) Surface map and (b) Atlantic section along 30°W. (c) Potential thermal expansion along 30°W, defined as the thermal expansion that would have a seawater parcel once it has been lifted up adiabatically to the surface. Contour line: 0.05 kg m⁻³ K⁻<fnref n="1" />.</figcaption></figure>

density variable (McDougall 1987). The most critical nonlinearity of the EOS arises from the large relative variations of the thermal expansion coefficient, while, at the same time, the haline contraction coefficient remains essentially constant (relative variations ≈2% in the global ocean). Assuming a seawater Boussinesq and hydrostatic model for the ocean,<fnref n="1" /> we define the thermal expansion a and haline contraction b coefficients as follows:

$$a(\Theta, S_A, Z) = -\left.\frac{\partial \rho}{\partial \Theta}\right|_{S_A,Z}, \quad b(\Theta, S_A, Z) = \left.\frac{\partial \rho}{\partial S_A}\right|_{\Theta,Z}. \quad (2)$$

Note that these definitions have units of kilograms per cubic meter per kelvin and kilograms per cubic meter per gram of salt per kilogram, respectively. They differ from the usual non-Boussinesq definitions by a quasi-constant factor 1/ρ (e.g., IOC et al. 2010).

The thermal expansion coefficient depends to first order on temperature, an effect that is commonly referred to as cabbeling (cab). This is clearly apparent in climatological maps of thermal expansion (Fig. 1; see also Table 1), with the largest values found in the warm thermocline (up to 0.32 kg m⁻³ K⁻<fnref n="1" />) and near-zero values in freezing cold polar regions. The thermal

<footnote n="1">¹ The Boussinesq hydrostatic model is preferred here because it is dynamically equivalent to the compressible model in the limit of vanishing perturbations of density while simplifying drastically all derivations (Young 2010; Roquet 2013). See also the discussion in section 4.</footnote>
<table><caption>TABLE 1. Mean properties of seawater in the global ocean. Polar regions are taken poleward of 60° of latitude (north and south). a and b are the Boussinesq form of thermal expansion and haline contraction, respectively. The buoyancy frequency N is given in cycles h⁻¹. The ||∇ₕρ||-weighted global values roughly corresponds to the average thermocline values, where the thermal wind is largest. Thermal and haline contributions to the HDG are given separately.</caption><thead><tr><th></th><th>Θ (°C)</th><th>S<sub>A</sub> (g kg<sup>-1</sup>)</th><th>Depth (m)</th><th>a (kg m<sup>-3</sup> K<sup>-1</sup>)</th><th>b [kg m<sup>-3</sup> (g kg<sup>-1</sup>)<sup>-1</sup>]</th><th>N (cycles h<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Global</td><td>3.90 ± 4.46</td><td>34.89 ± 0.35</td><td>2055 ± 1322</td><td>0.16 ± 0.04</td><td>0.77 ± 0.01</td><td>1.09 ± 1.28</td></tr><tr><td>Global (||∇ₕρ|| weighted)</td><td>9.38 ± 8.10</td><td>34.93 ± 0.80</td><td>789 ± 975</td><td>0.18 ± 0.07</td><td>0.77 ± 0.02</td><td>2.96 ± 2.81</td></tr><tr><td>Surface</td><td>18.15 ± 9.80</td><td>34.88 ± 1.28</td><td>0</td><td>0.24 ± 0.10</td><td>0.76 ± 0.02</td><td>3.28 ± 2.61</td></tr><tr><td>Surface, polar regions</td><td>0.03 ± 2.26</td><td>33.25 ± 2.04</td><td>0</td><td>0.05 ± 0.03</td><td>0.80 ± 0.01</td><td>5.46 ± 6.51</td></tr><tr><td>HDG (×10<sup>-6</sup> kg m<sup>-4</sup>)</td><td>||a∇ₕΘ||</td><td>||b∇ₕS<sub>A</sub>||</td><td>||∇ₕρ||</td><td></td><td></td><td></td></tr><tr><td>Global</td><td>0.14 ± 0.32</td><td>0.08 ± 0.19</td><td>0.15 ± 0.32</td><td></td><td></td><td></td></tr><tr><td>Global (||∇ₕρ|| weighted)</td><td>0.70 ± 0.79</td><td>0.37 ± 0.71</td><td>0.79 ± 0.91</td><td></td><td></td><td></td></tr><tr><td>Surface</td><td>0.83 ± 0.63</td><td>0.67 ± 1.02</td><td>1.09 ± 1.02</td><td></td><td></td><td></td></tr><tr><td>Surface, polar regions</td><td>0.20 ± 0.29</td><td>1.37 ± 2.61</td><td>1.42 ± 2.58</td><td></td><td></td><td></td></tr></tbody></table>

expansion coefficient can occasionally be negative when the freezing water is sufficiently fresh, like in parts of the Arctic Ocean or the Baltic Sea.

The thermal expansion coefficient also depends critically on depth, through the thermobaric effect. This is evident when comparing the thermal expansion values (e.g., Fig. 1b) with the “potential thermal expansion” (Fig. 1c), defined as the thermal expansion coefficient of a water parcel lifted adiabatically to the surface. Thermobaricity (therm) gives a substantial increase of thermal expansion with depth, which acts to confine regions of low thermal expansion (less than 0.1 kg m⁻³ K⁻¹) to the upper 1500 m of the water column.

Implications of the nonlinearity of the EOS on the ocean circulation are well understood at the local scale. There are primarily four dynamical effects associated with the nonlinear EOS, through

1) the definition of the horizontal density gradient and its role on the hydrostatic pressure gradient through the “thermal wind balance”

$$ \nabla_z \rho = -a \nabla_z \Theta + b \nabla_z S_A, \qquad (3) $$

where $\nabla_z$ is the horizontal gradient operator (i.e., computed at fixed depth z);

2) the definition of the squared buoyancy frequency, which gives the strength of the local stratification, a major control on the mixed layer depth and on the interior vertical mixing

$$ N^2 = g \left( a \frac{\partial \Theta}{\partial z} - b \frac{\partial S_A}{\partial z} \right); \qquad (4) $$

3) the slope of the neutral tangent plane, which defines the local “horizontal” (or lateral) direction of mixing, believed to also be the preferred direction for

advection in the quasi-adiabatic interior (e.g., McDougall and Jackett 2005)

$$ s = \frac{a \nabla_z \Theta - b \nabla_z S_A}{a \partial_z \Theta - b \partial_z S_A} = - \frac{g \nabla_z \rho}{N^2}; \quad \text{and} \quad (5) $$

4) the creation of local sources and sinks of density in the interior by isoneutral mixing (e.g., Iudicone et al. 2008; Klocker and McDougall 2010) or dianeut mixing (e.g., Killworth 1983). Cabbeling always acts to increase local density when temperature and salinity properties are mixed along a neutral plane, while thermobaricity can either increase or decrease the local density.

On the global scale though, our understanding of the nature and importance of the effects induced by the nonlinear EOS on the ocean circulation remain elusive, and it still remains to build a general theory of ocean circulation that rely on a nonlinear EOS. The importance of having a proper representation of the EOS has long been recognized in the isopycnal coordinate modeling community (Griffies 2004, p. 127) as the difficulty in defining a suitable density coordinate is a critical issue for this type of model (e.g., Sun et al. 1999; Adcroft et al. 2008). Early simulations used potential density referenced at the surface as the density coordinate, which led to a very unrealistic stratification of the simulated ocean—the simulated Antarctic Bottom Water was found on top of the North Atlantic Deep Water.

This paper aims at proposing a hierarchy of simplified EOS that are optimal in a least squares sense to be defined later and to show how realistic simulations of ocean circulation can be when these simplified EOS are used in an ocean general circulation model (OGCM) instead of the TEOS-10 relation. The comparison between simulated circulations will allow a better grasp on global effects of the nonlinearities of the EOS. Two major conclusions will be drawn from this study: 1) the
nonlinear nature of the EOS represents a major constraint on the general circulation of the ocean, mainly through its influence on the local stratification. A linearized EOS leads to a simulated ocean (with realistic surface forcing fields) that is 2°C colder than the real ocean; and 2) a simplistic linear representation of both cabbeling and thermobaricity effects (as first proposed by Vallis 2006) is enough to recover a realistic ocean simulation, paving the way for a proper treatment of the EOS in theoretical studies.

In a companion paper (Nycander et al. 2015, manuscript submitted to Geophys. Res. Lett.), it is shown that cabbeling plays a crucial role for the Antarctic Intermediate Water formation, while thermobaricity is more important for the Antarctic Bottom Water, using a similar range of model simulations with varying equations of state.

We want to warn the reader that our intention is not, in any way, to propose an “improved” equation of state to be used for realistic ocean simulations. We believe that the set of simplified equations of state that we will present here can be very useful in studies that try to gain understanding on the role of the equation of state in the setting of the large-scale circulation, as the different nonlinear effects can be readily distinguished unlike for realistic equations of state. Yet, understanding and simulating a geophysical system are two distinct goals (Held 2005), and an accurate EOS based on either the 1980 equation of state (EOS-80) or TEOS-10 standards should be used anytime a simulation of the real ocean circulation is attempted (Roquet et al. 2015). It is to avoid this kind of confusion that the word “realistic” has been placed inside quotation marks in the title of the paper.

The hierarchy of simplified EOS and the ocean model will be presented in section 2. The accuracy of the different simplified EOSs and model runs will be discussed in section 3. As a conclusion, an optimal form of simplified, yet “realistic,” seawater equations of state with four tunable parameters only will be proposed and discussed in section 4.

## 2. Method

### a. Defining a hierarchy of simplified EOS

A hierarchy of simplified EOS can be constructed systematically by defining a set of polynomial EOS functions of Absolute Salinity, Conservative Temperature, and depth as follows:

$$ \rho^{\text{poly}} = \sum_{ijk} R_{ijk} S_A^i \Theta^j Z^k. \qquad (6) $$

Indexes i, j, and k vary from 0 to a maximum value ($N_i$, $N_j$, and $N_k$, respectively) that defines the order of

the simplified EOS. The associated thermal expansion and haline contraction functional can easily be derived analytically:

$$ a^{\text{poly}} = -\sum_{ijk} i R_{ijk} S_A^i \Theta^{j-1} Z^k, \quad \text{and} \qquad (7) $$

$$ b^{\text{poly}} = \sum_{ijk} i R_{ijk} S_A^i \Theta^j Z^k. \qquad (8) $$

All polynomial terms independent of both salinity and temperature in Eq. (6) (i.e., the $R_{00k}Z^k$ terms) have no effect on the simulated circulation and tracer distribution, as they have no contribution to the horizontal density gradient. Consequently, these terms are absent from both a and b definitions, while all dynamical effects of the equation of state in a Boussinesq fluid are carried through these values (as argued in the introduction). In reality, they can still have a small effect in a non-z-coordinate numerical model due to discretization errors, but it is safe to assume that these effects remain sufficiently small to leave our conclusions unaffected.

Note that the TEOS-10 approximations proposed by Roquet et al. (2015) are also trivariate polynomials, but they differ from Eq. (6) as they take a monotonic function of Absolute Salinity as an argument instead of the Absolute Salinity itself.

Separating the depth-only dependent terms from the other, we define a density anomaly variable $\rho'$ as follows:

$$ \rho = \bar{\rho}(z) + \rho', \qquad (9) $$

with $\bar{\rho} = \sum_k R_{00k}Z^k$. The simplified EOSs considered in the following will be a representation of $\rho'$ only, as the focus is on dynamical effects.

Two nonlinear terms will be of particular importance in the following, as they are able to represent the bulk of nonlinear variations of density. The first nonlinear term varies quadratically with temperature $R_{020}\Theta^2$ inducing a linear dependence in temperature for thermal expansion. This term is the simplest possible representation of the cabbeling effect. The second nonlinear term of importance involves a temperature-depth product $R_{011}\Theta Z$ and thus a linear dependence in depth for thermal expansion that induces the thermobaricity effect. Note that the simplified equation of state proposed by Vallis (2006) incorporates the same two nonlinear terms.

### b. Determination of simplified EOS constants

A hierarchy of polynomial EOS has then been determined by minimizing in a least squares sense the global ocean error made on horizontal density gradients when $a^{\text{poly}}$ and $b^{\text{poly}}$ [see Eqs. (7) and (8)] were used instead of their true oceanic values. More specifically,
<table><caption>TABLE 2. RMS error statistics on thermal expansion *a*, haline contraction *b*, and on the horizontal density gradients, separating the thermal contribution *a*∇<sub>*h*</sub>Θ, the haline contribution *b*∇<sub>*h*</sub>S<sub>*A*</sub>, and the resulting total error on ∇<sub>*h*</sub>ρ. Relative errors normalized by their respective global-mean reference value (see Table 1) are also provided within parentheses. Coefficients for each of the simplified EOSs are given in Table 3.</caption><thead><tr><th>EOS</th><th>a (10<sup>-3</sup> kg m<sup>-3</sup> K<sup>-1</sup>)</th><th>b [10<sup>-3</sup> kg m<sup>-3</sup> (g kg<sup>-1</sup>)<sup>-1</sup>]</th><th>a∇<sub>h</sub>Θ (10<sup>-9</sup> kg m<sup>-4</sup>)</th><th>b∇<sub>h</sub>S<sub>A</sub> (10<sup>-9</sup> kg m<sup>-4</sup>)</th><th>∇<sub>h</sub>ρ (10<sup>-9</sup> kg m<sup>-4</sup>)</th></tr></thead><tbody><tr><td>lin</td><td>42.5 (27.3%)</td><td>10.0 (1.3%)</td><td>75.6 (54.2%)</td><td>3.6 (4.5%)</td><td>76.6 (50.5%)</td></tr><tr><td>cab</td><td>48.0 (30.8%)</td><td>10.0 (1.3%)</td><td>15.8 (11.3%)</td><td>3.6 (4.5%)</td><td>16.6 (10.9%)</td></tr><tr><td>cab–therm</td><td>3.0 (1.9%)</td><td>10.0 (1.3%)</td><td>5.7 (4.1%)</td><td>3.6 (4.5%)</td><td>7.0 (4.6%)</td></tr><tr><td>freez</td><td>13.3 (8.5%)</td><td>9.9 (1.3%)</td><td>10.0 (7.1%)</td><td>3.6 (4.5%)</td><td>10.3 (6.8%)</td></tr><tr><td>2order</td><td>2.9 (1.9%)</td><td>1.4 (0.2%)</td><td>5.5 (3.9%)</td><td>0.3 (0.4%)</td><td>5.6 (3.7%)</td></tr></tbody></table>

the following cost function has been minimized using the MATLAB optimization toolbox:

$$ \Pi = \frac{1}{V} \sum_{i \in \{1,N\}} v_i \{(D_{\Theta,i}^x)^2 + (D_{S_{A,i}}^x)^2 + (D_{\Theta,i}^y)^2 + (D_{S_{A,i}}^y)^2\}, \quad (10) $$

where *N* is the total number of grid points and *V* = Σ<sub>*N*</sub>v<sub>*i*</sub> is the total volume of the ocean domain, with v<sub>*i*</sub> as the volume of the *i*th grid point. The equations D<sup>x</sup><sub>Θ,i</sub> = (a<sup>poly</sup><sub>*i*</sub> - a<sup>ref</sup><sub>*i*</sub>)(δ<sub>*x*</sub>Θ)<sub>*i*</sub> and D<sup>x</sup><sub>S<sub>A</sub>,*i*</sub> = (b<sup>poly</sup><sub>*i*</sub> - b<sup>ref</sup><sub>*i*</sub>)(δ<sub>*x*</sub>S<sub>A</sub>)<sub>*i*</sub> are the temperature and salinity contributions to the zonal density gradient at the *i*th grid cell, respectively. Similarly, D<sup>y</sup><sub>Θ,i</sub> and D<sup>y</sup><sub>S<sub>A</sub>,*i*</sub> are the meridional contributions. Defining the cost function based on horizontal density gradients has the advantage of objectively scaling the relative importance of temperature and salinity gradients.

We used a climatology of ocean hydrographic properties to obtain a realistic distribution of horizontal gradients of temperature and salinity as well as the thermal expansion and haline contraction coefficients. The climatological fields of reference were taken from the Polar Science Center Hydrographic Climatology (PHC3.0; available online at http://psc.apl.washington.edu; Steele et al. 2001), a product obtained by refining the World Ocean Atlas (Locarnini et al. 2013) in the Arctic region. Climatological fields were remapped on an approximately 2° resolution, tripolar, mesh (the so-called ORCA2 mesh; Timmermann et al. 2005), which is nearly isotropic and avoids the North Pole singularity by moving the mesh poles to land points (Madec and Imbard 1996). The distribution of thermodynamic coefficients a<sup>ref</sup><sub>*i*</sub> and b<sup>ref</sup><sub>*i*</sub> was determined using the 55-term Boussinesq polynomial approximation of the TEOS-10 equation of state (Roquet et al. 2015), while (δ<sub>*x*</sub>Θ, δ<sub>*y*</sub>Θ)<sub>*i*</sub> and (δ<sub>*x*</sub>S<sub>A</sub>, δ<sub>*y*</sub>S<sub>A</sub>)<sub>*i*</sub> were estimated from discretized gradients of Θ and S<sub>A</sub>, respectively, on the ORCA2 grid.

The remapping of climatological fields was done mainly for practical reasons, as numerical simulations were performed on the same grid this facilitated subsequent

analysis. Note also that the optimization is nearly insensitive to the particular choice of discretization grid owing to the small number of parameters in the simplified EOS under consideration.

An alternative option for determining simplified EOS constants, which could be interesting in theory, is to minimize the error made on the squared buoyancy frequency [see Eq. (4); proportional to the vertical potential density gradient instead of the horizontal gradient] in the cost function. However, in practice, it gives highly unsatisfactory results so it was not pursued further. The reason is that the stratification is extremely inhomogeneous in the ocean, so the least squares minimization became entirely controlled by properties in the limited area right below the subtropical mixed layer and produces large residual errors virtually everywhere else.

The accuracy of the different simplified EOS can be quantified by the root-mean-square error (RMSE) made on climatological horizontal density gradients (see Table 2):

$$ \text{RMSE} = \sqrt{\frac{1}{V} \sum_{i} v_i \| (\nabla_z \rho)_i^{\text{poly}} - (\nabla_z \rho)_i^{\text{ref}} \|^2}, \quad (11) $$

which depends to the first order on the error made on thermal expansion values, following Eq. (3). Note that the cost function definition is similar to, but is not exactly, the squared RMSE. We preferred to separate thermal and haline contributions to the horizontal density gradient in the cost function to avoid possible compensation effects between thermal expansion and haline contraction coefficients during the optimization.

In this study, we will focus on terms in the simplified EOS of order two or less, that is, for which *i* + *j* + *k* ≤ 2 in Eq. (6). It is clear that using higher-order terms would improve the accuracy of the considered EOS, but at the expense of the physical insight that can be obtained from their study. Hence, only seven different terms will be considered: R<sub>100</sub>, R<sub>010</sub>, R<sub>110</sub>, R<sub>200</sub>, R<sub>020</sub>, R<sub>101</sub>, and R<sub>011</sub>. Values from these different parameters are compiled in Table 3.
<table><caption>TABLE 3. Coefficients for the different polynomial approximations to the equation of state.</caption><thead><tr><th></th><th>R<sub>010</sub></th><th>R<sub>100</sub></th><th>R<sub>020</sub></th><th>R<sub>011</sub></th><th>R<sub>200</sub></th><th>R<sub>101</sub></th><th>R<sub>110</sub></th></tr></thead><tbody><tr><td>lin</td><td>-1.775 × 10<sup>-1</sup></td><td>7.718 × 10<sup>-1</sup></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>cab</td><td>-8.44 × 10<sup>-2</sup></td><td>7.718 × 10<sup>-1</sup></td><td>-4.561 × 10<sup>-3</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>cab–therm</td><td>-6.51 × 10<sup>-2</sup></td><td>7.718 × 10<sup>-1</sup></td><td>-5.027 × 10<sup>-3</sup></td><td>-2.5681 × 10<sup>-5</sup></td><td></td><td></td><td></td></tr><tr><td>freez</td><td>-4.91 × 10<sup>-2</sup></td><td>7.718 × 10<sup>-1</sup></td><td>-5.539 × 10<sup>-3</sup></td><td>-3.4977 × 10<sup>-5</sup></td><td></td><td></td><td></td></tr><tr><td>2order</td><td>1.82 × 10<sup>-2</sup></td><td>8.078 × 10<sup>-1</sup></td><td>-4.937 × 10<sup>-3</sup></td><td>-2.4677 × 10<sup>-5</sup></td><td>-1.115 × 10<sup>-4</sup></td><td>-8.241 × 10<sup>-6</sup></td><td>-2.446 × 10<sup>-3</sup></td></tr></tbody></table>

### c. Model runs

Numerical simulations are performed using the ocean model Nucleus for European Modeling of the Ocean (NEMOv3.6; Madec et al. 2014), using the ORCA2-Louvain-la-Neuve Sea Ice Model (LIM) configuration that runs on a global domain with a 2° mean resolution (Timmermann et al. 2005; Mignot et al. 2013). The Louvain-la-Neuve sea ice model LIM2 (Fichefet and Maqueda 1997) is coupled to the ocean component. The configuration uses 31 z-coordinate levels with partial cells at the bottom and a filtered linear free surface (Roulet and Madec 2000). Horizontal mixing is evaluated along the isoneutral tangent plane, and the mean effect of eddy variability on tracer is represented with an eddy-induced velocity parameterization, following Gent et al. (1995). A background vertical diffusivity of $1.2 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$ was applied. In case of static instability, vertical diffusivity is enhanced to mimic convective processes. A turbulent eddy kinetic energy (TKE) scheme is used for the mixing effects of surface turbulence, following Blanke and Delecluse (1993). Also, a bottom boundary layer parameterization mimics deep overflows, double diffusion is parameterized, and a geothermal heating flux is applied at the ocean bottom.

The simulations were forced with the normal year forcing (NYF; Griffies et al. 2009), which has been constructed based on the NCEP-NCAR interannual forcing fields (Large and Yeager 2009) to retain synoptic variability (i.e., atmospheric storms), with a seamless transition from 31 December to 1 January. Forcings are applied on the ocean surface following the CORE bulk formulas (Large and Yeager 2004). A restoration to climatological surface salinities is also added in the form of freshwater flux to prevent the model from drifting too dramatically. The interannual variability in our simulations is thus necessarily related to internal variability in the coupled ocean-sea ice system, such as deep property drifts or abrupt changes in convective patterns in ice-covered areas.

The model uses the seawater Boussinesq and hydrostatic approximations. It incorporates a new implementation of the equation of state and associated thermodynamic potentials (Roquet et al. 2015), based on the TEOS-10 international standard (IOC et al.

2010). In practice, it means that the model equations are derived for Conservative Temperature and Absolute Salinity (instead of potential temperature and practical salinity, respectively). Moreover, the new thermodynamic module has been coded in a way that makes it straightforward to modify the equation of state in a consistent way. In particular, four forms of simplified EOS have been implemented (see Table 3 for coefficient values):

1) the linear EOS

$$ \rho^{\text{lin}} = R_{100}S_A + R_{010}\Theta; \quad (12) $$

2) the cabbeling EOS

$$ \rho^{\text{cab}} = R'_{100}S_A + R'_{010}\Theta + R'_{020}\Theta^2; \quad (13) $$

3) the cabbeling-thermobaricity EOS

$$ \rho^{\text{cab-therm}} = R''_{100}S_A + R''_{010}\Theta + R''_{020}\Theta^2 + R''_{011}\Theta Z; \quad (14) $$

4) the second-order EOS

$$ \rho^{2\text{order}} = R'''_{100}S_A + R'''_{010}\Theta + R'''_{020}\Theta^2 + R'''_{011}\Theta Z + R''_{200}S_A^2 + R''_{011}S_A Z + R''_{110}S_A \Theta. \quad (15) $$

A fifth simplified EOS has also been tested to investigate the sensitivity of the global circulation to the accuracy of the EOS in polar regions (i.e., near the freezing point). The so-called freezing EOS has the same form as the cabbeling-thermobaricity EOS but uses the additional constraint of a realistic value of the thermal expansion near the freezing point (thus reducing the number of optimized parameters to 3):

5) the freezing EOS

$$ \rho^{\text{freez}} = R'''_{100}S_A + R'''_{010}\Theta + R'''_{020}\Theta^2 + R'''_{011}\Theta Z, \quad (16) $$

with $R'''_{010} = -a^{\text{freez}} - 2\Theta^{\text{freez}} R'''_{020}$, using the numerical values $a^{\text{freez}} = 0.028 \text{ kg m}^{-3} \text{ K}^{-1}$ and $\Theta^{\text{freez}} = -1.9^{\circ}\text{C}$.

Each model run was initialized with the PHC3.0 climatological fields (Steele et al. 2001). They were then run for 1000 yr until a pseudosteady state was reached. The analyses are based on comparisons of the model
outputs of the year 1000 of each of the simulations that
differed only by their equation of state.

**3. Results**

a. Accuracy of polynomial EOS

The RMSE obtained when using the linear EOS is large, representing about 50% of the global-mean horizontal density gradient (HDG) norm (see Table 2). Generally speaking, density gradients are too large in cold regions because the thermal expansion *a* of the linear EOS is too large there and vice versa in warm regions. Adding the cabbeling term reduces the RMSE on the HDG by a factor of 5. Adding the thermobaricity term on top of the cabbeling term further reduces the RMSE by a further factor of 2. Hence, the addition of these two simple terms is enough to reduce by one order of magnitude the RMSE made on the HDG. Interestingly, adding the thermobaricity term alone, without the cabbeling term, does not improve the accuracy of the simplified EOS, indicating that the two terms are not playing a symmetrical role. Adding salinity-dependent nonlinear terms only marginally improves the overall accuracy of the simplified EOS, albeit it has a major effect on the accuracy of the haline contribution on the HDG. Higher-order temperature-dependent terms would be required to continue reducing the RMSE on the HDG. However, we will focus in this paper on simplified EOS with the thermal expansion and haline contraction coefficients having only linear dependencies, that is, the simplest possible expression of EOS nonlinearities.

The spatial distribution of thermal expansion errors is shown for three of the simplified EOS, both at the surface and across the Atlantic basin (section along 30°W longitude) in Fig. 2. For the linear EOS, the thermal expansion is simply constant, approximately equal to the global-mean thermal expansion value. The error distributions are not shown because they are simply offset versions of Figs. 1a and 1b. In this case, the error is minimal around the subtropical fronts, that is, around ±40° of latitude, for temperatures around 10°C (see Fig. 3, black contour). The error on the thermal expansion is largely reduced everywhere near the surface with the cabbeling EOS, although values remain significantly too large in polar regions. However, the error is made larger at depth than in the case of the linear EOS. This error increases nearly linearly with depth, reaching values as large as −0.1 kg m⁻³ K⁻¹ at the ocean bottom.

Adding the thermobaricity term dramatically im-
proves the situation not only at depth but also near the

surface. The error on the thermal expansion is reduced almost everywhere, except in the tropical regions near the surface, where thermal expansion values are slightly more overestimated. A band of negative errors follows the bottom of the thermocline, while positive errors are found both in tropical and polar regions. The addition of the salinity-dependent terms (not shown) does not much change this picture, except that errors in low salinity regions are significantly reduced, which includes the equatorial band, the Arctic Ocean, and the Baltic Sea. These improvements might be of large regional importance, but they are unlikely to be of primary importance at the global scale.

The freezing EOS has been forced to have a realistic thermal expansion value for surface freezing waters (Fig. 3). This added constraint prevents the simplified EOS from performing as well as the cabbeling-thermobaricity EOS on average (Table 2). This led to comparatively larger errors on the thermal expansion at great depth or in the subtropical near-surface regions but nearly no error in the convectively active Southern Ocean subpolar region (see Fig. 2).

b. *Analysis of ORCA2 runs*

The level of realism of the simulated ocean circulation obtained when the different simplified EOS are used will now be discussed, and it is thought that this should help the understanding of how and where nonlinearities of the EOS matter. Differences relative to the reference run, which uses the TEOS-10 equation of state, are presented for a few key property distributions obtained once a near-steady state was reached (i.e., after a 1000-yr-long integration).

1) OCEAN CIRCULATION

The distribution of the mean sea surface height (SSH) is shown in Fig. 4a for the reference run. It features classical large-scale patterns of the upper-ocean circulation, including the subtropical and subpolar gyre systems, and the zonally oriented Antarctic Circumpolar Current (ACC) system, seen as a ~2-m SSH meridional drop in the Southern Ocean. The ~0.5-m sea level difference between the Atlantic and Indo-Pacific basins is also observed consistently with observations, primarily related to differences in salt content.

Although the overall SSH pattern is not funda-
mentally different for the run using a linear EOS (see
Fig. 4b), the magnitude of SSH variations is signifi-
cantly larger, implying in particular ~60% stronger
surface geostrophic currents for the ACC. The SSH
distribution is dramatically improved when the cabb-
eling term is added, and SSH differences with the
reference run rarely exceed 0.2 m. Further adding the
<figure><img src="image_3.png" /><figcaption>FIG. 2. Residual error on the thermal expansion coefficient for (upper) the cabbeling EOS, (middle) cabbeling-thermobaricity EOS, and (lower) freezing EOS, showing for each simplified EOS (left) the surface map and (right) the Atlantic section along 30°W. Contour line: 0.02 kg m⁻³ K⁻¹. The linear EOS is not shown as it has a constant thermal expansion value, so the residual error is an offset version of Fig. 1. Also, the 2order EOS is not shown because it was found to be very similar to the cabbeling-thermobaricity EOS case, except some marginal improvements in low salinity areas such as the Arctic or the tropical Pacific.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

FIG. 3. Position in temperature–salinity (TS) space of the con-
tour where the error on surface thermal expansion coefficient is
zero (i.e., where $a^{\text{poly}} = a^{\text{ref}}$) for each simplified EOS. The freezing
line (thin solid blue line), corresponding to the temperature at
which seawater freezes, and contours of surface thermal expansion
coefficient (thin dash color lines) are superimposed.

thermobaricity term does not significantly improve
the SSH distribution, although some local differences
can be observed most clearly along the path of the
ACC (i.e., in the latitude band 60°–40°S).

Perhaps, not surprisingly, the addition of the ther-
mobaricity term has a much stronger impact on the
deep circulation. This is best viewed from the distri-
bution of Atlantic meridional overturning circulation
(AMOC, see Fig. 5), defined as the zonally integrated
streamfunction in the Atlantic basin for latitudes north
of the southern tip of Africa at 30°S. The reference
model AMOC is made of two overturning cells, with
a mean separation depth of 2600 m. The shallow cell
in the reference run has a 9-Sverdrup (Sv; 1 Sv ≡
10⁶ m³ s⁻¹) maximum transport, and the bottom cell
has a 2-Sv transport, which is somewhat weaker in both
cases than what observations tend to indicate. The
AMOC simulated with the linear EOS is qualitatively
similar to the reference one but with large quantitative
differences, mainly as a result of a ~500 m upward
shift of the southward return branch of the shallow
overturning cell.

The AMOC simulated with the cabbeling-only EOS is
significantly better, biased toward positive values (i.e.,
with a stronger shallow cell and a very weak bottom
cell), while the AMOC simulated with both cabbeling
and thermobaricity terms is biased toward negative

values (i.e., with a weaker shallow cell and a rather
strong bottom cell). These two simulations have rather
different AMOC, yet none is clearly closer to the ref-
erence AMOC than the other. We will return to this
point later, when presenting and interpreting results
based on the freezing EOS.

2) DEEP CONVECTION AND TEMPERATURE DISTRIBUTION

Discrepancies between the AMOC of the different
runs are related to dramatic changes in the depth and
spatial extent of convective regions both in the Southern
Ocean and in the North Atlantic polar regions. This can
be clearly seen from the distribution of maximum mixed
layer depth (MLD; see Fig. 6). Here, the maximum
MLD is estimated at each grid point as the maximum
value of monthly mean MLD.

The reference run compares reasonably well with
observations, although it has some significant convective
biases, as in the Labrador Sea where convection is ab-
sent. Note that errors in the distribution of deep con-
vection between the different possible sites in the North
Atlantic are (unfortunately) very common in ocean
simulations. With the linear EOS, deep convection be-
comes widespread in the Southern Ocean, and large
changes in the position of convection sites also appear in
the North Atlantic. As a result, the global-mean tem-
perature is dramatically reduced by nearly 2°C (from
3.8° to 1.9°C; see Table 4 and Fig. 7b).

Most of this convective overactivity is fixed by adding
the cabbeling term in the EOS. There is still slightly too
much convection happening in the Weddell Sea and in
the North Atlantic, resulting in a global-mean temper-
ature of 3.5°C, that is, 0.3°C colder than in the reference
run. However, the temperature difference with the ref-
erence run is not spatially homogeneous. Large cold
anomalies are confined in the Southern Ocean (Fig. 7c),
while warm anomalies are observed elsewhere. This
distribution reflects a barotropization of the ACC that
can also be seen in the SSH pattern (Fig. 4c) that follows
more closely the bottom topography along the ACC
pathway.

As for the AMOC patterns, the cabbeling–thermobaricity
EOS does not clearly perform better than the cab-
beling EOS in terms of MLD or temperature distribu-
tion. It features a similar region of abnormally deep
mixed layer in the Weddell Sea and an associated
–0.5°C global-mean temperature difference with the
reference. Yet, the addition of the thermobaricity
term does improve the bottom stratification, allowing
the cold tongue of Antarctic Bottom Water to slip
northward at great depth (as indicated by the reduced
horizontal gradient of temperature anomalies below
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Distribution of SSH for the reference run based on the TEOS-10 equation of state (solid contour interval: 50 cm, color contour: 10 cm). The anomaly relative to the reference distribution is then plotted for the runs based on (b) the linear EOS, (c) the cabbeling EOS, and (d) the cabbeling-thermobaricity EOS (solid contour interval: 25 cm, color contour: 5 cm).</figcaption></figure>

2000-m depth in Fig. 7d compared to Fig. 7c). This in turn improves the SSH distribution in the Southern Ocean, an important yet unexpected effect of thermobaricity.

Because the linear free-surface implementation is unable to conserve tracers exactly (Roullet and Madec 2000), global-mean salinity is drifting in our simulations by about −0.2 g kg⁻¹ during the transitory regime. This happens despite the applied sea surface salinity (SSS) restoration used here to avoid any dramatic drift of the model compared to observations that could be induced by misrepresentation of surface freshwater fluxes resulting from forcings and ocean-ice interaction. It is clearly a limitation of this study, and it could be indirectly responsible for the observed −0.15°C discrepancy between the global-mean temperature of the reference run and the observed global-mean temperature (Table 4). The use of a nonlinear free surface could help reducing this issue, yet only a coupled ocean-atmosphere

simulation could fully grasp the climatic implications of changing the EOS.

### 3) THE FREEZING EOS SIMULATION

So far, this comparative study has highlighted the prominent importance of deep convection processes happening in highly localized ocean sites through its major impact on the large-scale stratification. This suggests that it is critical to consider a realistic distribution of thermal expansion coefficients in the model especially for surface water masses near the freezing point temperature. To test this idea, an additional simplified EOS has been designed and tested using the ORCA2-LIM configuration, with parameters constrained in a way that enforced a minimized error on the surface thermal expansion coefficients near the freezing temperature for salinities around 35 g kg⁻¹ (see Fig. 3).

Errors on mean SSH, AMOC, MLD, and Atlantic section of temperature are shown in Fig. 8, indicating
<figure><img src="image_6.png" /><figcaption>FIG. 5. As in Fig. 4, but for AMOC. Solid line contour interval is 1 Sv; color contours are 0.2 Sv.</figcaption></figure>

marked improvements compared to other simulations in every respect. The error on the maximum MLD has been particularly well reduced in the Southern Ocean, confirming that large MLD errors in the cab–therm simulation were mainly a result of the overestimation of the thermal expansion values there. As a result, the temperature distribution is dramatically improved, with a global-mean difference of only 0.03°C with the reference run.

The upper thermocline is too warm, which is probably a result of thermal expansion values being too large in a very thin surface layer, acting to increase locally the stratification and in turn to isolate more surface subtropical waters from the thermocline waters. A similar positive feedback is at play in all simulations with a simplified equation of state, resulting in a more or less marked warm bias in subtropical near-surface waters. Higher-order terms in the simplified EOS would be required to reduce this warm bias acting on waters above 25°C.

4. Conclusions and discussion

In this study, the global-scale impacts of nonlinearities
in the EOS have been assessed using a holistic approach.

The NEMO-based ORCA2 ocean configuration has been
run in forced mode to near equilibrium using a hierarchy
of increasingly complex nonlinear equations of states.
Coefficients of these simplified EOSs were determined
using a least squares method aimed at minimizing the
RMS error on the horizontal density gradients. This
strategy has been followed because the horizontal density
gradient is the only term involving density in the hori-
zontal momentum equations for a Boussinesq hydrostatic
ocean model. It was found that the most important non-
linearity of the equation of state relates to cabbeling,
which is well captured by a quadratic temperature de-
pendence of density, or equivalently, a linear variation of
thermal expansion as a function of temperature. This
nonlinearity greatly reduces the sensitivity of ocean dy-
namics to temperature variations near the freezing point,
thus increasing the relative importance of salinity in polar
regions. The second nonlinearity relates to thermobar-
icity through a linear increase of thermal expansion with
depth, which makes polar ocean dynamics sensitive again
to temperature variations at depth. This has important
effects on the bottom stratification of the ocean and on
the weakly stratified ACC structure and pathway.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Distribution of maximum monthly mean mixed layer depth (maxMLD) for (a) the reference run, (b) the linear EOS, (c) the cabbeling EOS, and (d) the cabbeling-thermobaricity EOS (solid contour interval: 100 m, color scale saturated above 500-m depth).</figcaption></figure>

One secondary goal of this study has been to assess the overall realism that can be obtained for a simulation using a simplified EOS with a form similar to the one proposed by Vallis (2006). Indeed, this EOS form is the simplest possible one that includes representation of both cabbeling and thermobaricity effects in a way that makes it straightforward to distinguish their relative contribution. It was found that ocean properties are very satisfactorily reproduced with such simplified EOS, especially when the thermal

expansion coefficient is accurately reproduced at low temperature.

Two polar regions, although small in spatial extent, appear particularly important in setting the large-scale ocean stratification: the Weddell Sea, where most of Antarctic Bottom Water is formed, and the North Atlantic Bottom Water is formed, and the North Atlantic Deep Water occurs. Our model results must be taken with some caution as the representation of deep water formation in a numerical model is known to be problematic.

<table><caption>TABLE 4. Mean properties of seawater in ORCA2-LIM ocean simulations. Climatological values are based on the PHC3.0 product (Steele et al. 2001).</caption><thead><tr><th></th><th>Θ</th><th>S<sub>A</sub></th><th>a</th><th>b</th><th>N (cycles h<sup>-1</sup>)</th></tr></thead><tbody><tr><td>ref (TEOS-10)</td><td>3.77 ± 4.72</td><td>34.74 ± 0.42</td><td>0.16 ± 0.04</td><td>0.77 ± 0.01</td><td>1.15 ± 1.32</td></tr><tr><td>linear</td><td>1.94 ± 5.73</td><td>34.53 ± 0.46</td><td>0.18</td><td>0.77</td><td>1.20 ± 1.24</td></tr><tr><td>cab</td><td>3.52 ± 5.14</td><td>34.67 ± 0.43</td><td>0.12 ± 0.05</td><td>0.77</td><td>1.15 ± 1.40</td></tr><tr><td>cab-therm</td><td>3.31 ± 5.17</td><td>34.67 ± 0.42</td><td>0.15 ± 0.04</td><td>0.77</td><td>1.21 ± 1.40</td></tr><tr><td>freez</td><td>3.74 ± 4.98</td><td>34.72 ± 0.41</td><td>0.17 ± 0.05</td><td>0.77</td><td>1.21 ± 1.41</td></tr><tr><td>2order</td><td>3.36 ± 5.15</td><td>34.68 ± 0.42</td><td>0.15 ± 0.04</td><td>0.77 ± 0.01</td><td>1.21 ± 1.40</td></tr><tr><td>climatology</td><td>3.90 ± 4.46</td><td>34.89 ± 0.35</td><td>0.16 ± 0.04</td><td>0.77 ± 0.01</td><td>1.09 ± 1.28</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>FIG. 7. As in Fig. 4, but for Conservative Temperature. Solid line contour interval is (a) 4°C for the reference run, and (b)–(d) 1°C for anomaly plots.</figcaption></figure>

However, considering that around 80% of the ocean volume is formed by deep convection in polar regions, it is reasonable to conclude that the large-scale stratification strongly depends on the level of accuracy of the EOS in polar regions.

Based on the findings of this study, we propose the following equation of state as the simplest, yet “realistic,” for seawater:

$$ \rho' = -\frac{C_b}{2}(\Theta - \Theta_o)^2 - T_h Z \Theta + b_o S_A, \quad (17) $$

with

$$ \left\{ \begin{array}{l} C_b = 0.011 \text{ kg m}^{-3} \text{ K}^{-2} \\ T_h = 2.5 \times 10^{-5} \text{ kg m}^{-4} \text{ K}^{-1} \\ b_o = 0.77 \text{ kg m}^{-3} (\text{g kg}^{-1})^{-1} \\ \Theta_o = -4.5^{\circ}\text{C} \end{array} \right. $$

Note that this simplified form of EOS is formally equivalent to the polynomial form in Eq. (6) with $R_{000} = -C_b\Theta_o^2/2$ (whose value has no effect on ocean

dynamics), $R_{100} = b_o$, $R_{010} = \Theta_o C_b$, $R_{020} = -C_b/2$, and $R_{011} = -T_h$.

One difficulty with the polynomial model of simplified EOS proposed in Eq. (6) is that the different terms are not easy to interpret. In Eq. (17), we propose a different form of simplified EOS, with a better design because each parameter has a straightforward interpretation:

1) $C_b$ sets the sensitivity of thermal expansion to temperature;

2) $\Theta_o$ corresponds to the temperature at which surface thermal expansion is zero;

3) $T_h$ gives the sensitivity of thermal expansion to depth; and

4) $b_o$ is the constant haline contraction.

The second parameter is of utmost importance because it defines where salinity will become comparatively important in determining the stratification and the circulation. The cabbeling property depends on the two first parameters taken together, while thermobaricity is set by the third parameter.
<figure><img src="image_9.png" /><figcaption>FIG. 8. Anomaly relative to the reference run of properties in the model run using the freezing EOS. Anomalies are shown for (a) the SSH, (b) the AMOC, (c) the maxMLD, and (d) the Conservative Temperature using the same color scale and contours as in Figs. 4–7, respectively.</figcaption></figure>

The actual value of $\Theta_o$ refers to the fact that, for average oceanic salinity conditions (i.e., $S_A \approx 35 \text{ g kg}^{-1}$), the thermal expansion coefficient would have vanished at around $\Theta = -4.5^\circ\text{C}$ if the seawater was not yet frozen. Note that in reality the temperature at which thermal expansion vanishes follows a nearly linear law in salinity, having a positive value of about $4^\circ\text{C}$ for freshwater. A slightly more complex EOS can be defined to allow such a weak salinity dependence, adding a new parameter $\epsilon$ as follows:

$$ \rho' = -\frac{C_b}{2}(\Theta - \Theta_o - \epsilon S_A)^2 - T_h Z \Theta + b_o S_A, \quad (18) $$

with

$$ \left\{ \begin{array}{l} C_b = 0.011 \text{ kg m}^{-3} \text{ K}^{-2} \\ T_h = 2.5 \times 10^{-5} \text{ kg m}^{-4} \text{ K}^{-1} \\ b_0 = 0.77 \text{ kg m}^{-3} (\text{g kg}^{-1})^{-1} \\ \Theta_o = 4^\circ\text{C} \\ \epsilon = -0.25 \text{ K(g kg}^{-1})^{-1} \end{array} \right. $$

This form of simplified EOS has not been tested with the ORCA2 configuration, but it is very similar to the 2order simplified EOS (see Table 3). Owing to the very small differences obtained when implementing the cab-therm or 2order forms (see Table 4), we do not expect any major global-scale differences between simulations based on Eqs. (17) or (18). However, it is a central feature distinguishing brackish seas such as the Baltic Sea from the large-scale salty ocean, and it may also play a more substantial role in relatively fresh regions such as the tropical eastern Pacific or the Arctic Ocean.

In this study, we have used the Boussinesq and hydrostatic approximations (Vallis 2006). The Boussinesq approximation mainly consists of replacing the equation of mass conservation with a simpler equation of volume conservation (incompressibility), while the hydrostatic approximation neglects vertical accelerations. This set of approximations is often used in ocean modeling because it greatly simplifies the equations of motion while introducing negligible errors with respect to other model uncertainties (Losch et al. 2004). We therefore argue
that our conclusions are essentially insensitive to these assumptions, in the sense that the model sensitivity due to hydrostatic Boussinesq assumptions (introducing ~1% differences on the tracer fields) is well below the sensitivity to the use of a simplified EOS (typically 10%–50% differences).

Yet, an increasing number of ocean models are relaxing the Boussinesq approximation because they then have a more transparent treatment of sea level variations by directly accounting for the steric contribution (e.g., Griffies and Adcroft 2008). For such models, compressible equivalents of these simplified equations of state can be easily derived, using the approximation that specific volume $v$ has small relative variations $v = v_o + v'$, with $v_o \approx 10^{-3} \text{ m}^3\text{kg}^{-1}$ and

$$v'(\Theta, S_A, p) \approx -v_o^2 \rho'(\Theta, S_A, Z), \quad (19)$$

where $\rho'$ is defined by Eq. (17).

This work provides useful insights on the strengths and limitations of existing studies using an idealized equation of state (e.g., Hogg et al. 2013; Abernathey et al. 2011) by presenting basic characteristics of simulated ocean circulations induced with a linear or weakly nonlinear equation of state. In particular, it is found that the “bottom water” overturning cell in simple models with a linear EOS is expected to be very intense, albeit for the wrong reasons related to the heat forcing rather than the salt forcing in the real world. Because nonlinearities of the EOS have such a profound effect on the large-scale stratification, close to the surface but also at depth, it is simply impossible to obtain a realistic thermohaline circulation from a linear EOS. It emphasizes the need for a sufficiently accurate, fundamentally nonlinear EOS in ocean theories and models when one wants to explain and simulate the large-scale stratification.

It is rather astonishing to notice that although nonlinear terms must be accounted for, an EOS as simple as the one proposed in Eq. (17) suffices to reproduce with good accuracy observed large-scale properties of the ocean. Furthermore, the coefficients found here using the optimal least squares method do not differ so much from the simple estimates first proposed by Vallis (2006). In his book, Vallis (2006) simply took the values of the different EOS parameters at salinity $S_A = 35 \text{ g kg}^{-1}$, temperature $\theta = 10^\circ\text{C}$, and surface pressure, equivalent to a first-order Taylor expansion of the real EOS around these reference values. We were able to demonstrate herein that Vallis's (2006) simplified EOS is quite accurate in practice, although our proposed simplified EOS should be preferred in the future.

The availability of a simplified yet realistic EOS will be particularly useful for theoretical studies trying to

gain an understanding of the nature and variability of the meridional overturning circulation. For example, simple one dimensional models such as thermohaline loops (Wunsch 2005) can serve as metaphors for the circulation, potentially useful for the development of theoretical concepts and as a pedagogical tool, owing to their drastically simplified dynamics and yet rich thermodynamics. Despite their formal simplicity, thermohaline loops can feature chaotic, bistable, or oscillating behaviors. A recent study by Pollmann et al. (2015) shows that the use of a nonlinear EOS can substantially modify the behavior of the thermohaline loop, suggesting that the nonlinear EOS might have a profound impact on the nature of the meridional overturning circulation.

Important effects have been observed in ocean simulations by changing the form of the equation of state. Yet, the magnitude of effects was limited by the use of fixed forcing fields at the ocean surface. It would be interesting to perform a similar study with a coupled ocean-atmosphere climate model, assessing what would be the global climate sensitivity to a change in the seawater equation of state. It is very likely that the climate state, including the mean state of the atmosphere, would be profoundly modified by a change of the equation of state. In this sense, we argue that the form of the seawater equation of state is a major determinant of the Earth climate system.

*Acknowledgments.* We thank Trevor McDougall and an anonymous reviewer for their helpful and constructive comments. This work has been funded by the Bolin Centre for Climate Research, research area “Oceans, atmosphere dynamics and climate.” The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC, Linköping.

## REFERENCES

* Abernathey, R., J. Marshall, and D. Ferreira, 2011: The dependence of Southern Ocean meridional overturning on wind stress. *J. Phys. Oceanogr.*, **41**, 2261–2278, doi:10.1175/JPO-D-11-023.1.

* Adcroft, A., R. Hallberg, and M. Harrison, 2008: A finite volume discretization of the pressure gradient force using analytic integration. *Ocean Modell.*, **22**, 106–113, doi:10.1016/j.ocemod.2008.02.001.

* Blanke, B., and P. Delecluse, 1993: Variability of the tropical Atlantic Ocean simulated by a general circulation model with two different mixed-layer physics. *J. Phys. Oceanogr.*, **23**, 1363–1388, doi:10.1175/1520-0485(1993)023<1363:VOTTAO>2.0.CO;2.

* Fichefet, T., and M. A. M. Maqueda, 1997: Sensitivity of a global sea ice model to the treatment of ice thermodynamics and dynamics. *J. Geophys. Res.*, **102**, 12 609–12 646, doi:10.1029/97JC00480.

* Gent, P. R., J. Willebrand, T. J. McDougall, and J. C. McWilliams, 1995: Parameterizing eddy-induced tracer transports in ocean
circulation models. *J. Phys. Oceanogr.*, **25**, 463-474, doi:10.1175/1520-0485(1995)025<0463:PEITTI>2.0.CO;2.

Griffies, S. M., 2004: *Fundamentals of Ocean Climate Models*. Princeton University Press, 518 pp.

—, and A. J. Adcroft, 2008: Formulating the equations of ocean models. *Ocean Modeling in an Eddying Regime, Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 281-317.

—, and Coauthors, 2009: Coordinated Ocean-Ice Reference Experiments (CORES). *Ocean Modell.*, **26**, 1-46, doi:10.1016/j.ocemod.2008.08.007.

Held, I. M., 2005: The gap between simulation and understanding in climate modeling. *Bull. Amer. Meteor. Soc.*, **86**, 1609-1614, doi:10.1175/BAMS-86-11-1609.

Hogg, A. M., H. A. Dijkstra, and J. A. Saenz, 2013: The energetics of a collapsing meridional overturning circulation. *J. Phys. Oceanogr.*, **43**, 1512-1524, doi:10.1175/JPO-D-12-0212.1.

IOC, SCOR, and IAPSO, 2010: The International Thermodynamic Equation of Seawater—2010: Calculation and use of thermodynamic properties. Intergovernmental Oceanographic Commission, Manuals and Guides 56, 220 pp. [Available online at http://www.teos-10.org/pubs/TEOS-10_Manual.pdf.]

Iudicone, D., G. Madec, B. Blanke, and S. Speich, 2008: The role of Southern Ocean surface forcings and mixing in the global conveyor. *J. Phys. Oceanogr.*, **38**, 1377-1400, doi:10.1175/2008JPO3519.1.

Killworth, P. D., 1983: Deep convection in the World Ocean. *Rev. Geophys.*, **21**, 1-26, doi:10.1029/RG021i001p00001.

Klocker, A., and T. J. McDougall, 2010: Influence of the nonlinear equation of state on global estimates of dianeutral advection and diffusion. *J. Phys. Oceanogr.*, **40**, 1690-1709, doi:10.1175/2010JPO4303.1.

Large, W., and S. Yeager, 2004: Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies. NCAR Tech. Note NCAR/TN-460+STR Tech. Rep., 105 pp., doi:10.5065/D6KK98Q6.

—, and —, 2009: The global climatology of an interannually varying air-sea flux data set. *Climate Dyn.*, **33**, 341-364, doi:10.1007/s00382-008-0441-3.

Locarnini, R. A., and Coauthors, 2013: *Temperature*. Vol. 1, *World Ocean Atlas 2013*, NOAA Atlas NESDIS 73, 40 pp.

Losch, M., A. Adcroft, and J.-M. Campin, 2004: How sensitive are coarse general circulation models to fundamental approximations in the equations of motion? *J. Phys. Oceanogr.*, **34**, 306-319, doi:10.1175/1520-0485(2004)034<0306:HSACGC>2.0.CO;2.

Madec, G., and M. Imbard, 1996: A global ocean mesh to overcome the North Pole singularity. *Climate Dyn.*, **12**, 381-388, doi:10.1007/BF00211684.

—, and Coauthors, 2014: Nemo ocean engine, version 3.6. Note du Pole de Modelisation de l'Institut Pierre-Simon Laplace Tech. Rep. 27, 392 pp. [Available online at http://www.nemo-ocean.eu/content/download/164969/671898/file/NEMO_book_V3_6.pdf.]

McDougall, T. J., 1987: Neutral surfaces. *J. Phys. Oceanogr.*, **17**, 1950-1964, doi:10.1175/1520-0485(1987)017<1950:NS>2.0.CO;2.

—, 2003: Potential enthalpy: A conservative oceanic variable for evaluating heat content and heat fluxes. *J. Phys. Oceanogr.*, **33**, 945-963, doi:10.1175/1520-0485(2003)033<0945:PEACOV>2.0.CO;2.

—, and D. R. Jackett, 2005: The material derivative of neutral density. *J. Mar. Res.*, **63**, 159-185, doi:10.1357/0022240053693734.

—, —, F. J. Millero, R. Pawlowicz, and P. M. Barker, 2012: A global algorithm for estimating Absolute Salinity. *Ocean Sci.*, **8**, 1123-1134, doi:10.5194/os-8-1123-2012.

Mignot, J., D. Swingedouw, J. Deshayes, O. Marti, C. Talandier, R. Séférían, M. Lengaigne, and G. Madec, 2013: On the evolution of the oceanic component of the IPSL climate models from CMIP3 to CMIP5: A mean state comparison. *Ocean Modell.*, **72**, 167-184, doi:10.1016/j.ocemod.2013.09.001.

Pollmann, F., F. Roquet, and G. Madec, 2015: Effects of the asymmetry between surface and interior flow on the dynamics of a thermohaline loop. *J. Phys. Oceanogr.*, doi:10.1175/JPO-D-15-0022.1, in press.

Roquet, F., 2013: Dynamical potential energy: A new approach to ocean energetics. *J. Phys. Oceanogr.*, **43**, 457-476, doi:10.1175/JPO-D-12-098.1.

—, G. Madec, T. J. McDougall, and P. M. Barker, 2015: Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard. *Ocean Modell.*, **90**, 29-43, doi:10.1016/j.ocemod.2015.04.002.

Roullet, G., and G. Madec, 2000: Salt conservation, free surface, and varying levels: A new formulation for ocean general circulation models. *J. Geophys. Res.*, **105**, 23 927-23 942, doi:10.1029/2000JC900089.

Steele, M., R. Morley, and W. Ermold, 2001: PHC: A global ocean hydrography with a high-quality Arctic Ocean. *J. Climate*, **14**, 2079-2087, doi:10.1175/1520-0442(2001)014<2079: PAGOHW>2.0.CO;2.

Sun, S., R. Bleck, C. Rooth, J. Dukowicz, E. Chassignet, and P. Killworth, 1999: Inclusion of thermobaricity in isopycnic-coordinate ocean models. *J. Phys. Oceanogr.*, **29**, 2719-2729, doi:10.1175/1520-0485(1999)029<2719:IOTIIC>2.0.CO;2.

Timmermann, R., H. Goosse, G. Madec, T. Fichefet, C. Ethe, and V. Duliere, 2005: On the representation of high latitude processes in the ORCA-LIM global coupled sea ice-ocean model. *Ocean Modell.*, **8**, 175-201, doi:10.1016/j.ocemod.2003.12.009.

Vallis, G. K., 2006: *Atmospheric and Oceanic Fluid Dynamics*. Cambridge University Press, 745 pp.

Wunsch, C., 2005: Thermohaline loops, Stommel box models, and the Sandström theorem. *Tellus*, **57A**, 84-99, doi:10.1111/j.1600-0870.2005.00093.x.

Young, W. R., 2010: Dynamic enthalpy, Conservative Temperature, and the seawater Boussinesq approximation. *J. Phys. Oceanogr.*, **40**, 394-400, doi:10.1175/2009JPO4294.1.