Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard

Ocean Modelling 90 (2015) 29–43

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard F. Roquet a,∗, G. Madec b,c, Trevor J. McDougall d, Paul M. Barker d a

Department of Meteorology (MISU), Stockholm University, S-106 91 Stockholm, Sweden LOCEAN Sorbonne Universits (UPMC, Univ Paris 6)/CNRS/IRD/MNHN, Paris 75252, France c National Oceanography Centre, Southampton, Marine Systems Modelling Group, European Way, Southampton SO14 3ZH, UK d School of Mathematics and Statistics, University of New South Wales, Sydney, New South Wales, Australia b

a r t i c l e

i n f o

Article history: Received 14 October 2014 Revised 7 April 2015 Accepted 8 April 2015 Available online 17 April 2015 Keywords: Physical oceanography Equation of state TEOS-10 Ocean thermodynamics Ocean modelling Seawater density

a b s t r a c t A new set of approximations to the standard TEOS-10 equation of state are presented. These follow a polynomial form, making it computationally efficient for use in numerical ocean models. Two versions are provided, the first being a fit of density for Boussinesq ocean models, and the second fitting specific volume which is more suitable for compressible models. Both versions are given as the sum of a vertical reference profile (6th-order polynomial) and an anomaly (52-term polynomial, cubic in pressure), with relative errors of 0.1% on the thermal expansion coefficients. A 75-term polynomial expression is also presented for computing specific volume, with a better accuracy than the existing TEOS-10 48-term rational approximation, especially regarding the sound speed, and it is suggested that this expression represents a valuable approximation of the TEOS-10 equation of state for hydrographic data analysis. In the last section, practical aspects about the implementation of TEOS-10 in ocean models are discussed.

1. Introduction The equation of state is an empirical thermodynamic relationship linking seawater density to a number of state variables, most typically temperature, salinity and pressure. Because density gradients control the pressure gradient force through the hydrostatic balance, the equation of state provides a fundamental bridge between the distribution of active tracers and the fluid dynamics. Recently, TEOS-10 (International Thermodynamic Equation Of Seawater – 2010, IOC et al. (2010)) has been released as the official replacement for the previous EOS-80 standard. Apart from an improved accuracy of the formulation owing to considerations of a greater number of refined laboratory measurements, this standard introduced two new concepts to be used as the conservative state variables: Absolute Salinity (unit: g/kg, notation: SA ) and Conservative Temperature (unit: o C, notation: ). All thermodynamical variables, including seawater density, are derived analytically from the seawater Gibbs function Feistel (2008), and are given as a function of SA ,  and the gauge pressure p (unit: dbar). Because the exact expression for the TEOS-10 seawater density (or equivalently, specific volume) is complex and rather computing expensive, it is not well suited for practical applications such as ∗

Corresponding author. Tel.: +46 8 16 4338. E-mail address: [email protected] (F. Roquet).

http://dx.doi.org/10.1016/j.ocemod.2015.04.002 1463-5003/© 2015 Elsevier Ltd. All rights reserved.

© 2015 Elsevier Ltd. All rights reserved.

hydrographic data analysis or ocean modeling. This is why an approximation has been proposed in TEOS-10, which consists of a 48-term rational function, with an accuracy better than field measurement accuracy over the whole oceanographic range of (SA ,, p) values. Although this expression is very well suited for most practical applications, it is not so well adapted for ocean modeling applications because it remains difficult to compute analytically its partial derivatives or integrals. In particular, two quantities are key and need to be computed at least as often as the seawater density in an ocean model, namely the thermal expansion coefficient α and the haline contraction coefficient β , which are proportional to the  and SA partial derivatives, respectively,

α(SA , , p) =

 1 ∂ v  1 =− v ∂  SA ,p ρ

∂ρ  , ∂  SA ,p

 1 ∂ v  1 = v ∂ SA ,p ρ

∂ρ  , ∂ SA ,p

β(SA , , p) = −



(1)



(2)

where v is the specific volume, and ρ = 1/v is the density. The thermal expansion and haline contraction coefficients are useful to estimate the pressure gradient force term (through computation of the isobaric gradient of specific volume, p v = (v/x, v/y) at constant p), and the buoyancy frequency N,

∇p v = v [α∇p  − β∇p S] ,

(3)

30

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

 N2 = g

α



∂ ∂ SA −β . ∂z ∂z

(4)

These two coefficients are also important in isoneutral and vertical mixing parameterizations. The equation of state is also important to determine the sound speed (unit: m s−1 ), which is a function of the local compressibility,

 cs =



∂ρ   ∂ P SA ,

−1 ,

(5)

where P is the absolute pressure, P = Po + p (preferably expressed in Pa = 10−4 dbar), with Po = 101, 325 Pa. Finally, one quantity of great interest for the analysis of ocean energetics is enthalpy,

h(SA , , P ) = cpo  + h† ,

(6)

where the first right hand side (r.h.s.) term is the potential enthalpy (McDougall, 2003), and the second is dynamic enthalpy (Young, 2010) P defined as h† = Po v(SA , , P  )dP  . Potential enthalpy is an excellent measure of heat content. Dynamic enthalpy is central to the ocean potential energy budget, as it represents the amount of total potential energy that would be needed to lift the water parcel adiabatically to the level P = Po . In this paper, we present a set of TEOS-10 polynomial approximations which are specifically designed for ocean models. The proposed TEOS-10 approximations are suitable for ocean models using z- or sigma-vertical coordinates (e.g. Griffies, 2004), and could potentially be adapted to isopycnic coordinate models following standard procedures (Adcroft et al., 2008). A fundamental distinction must be made however between Boussinesq and non-Boussinesq models, as the former requires the fast computation of density and its partial derivatives, while the latter rather depends on specific volume. Various TEOS-10 approximations for Boussinesq or compressible (i.e. non-Boussinesq) models are proposed in this paper. The accuracy of these approximations is shown to be comparable to the TEOS-10 rational function approximation, but the polynomial expressions have simpler and more computationally efficient expressions for their derived quantities which make them more adapted for use in ocean models. Also, a high precision form is presented, that we suggest as a replacement of the TEOS-10 rational function approximation for hydrographic data analysis. The paper is organized as follows. Elements related to the implementation of Boussinesq ocean models will be provided in Section 2, then polynomial approximations for Boussinesq models Section 3 and compressible models (Section 4) will be described. In Section 5, practical aspects on the implementation of TEOS-10 in ocean models will be discussed, and conclusions will finally be drawn in Section 6.

2. The Boussinesq approximation in ocean models Currently, a majority of ocean general circulation models (OGCMs) are making both the hydrostatic and Boussinesq approximations (e.g. Vallis, 2006). The Boussinesq approximation is well adapted to the ocean case because seawater is nearly incompressible, and all dynamical effects of specific volume variations can be safely neglected except for the contribution to buoyancy variations. Ocean dynamics is sensitive at leading order to the horizontal gradient of density instead of density itself, through its effect on the horizontal pressure gradient. In the Boussinesq approximation, the only dynamical term where seawater density plays a role is precisely the horizontal density gradient term (z ρ = (ρ /x, ρ /y) at constant depth z), as other terms are proportional to a constant reference density ρ o instead of density itself,

∇z p(z) = ∇z p(z = 0) + g



0 z

∇z ρ dz ,

(7)

where p is the hydrostatic pressure, ρ is the in situ density, and g is the gravitational acceleration assumed constant. By convention, the vertical axis is pointing upward, so that depth is zero at the surface, and negative in the interior of the ocean. In a Boussinesq fluid, the equation of state is a function of the depth z (as it must be a function of the static pressure − ρ o gz), instead of the in situ pressure (Vallis, 2006; Shchepetkin and McWilliams, 2011). This is required to get a consistent description of the energetics of the Boussinesq fluid model (Young, 2010). This can be easily understood by the fact that in a Boussinesq fluid, it is the volume which is conservative, in contrast with the compressible case where it is mass which −1 is conservative. In the following, we simply used z = −p × 1 m dbar for the polynomial fits, equivalent to using a value ρ 0  1020 kg m−3 for the conversion. A simple relationship exists between horizontal gradients of ρ ,  and SA ,

∇z ρ = −a∇z  + b∇z S,

(8)

where a and b are the “thermal expansion” and “haline contraction” coefficients, respectively,

a(SA , , z) = −

∂ρ  ∂ρ  , b(SA , , z) =   . ∂  SA ,z ∂ SA ,z

(9)

Note that we use here definitions that differ slightly from the usual ones, as they are not divided by density. This form is indeed more suitable for the Boussinesq model, as is apparent in the previous Eq. (8). As a consequence, a and b have here units of [kg m−3 ◦ C−1 ] and [kg m−3 (g/kg)−1 ], respectively. Typical values of these coefficients in the world ocean are given in Table 1.

Table 1 Mean properties of seawater in the global ocean, based on the WOA13 climatology (Locarnini et al., 2013). Polar regions are taken poleward of 60o of latitude, while tropical regions are equatorward of 20o . ρ is the in situ density, while r is the density anomaly as defined in Section 3.1. a, b and Nbous are the Boussinesq forms of thermal expansion, haline contraction, and buoyancy frequency, respectively (see Eqs. (9) and (10)). Note the large reduction of thermal expansion, and to a lesser degree, the increase in haline contraction, as we move poleward. The ||z ρ ||-weighted global values correspond roughly to the average thermocline values.

 (o C)

SA (g/kg)

Depth (m)

ρ (kg m−3 )

Global Global (||z ρ ||-weighted) Surface Surface, polar regions Surface, tropical regions

3.9 9.4 18.2 0.0 26.9

34.9 34.9 34.9 33.3 35.2

2050 780 0 0 0

1036.8 1030.1 1024.6 1026.6 1022.8

r (kg m−3 )

a (kg m−3 K−1 )

b (kg m−3 (g/kg)−1 )

Nbous (cycle/h)

Global Global (||z ρ ||-weighted) Surface Surface, polar regions Surface, tropical regions

1027.6 1026.6 1024.6 1026.5 1022.7

0.16 0.18 0.24 0.05 0.32

0.77 0.77 0.76 0.80 0.74

1.08 3.30 5.49 3.42 2.96

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

8

×10

2

2

residual [m3 kg-1 ]

4

0

-2

0

-2

-4

-4

-6

-6

0

-2000

-4000

-6000

-8000

-10000

(b) standard error on v0: 3.3e-12 m3 kg-1

-12

v0 = π . P1(π), P1 5-order polynomial

6

4

-8

×10

8

r0 = ζ . P0(ζ ) , P0 5-order polynomial

6

residual [kg m-3 ]

(a) standard error on r0(Z): 3.3e-06 kg m-3

-6

31

-8

0

2000

4000

6000

8000

10000

pressure [dbar]

depth [m]

Fig. 1. Residual error of the polynomial fit for (a) the density reference profile ro (z) used in polyTEOS10-bsq, and (b) the specific volume reference profile vo (p) used in polyTEOS1055t and polyTEOS10-75t. The coefficients of these polynomials are given in Appendix A for density and Appendix C.1 for specific volume.

The Boussinesq definition of the squared buoyancy frequency is defined as, 2 Nbous =

g





a

ρo

∂ ∂ SA −b . ∂z ∂z

(10)

Contrary to the definition for a compressible fluid, the Boussinesq definition of buoyancy frequency is here inversely proportional to the reference density instead of the density itself. This is another reason why the values a and b are better suited for a Boussinesq fluid than the usual definitions of α and β . The sound speed for the Boussinesq model is simply defined as,

 −ρo g cs,bous =  . ∂ρ  ∂z 

(11)

SA ,

The quantity of interest that plays a similar role to dynamic enthalpy for a Boussinesq fluid is the Boussinesq dynamic enthalpy (Young, 2010), Boussinesq potential energy (Roquet, 2013) or effective potential energy (Nycander, 2011),

EPb = −g

0 z

ρ(SA , , z )dz .

within the range of possible (SA ,, z) oceanic values. Here we present another approximation, referred to as polyTEOS10-bsq, which is particularly well suited for Boussinesq ocean models. This fit is written in the form of a 52-parameter polynomial, cubic in depth, with a similar degree of precision. The main advantage of this formulation is that derived quantities such as the thermal expansion and haline contraction coefficients can be obtained very efficiently without approximation, making it particularly useful for application in an OGCM. It is also possible to get a simple expression of its depth-integral to compute the Boussinesq potential energy. Noting that the dynamic of a Boussinesq fluid is insensitive to the addition of any given function of z, we decompose density into a vertical reference profile and a residual function r(SA ,, z), that we will call density anomaly hereinafter,

ρ(SA , , z) = ro (z) + r(SA , , z)

(13)

with

ro (z)  ρ(35.16504 g/kg, 4 ◦ C, z) − ρ(35.16504 g/kg, 4 ◦ C, 0 m). (14)

(12)

The Boussinesq potential energy has units of J m−3 , and would be strictly equal to the gravitational potential energy if the equation of state was independent of depth. 3. Polynomial approximation of the TEOS-10 equation of state for Boussinesq models

The vertical reference profile ro (z) is obtained by fitting the profile of density anomalies taken at (SA = 35.16504 g/kg,  = 4 ◦ C), values that are close to the global mean values found in the ocean, with a 6-order polynomial. The standard error of the fit is 3.2×10−6 kg m−3 over the range [0–8000 m] (Fig. 1), a value well below the threshold of the smallest density measurement accuracy (about 2×10−3 kg m−3 ). The form of the density anomaly fit is as follows:



3.1. polyTEOS10-bsq: fitting the density anomaly The existing form of the TEOS-10 expression for density that is written as a function of Conservative Temperature is a 48-coefficient rational function with the numerator quadratic in the depth z and the denominator cubic in z. This form is an approximation of the expression derived from the Gibbs function for seawater of Feistel (2008)

s=

SA + rSA SAu

,

τ=

 , u

ζ =−

z zu

(15)

r = Poly1,1 (s, τ ) . ζ 3 + Poly2,2 (s, τ ) . ζ 2 + Poly4,4 (s, τ ) . ζ + Poly6,6 (s, τ ),

(16)

32

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43 Table 2 Error statistics (mean ± std) for the TEOS10-48t rational function, the polyTEOS10-bsq fit of the TEOS10 equation of state and the polyTEOS10-stif stiffened form. The exact expression derived from the Gibbs function is used as the reference. The “funnel” range corresponds to the (SA --z) volume defined in McDougall et al. (2003). The “surface” range corresponds to surface values (z = 0) in the funnel. The “WOA13” range corresponds to the World Ocean Atlas 2013 climatology (annual, 1o −1 resolution) (Locarnini et al., 2013). Finally, the “cube” range corresponds to all values in the (SA = 0/42 g kg ,  = −2/40 ◦ C, z = 0/10, 000 dbar) cube. For WOA13 values, the accuracy is estimated as the ratio between the standard error and the associated global-ocean standard deviation. EOS

Funnel

Surface

TEOS10-48t polyTEOS10-bsq polyTEOS10-stif

0.08 ± 0.46 0.08 ± 0.32 0.08 ± 0.31

0.07 ± 0.30 −0.04 ± 0.33 −0.04 ± 0.33

TEOS10-48t polyTEOS10-bsq polyTEOS10-stif

−0.01 ± 0.07 −0.01 ± 0.09 −0.01 ± 0.09

−0.01 ± 0.06 −0.00 ± 0.05 −0.00 ± 0.05

TEOS10-48t polyTEOS10-bsq polyTEOS10-stif

0.00 ± 0.09 0.00 ± 0.10 0.00 ± 0.10

TEOS10-48t polyTEOS10-bsq polyTEOS10-stif

−0.92 ± 6.84 −1.21 ± 8.49 −1.17 ± 8.21

1 ∂r 1 ∂r . , b= u ∂τ 2sSAu ∂ s

cs,bous =

−ρo g

∂ ro 1 ∂r ∂ z − zu ∂ζ

.

Error on b [10−3 kg m−3 (g/kg)−1 ] 0.01 ± 0.07 0.07 ± 0.02 (0.1%) 0.01 ± 0.15 −0.02 ± 0.03 (0.1%) 0.01 ± 0.15 −0.02 ± 0.03 (0.1%)

−0.13 ± 0.78 −1.60 ± 3.79 −1.65 ± 3.85

−6.40 ± 14.22 −3.07 ± 18.66 −2.90 ± 18.00

(17)

(18)

The Boussinesq potential energy defined in Eq. (12) is the sum,

EPb

b b = EP,o (z) + EP,r (SA , , z)

1.27 ± 3.77 1.83 ± 27.45 −0.31 ± 30.77 −0.00 ± 0.49 1.14 ± 2.61 1.43 ± 3.20

Because both expressions involve reduced-order polynomials, they are faster to compute than the equation of state (only 35 polynomial coefficients each). This represents a significant computational advantage over the rational function form of the equation of state. For example, the corresponding calculations for the 48-term rational function approximation requires 80 multiplications and a division. The sound speed can also be obtained analytically,



Error on ρ [10−3 kg m−3 ] −0.26 ± 0.38 (0.02%) −0.14 ± 0.21 (0.01%) −0.13 ± 0.21 (0.01%)

Cube

Error on a [10−3 kg m−3 K−1 ] −0.05 ± 0.05 (0.1%) −0.02 ± 0.07 (0.1%) −0.02 ± 0.06 (0.1%)

where the Polyn, n are polynomials of nth-order in both s and τ (thus involving (n + 1)(n + 2)/2 coefficients). Reference values are chosen following TEOS-10 recommendations, SAu = 40 × (35.16504 g kg−1 /35), u = 40 ◦ C and zu = 104 m. Potential density with respect to the sea surface (reference pressure of 0 dbar) is simply given by σ o = Poly6, 6 (s, τ ). Coefficients have been determined with a nonlinear least-squares method, using the Optimization Toolbox available in Matlab (see Appendix A). The data points used in the minimization procedure were carefully chosen in the (SA ,, z) space to obtain a good fit in the entire oceanic range as defined in TEOS-10, while getting the best results for typical oceanic values. The parameter rS was first adjusted A simultaneously with the polynomials coefficients, and a round value close to its optimal value has then been selected (32 g/kg in this case) for which the final polynomials coefficients have been determined. Because this parameter is positive, haline contraction remains finite as Absolute Salinity tends toward zero, and density anomaly is well defined even for slightly negative salinity values, which can be advantageous when implemented in a numerical model. Note that the use of half-integer powers of salinity has been retained primarily for practical reasons, as it improved significantly the accuracy of the fit in the low salinity range (i.e., for SA < 10 g/kg). Expressions for a and b can be easily derived analytically,

a=−

WOA13

(19)

between a background profile of Boussinesq potential energy EP, o (z), simply obtained by integrating vertically g.ro (z), and the Boussinesq

Error on cs [cm s−1 ] 2.23 ± 6.71 (0.3%) −1.53 ± 8.68 (0.3%) −1.38 ± 8.27 (0.3%)

0.38 ± 71.02 −56.25 ± 178.77 −35.82 ± 184.52

b , potential energy anomaly EP,r



1 1 . Poly1,1 (s, τ ) . ζ 3 + . Poly2,2 (s, τ ) . ζ 2 4 3  1 . Poly4,4 (s, τ ) . ζ 1 + Poly6,6 (s, τ ) + 2

b EP,r = −gzu ζ .

(20)

which is easily derived analytically, as well as its partial derivatives with respect to SA and . This is another major advantage of using this approximation of TEOS-10 over the standard 48-term rational function formulation in a numerical model. Note that the factorized r.h.s. b is the prodterm in Eq. (20) is simply − gzu ζ = gz, meaning that EP,r uct of a density variable (given by the r.h.s. bracketed expression) with gz, similar to potential energy, except that this density variable is not in situ density in this case. The difference between ρ and this density variable reflects the nonlinear nature of the equation of state. Residual error statistics are provided in Table 2, using the exact TEOS-10 expressions as the reference (see also Fig. 2). More specifically, these exact values are obtained by first calculating the in situ temperature t using the function gsw_t_from_CT of the Gibbs SeaWater Oceanographic Toolbox (McDougall and Barker, 2011) and then calculating the in situ density ρ using gsw_rho_t_exact. It is seen that polyTEOS10-bsq has a similar accuracy to TEOS10-48t. As discussed in Appendix K of IOC et al. (2010), the root mean square (rms) uncertainty in the density data that underlies the TEOS-10 Gibbs function is 4×10−3 kg m−3 . It is clear from Table 2 that all three approximate expressions are comfortably inside this underlying uncertainty limit for seawater properties that are close to oceanic values. The corresponding underlying rms uncertainty in the thermal expansion coefficient is 0.73 ×10−6 K−1 , corresponding to an underlying uncertainly of a of 0.73 ×10−3 kg m−3 K−1 . For seawater properties in the oceanic range, the three approximate expressions match the TEOS-10 Gibbs function a factor of ten better than this (in rms). The error in the fit to the Gibbs function values of b are even smaller in terms of their contribution to density gradients. The rms underlying uncertainty in the speed of sound in seawater is 3.5 cm s−1 , and each of the three approximate expressions have rms errors in their sound speeds of approximately twice this amount, even for seawater values in the oceanic ranges.

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

Conservative Temperature [ C]

40

(a) TEOS10-48t 0 km (rms 3.0e-04)

40

35

35

30

30

25

25

20

20

15

15

10

10

5

5

0

0

33

(b) polyTEOS10-bsq 0 km (rms 3.4e-04)

o

10

Conservative Temperature [ C]

40

20

30

40

(c) TEOS10-48t 2 km (rms 5.6e-04)

10 40

35

35

30

30

25

25

20

20

15

15

10

10

5

5

0

0

20

30

40

(d) polyTEOS10-bsq 2 km (rms 3.8e-04)

o

10

Conservative Temperature [ C]

40

20

30

10

40

(e) TEOS10-48t 4 km (rms 3.7e-04)

40

20

30

40

(f) polyTEOS10-bsq 4 km (rms 3.3e-04)

10 -3 kg m-3 2

35

35

1.5

30

30

1

25

25

0.5

20

20

0

15

15

-0.5

10

10

-1

5

5

o

-1.5

0

0 10 20 30 Absolute Salinity [g/kg]

-2

10 20 30 Absolute Salinity [g/kg]

40

40

Fig. 2. TS diagrams of the difference with the reference TEOS10 density for (left) TEOS10-48t, and (right) polyTEOS10-bsq equations of state. From top to bottom, differences are shown at depths 0, 2 and 4 km deep. The contour interval is 2×10−4 kg m−3 . There is no significant difference in precision between TEOS10-48t and polyTEOS10-bsq within the funnel range.

3.2. polyTEOS10-stif: a stiffened equation of state An alternative approach to implement the equation of state in a Boussinesq ocean model relies on the so-called stiffened density (Dukowicz, 2001; Shchepetkin and McWilliams, 2011),

ρ = r1 (z) × r• (SA , , z)

(21)

where the reference ratio profile r1 (z) = ρ(35.16504g/kg, 4 ◦ C, z)/ ρ(35.16504g/kg, 4 ◦ C, z = 0) is used. Note that different notations are

used compared to Dukowicz (2001), where r1 was denoted r instead, and the stiffened density r• was denoted ρ • . The Boussinesq approximation consists in neglecting variations of the reference ratio profile, by making the approximation r1 (z) = 1. In other words, the implemented equation of state is r• (SA , , z) instead of ρ(SA , , z). In this formulation, an error of magnitude r1 − 1 is made on the thermal expansion and haline contraction coefficients, increasing from 0 at the surface to a few percents at the ocean bottom depth. Shchepetkin and McWilliams (2011) describe a modified-Boussinesq method where variations due to r1 are

34

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

retained, reducing significantly the Boussinesq-approximation errors while maintaining the simplicity of a Boussinesq code. To our knowledge, the stiffened approach of using a multiplicative function of height (as in Eq. (21)) is currently used in the σ -coordinate model ROMS (Shchepetkin and McWilliams, 2003, 2011), and it is implemented as an option in the model POP. It is also implemented in several isopycnic models such as MICOM, HYCOM or HIM, although not necessarily in a Boussinesq context (see e.g. Hallberg, 2005). We now propose the polyTEOS10-stif equation of state, giving polynomial approximations of r1 (z) and r• , suitable for ocean models based on the stiffened approach. The reference ratio profile is approximated with a 6th order polynomial in z. The same polynomial form is used here for the stiffened density r• as for the density anomaly in polyTEOS10-bsq (see Eqs. (15)–(16)). In this case, expressions for a, b and c simply become,

r1 ∂ r • ∂ r• a=− , b= , cs = u ∂τ 2sSu ∂ s r1



−ρo g

r• ∂∂rz1

−

r1 ∂ r • zu ∂ζ

the specific volume anomaly,

v(SA , , p) = vo (p) + δ(SA , , p).

Similarl to the density fit, the vertical reference profile is taken as the 6th-order polynomial fit of the anomaly of specific volume relative to the surface level (p = 0 dbar), for a water parcel with constant properties SA = 35.16504 g/kg and  = 4 ◦ C,

vo (p)  v(35.16504 g/kg, 4 ◦ C, p) − v(35.16504 g/kg, 4 ◦ C, 0 dbar) (24) and the specific volume anomaly is fitted using the following polynomial form,

δ(SA , , p) = Poly1,1 (s, τ ) . π 3 + Poly2,2 (s, τ ) . π 2 + Poly4,4 (s, τ ) . π + Poly6,6 (s, τ ) 

(22)

For the expression of the Boussinesq potential energy, if r1 has been neglected, the depth-integral of g.r• should be considered. The energetics of the modified-Boussinesq model is less trivial, and it is not the purpose of this paper to derive it in detail. Coefficients are provided in the Appendix B. Residual errors for the polyTEOS10-stif fit are almost identical to those of the polyTEOS10bsq. This is expected, because the distributions of density anomaly r and of the stiffened density r• are very similar. Note in particular that, in both cases, the density value is constant with respect to depth when SA = 35.16504 g/kg and  = 4 ◦ C.

4. Compressible case: approximations of specific volume 4.1. polyTEOS10-55t: a polynomial fit of specific volume In the compressible case, the important variable is specific volume, instead of density for the Boussinesq case. Here we apply exactly the same procedure to obtain a polynomial fit of TEOS-10, but applied to v. Hence, specific volume is decomposed into a vertical reference profile (function of gauge pressure, unit: dbar) and the residual, called

(23)

with s =

SA + δSA SAu

,

τ=

p  , π= u pu

(26)

where we recall that Polyn, n (s, τ ) are polynomials of nth-order whose parameters are computed using a nonlinear least-squares fitting method. Reference values are chosen following TEOS-10 recom−1 mendations, SAu = 40 × (35.16504 g kg /35), u = 40 ◦ C and pu = −1 104 dbar. In this case, δS = 24 g kg was determined to be the most A suitable value. A 6-order polynomial fit is first computed for the vertical pro−1 file vo (p), with a standard error of 3.1×10−12 m3 kg over the range [0–8000 dbar]. Then the fit of specific volume anomaly is computed, which follows the same polynomial form as for the density anomaly in the previous section, i.e. with 52 parameters (see Appendix C.2). The resulting equation of state is a polynomial with 55 parameters (52 for δ plus 3 for the higher-order terms of vo ). Residual error statistics are provided in Table 3, using the exact TEOS-10 expressions as the reference. It is seen that polyTEOS10-55t has a similar accuracy to TEOS10-48t. The distributions of specific volume difference with the reference are also shown at various depth in Fig. 3 and similarly differences in sound speed are shown in Fig. 4. Another useful test to assess the fit quality is to derive in situ temperature and compute its accuracy over the range of oceanic values. It can be shown that in situ temperature is related to potential

Table 3 Error statistics (mean ± std) for different approximations of the TEOS10 equation of state (non-Boussinesq form), using the exact expression derived from the Gibbs function as the reference. Definitions of the (SA --p) ranges are same as in Table 2. EOS

Funnel

(25)

Surface

WOA13

Cube

−1

TEOS10-48t polyTEOS10-55t polyTEOS10-75t

−0.08 ± 0.44 −0.08 ± 0.30 −0.02 ± 0.19

−0.07 ± 0.30 0.03 ± 0.32 0.01 ± 0.26

Error on v [10−9 m3 kg ] 0.24 ± 0.35 (0.02%) 0.10 ± 0.20 (0.01%) 0.10 ± 0.17 (0.01%)

−1.21 ± 3.56 3.44 ± 32.32 2.66 ± 17.54

TEOS10-48t polyTEOS10-55t polyTEOS10-75t

−0.01 ± 0.06 −0.02 ± 0.10 −0.00 ± 0.03

−0.01 ± 0.06 −0.00 ± 0.06 −0.00 ± 0.04

Error on α [10−6◦ C−1 ] −0.05 ± 0.05 (0.1%) −0.01 ± 0.06 (0.1%) −0.02 ± 0.04 (0.1%)

−0.01 ± 0.50 1.64 ± 3.61 0.03 ± 1.54

TEOS10-48t polyTEOS10-55t polyTEOS10-75t

0.00 ± 0.08 0.01 ± 0.09 0.01 ± 0.08

−0.00 ± 0.09 0.00 ± 0.14 0.00 ± 0.14

Error on β [10−6 (g/kg)−1 ] 0.07 ± 0.02 (0.1%) −0.01 ± 0.03 (0.1%) −0.01 ± 0.01 (0.1%)

−0.11 ± 0.75 −1.61 ± 3.58 0.58 ± 1.30

TEOS10-48t polyTEOS10-55t polyTEOS10-75t

−0.92 ± 6.84 −1.00 ± 7.86 −0.15 ± 2.37

−6.90 ± 14.50 −3.04 ± 17.09 −2.26 ± 5.37

TEOS10-48t polyTEOS10-55t polyTEOS10-75t

−0.03 ± 0.09 −0.03 ± 0.15 −0.01 ± 0.07

0.00 ± 0.00 0.00 ± 0.00 0.00 ± 0.00

Error on cs [cm s−1 ] 2.23 ± 6.71 (0.3%) −1.25 ± 8.09 (0.3%) −1.34 ± 3.96 (0.2%)

−0.34 ± 75.61 −3.35 ± 206.60 42.49 ± 185.70

Error on t [10−3◦ C] −0.06 ± 0.08 (0.02%) −0.03 ± 0.05 (0.01%) −0.03 ± 0.04 (0.01%)

−0.10 ± 0.77 3.19 ± 8.08 −0.13 ± 2.74

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

(a) TEOS10-48t 0 km (rms 2.9e-10)

Conservative Temperature [ C]

40

40

(b) polyTEOS10-55t 0 km (rms 3.2e-10)

40

35

35

35

30

30

30

25

25

25

20

20

20

15

15

15

10

10

10

5

5

5

0

0

35

(c) polyTEOS10-75t 0 km (rms 2.6e-10)

o

10

Conservative Temperature [ C]

40

20

30

40

(d) TEOS10-48t 2 km (rms 5.4e-10)

0 10

40

20

30

40

10

(e) polyTEOS10-55t 2 km (rms 3.8e-10)

40

35

35

35

30

30

30

25

25

25

20

20

20

15

15

15

10

10

10

5

5

5

20

30

40

(f) polyTEOS10-75t 2 km (rms 1.7e-10)

o

0

0 10

Conservative Temperature [ C]

40

20

30

(g) TEOS10-48t 4 km (rms 3.4e-10)

0 10

40 40

20

30

40

10

(h) polyTEOS10-55t 4 km (rms 3.0e-10)

40

20

30

40

(i) polyTEOS10-75t 4 km (rms 1.6e-10)

×10 -9 m3 kg -1 2

35

35

35

1.5

30

30

30

1

25

25

25

0.5

20

20

20

0

15

15

15

-0.5

10

10

10

5

5

5

0

0

0

o

10 20 30 Absolute Salinity [g/kg]

40

10 20 30 Absolute Salinity [g/kg]

-1 -1.5 -2

40

10 20 30 Absolute Salinity [g/kg]

40

−1

Fig. 3. Difference with the reference TEOS10 specific volume [unit: m3 kg ] for (left) TEOS10-48t, and (middle) polyTEOS10-55t, (right) polyTEOS10-75t. From top to bottom, −1 differences are shown at depths 0, 2000 and 4000 m. The contour interval is 2×10−10 m3 kg . There is no significant difference in precision between TEOS10-48t and polyTEOS1055t, while polyTEOS10-75t is clearly improved.

volume (at fixed values of SA and , see IOC et al. (2010)), so that,

temperature as follows:

t=θ

To + θ + cpo

h† = 104

∂ h†   , ∂  SA ,p

(27)



p 0

vo (p ) dp +



p 0



δ(SA , , p )dp .

(28)

The  partial derivative of dynamic enthalpy is then determined analytically,



−1

where cpo = 3991.86795711963 J kg K−1 is the reference value of specific heat capacity (Graham and McDougall, 2013), and To = 273.15 K is the Celsius zero point. Potential temperature θ is determined using the TEOS-10 algorithm which essentially depends on the definition of specific entropy. The polynomial expression for the dynamic enthalpy h† of seawater is the pressure integral of specific

∂ h†  pu ∂ Poly1,1 π 3 ∂ Poly2,2 π 2 = 104 π . . +  ∂  SA ,p u ∂τ 4 ∂τ 3  ∂ Poly4,4 π ∂ Poly6,6 + + . ∂τ 2 ∂τ

(29)

and the resulting in situ temperature can be compared to the true value obtained using the standard TEOS-10 algorithm. As for other

36

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

Conservative Temperature [ C]

40

(a) TEOS10-48t 0 km (rms 1.6e-01)

40

(b) polyTEOS10-55t 0 km (rms 1.8e-01)

40

35

35

35

30

30

30

25

25

25

20

20

20

15

15

15

10

10

10

5

5

5

(c) polyTEOS10-75t 0 km (rms 5.6e-02)

o

0

0 10

Conservative Temperature [ C]

40

20

30

(d) TEOS10-48t 2 km (rms 2.4e-02)

0 10

40 40

20

30

40

(e) polyTEOS10-55t 2 km (rms 6.5e-02)

10 40

35

35

35

30

30

30

25

25

25

20

20

20

15

15

15

10

10

10

5

5

5

20

30

40

(f) polyTEOS10-75t 2 km (rms 1.3e-02)

o

10

Conservative Temperature [ C]

40

0

0

0 20

30

10

40

(g) TEOS10-48t 4 km (rms 5.8e-02)

40

20

30

10

40

(h) polyTEOS10-55t 4 km (rms 4.3e-02)

40

20

30

40

(i) polyTEOS10-75t 4 km (rms 1.4e-02)

0.2

35

35

35

0.15

30

30

30

0.1

25

25

25

0.05

20

20

20

0

15

15

15

-0.05

10

10

10

5

5

5

0

0

0

o

10 20 30 Absolute Salinity [g/kg]

40

10 20 30 Absolute Salinity [g/kg]

40

-0.1 -0.15 -0.2

10 20 30 Absolute Salinity [g/kg]

40

Fig. 4. Difference with the reference TEOS10 sound velocity [unit: m s−1 ] for (left) TEOS10-48terms, and (middle) polyTEOS10-55t, (right) polyTEOS10-75t. From top to bottom, differences are shown at depths 0, 2000 and 4000 m. The contour interval is 2 cm s−1 .

parameters, the accuracy of in situ temperature with polyTEOS10-55t is similar to the one obtained with the TEOS10-48t rational approximation, although errors grow faster in the cubic range outside the funnel. 4.2. polyTEOS10-75t: an accurate polynomial fit of specific volume The second one, polyTEOS10-75t, is of 5th-order order in the pressure, and thus gets large improvements in the accuracy of the corresponding sound velocity,

δ(SA , , p) = Poly0,0 (s, τ ) . π 5 + Poly1,1 (s, τ ) . π 4 + Poly2,2 (s, τ ) . π 3 + Poly4,4 (s, τ ) . π 2 + Poly5,5 (s, τ ) . π + Poly6,6 (s, τ ). (30)

Note that the interest of separating the vertical reference profile and the specific volume anomaly is strongly reduced in this case, as both expressions are of similar order in pressure (although the highest order in pressure is simply multiplied by a constant and so can easily be absorbed into the vertical reference profile vo (p)). Finally, specific volume is obtained with a polynomial expression involving 75 terms, i.e. 74 for δ (see Appendix C.3) plus 1 for the highest-order term of vo (Appendix C.1). The polyTEOS10-75t fit has a much better accuracy than the two other approximations in almost every respect, except for low salinity values of haline contraction near the surface where TEOS10-48t performs slightly better. Although polyTEOS10-75t needs more parameters than TEOS10-48t, its use remains computationally

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

(a) In-situ density 30W 500

45

1500

2000

2000

40

2500 3000 4000

30

5000 5500

25

-50

0 50 Latitude [deg north]

27

1500 2000

26.5

2500

26.5

2500 26

3500

26

3000 3500

25.5

4000

4500

27.5

1000

27

3000 35

3500

28

500

27.5

1000

1500

(c) Horizontal density anomaly 30W

28

500

1000

Depth [m]

(b) Density anomaly 30W

50

37

25.5

4000

4500

25

4500

25

5000

24.5

5000

24.5

5500

5500 -50

24

0 50 Latitude [deg north]

-50

24

0 50 Latitude [deg north]

Fig. 5. Climatological section along 30o W of (a) the in situ density ρ , (b) the density anomaly r, and (c) the horizontal density anomaly ρ  , computed with the WOA13 climatology. The three density variables provide the same horizontal density gradient everywhere by construction. However, the third definition minimizes the magnitude of vertical gradients, making it the best choice as model density variable for global oceanic simulations.

advantageous because its associated partial derivatives and integrals are easily derived analytically. 5. About the implementation of TEOS-10 in a Boussinesq model In a Boussinesq model, the choice of the ρ o value is important, yet it varies significantly among OGCMs, as it is a matter of personal preference. A usual value is 1035 kg m−3 , the ocean average density. This choice minimizes the error of the pressure gradient at mid-depth (i.e. around 2000 m depth). Another sensible choice would be to take 1030 kg m−3 , which is a typical value of density in the thermocline, minimizing the mean error of the thermal wind. Finally, values around 1025 kg m−3 allow minimizing the pressure gradient error near the surface, where currents are most intense and air-sea fluxes are taking place. Note also that the salinity forcing is generally applied as a freshwater flux, i.e. it is proportional to ρ o , so that a value typical of surface density should minimize the associated error. Vertical variations of density anomalies r are reduced by about one order of magnitude compared to those of the in situ density ρ , however horizontal density gradients are strictly identical for both expressions. This is a much desirable property when horizontal density gradients must be computed between adjacent grid points at two different depths. This is most often the case for terrain-following or σ -coordinate, or even for z-coordinate models with a nonlinear free surface and/or partial cells. The horizontal gradient is indeed obtained as follows:

∇z ρ model = ∇s ρ model −

∂ρ model ∇s z ∂z

(31)

where s is a non-horizontal vertical coordinate. The discretization error due to the non-horizontal alignment of the vertical coordinate will also scale as the vertical gradient of the model density variable, as can be seen from the second r.h.s. term. This is why it is desirable to have a model density variable with the smallest possible vertical gradient. Using density anomaly r instead of in situ density ρ as the model density variable already reduces significantly the mean vertical gradient (compare Figs. 5a and 5b), however it is possible to reduce it even more by removing the mean density profile computed over the considered model domain. For a global configuration of the world ocean, a 6-order polynomial fit of ρ(z) is proposed (see Appendix D.1) based on WOA13 climatological data. It is then preferable to use ρ  = ρ − (ρ(z) − ρ(z = 0)) as the model density variable (see Fig. 5c). As different as the distribution of the three density variables ρ , r, and ρ  might look, it must be kept in mind that their horizontal density gradients are strictly equal everywhere. In other words, their use in an ocean model would produce the same ocean circulation in the limit of vanishing discretization errors.

At the same time, the mean vertical gradient is progressively reduced −1 −1 from −5.3 ± 3.5 kg m−3 km for ρ , to −0.9 ± 3.4 kg m−3 km for −1 r, and finally −0.2 ± 3.0 kg m−3 km for ρ  . For a regional configuration, the horizontal-mean density computed over the considered domain should be used instead of ρ(z) to minimize discretization errors. Alternatively, a generic method of reduction of discretization error on the horizontal density gradients can be implemented to treat cases where the temperature and salinity distribution differ significantly from the actual oceanic climatology, such as in climate models, paleo-simulations or idealized ocean domains. In this case, a mean horizontal density profile can be directly estimated using the temperature and salinity distribution of the considered model simulation, either at time t = 0, or at regular time intervals, and consecutively removed from the equation of state. One might wonder why the reference vertical profile ro (z) used in the polyTEOS10-bsq approximation was defined as the vertical density anomaly relative to the surface taken at constant reference temperature and salinity values, rather than using directly the profile of mean horizontal density as the reference. The reason is simply that the quality of the fit was found to be much better in this way. By construction, the density anomaly variable does not vary at all with depth at reference temperature and salinity values, and it varies rather smoothly with depth in the entire salinity/temperature oceanic range. Fast depth variations have been captured efficiently in the 6-order polynomial fit of the reference vertical density profile, and the remaining depth variations of density can be well reproduced with a cubic fit. This would not be true if the reference vertical density profile had been the horizontal mean density profile. For a non-Boussinesq model, gradients of specific volume along isobars must be computed. A similar error reduction procedure can be followed, based on the mean vertical (i.e. function of pressure) profile of specific volume. A fit for the global-mean profile of specific volume is also provided in Appendix D.2. The dimensional variable σ = ρ model /ρo − 1 should be used in the ocean model code, instead of ρ model , in order to achieve the best possible numerical accuracy. Indeed, values of σ are centered, typically ranging between −10−2 and 10−2 . Similarly, it is advisable to manipb /(ρo gz) − 1 instead of the Boussinesq potential ulate σEb = EP,model P

energy itself, as vertical variations of σEb are comparatively small. It P

is easy to verify that σEb is always well defined for z = 0, with σEb = σ P

P

there. Note that for a linear equation of state, σ = σEb at any depth, P

so that the discrepancy between σ and σEb is solely related to therP

mobaric effects (i.e. the depth-dependence of thermal expansion).

38

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

Another important aspect of the TEOS-10 standard is that, unlike potential temperature (notation: θ ), the Conservative Temperature is not equal to in situ temperature at the surface. It typically differs by a few tenth of ◦ C in the world oceans, especially in the low salinity tropical band, the high salinity subtropical Atlantic regions or the closed Mediterranean and Baltic seas (McDougall, 2003). This implies that model temperatures must be converted into in situ temperature in order to compute heat flux at the surface when using a bulk forcing method (e.g. Brodeau et al., 2010), or when the ocean model is coupled to an atmospheric model. In TEOS-10, the algorithm to compute  from θ uses a straightforward polynomial, while the inverse algorithm to compute θ from  is slightly more computationally expensive (about 3 times), using a modified Newton–Raphson iterative solution procedure (McDougall and Wotherspoon, 2014).

•

•

•

6. Conclusions and discussion The TEOS-10 equation of state represents the new standard, and replaces the previous EOS-80 standard for all oceanographic applications. This includes both hydrographic data analysis and ocean modeling applications. Apart from the necessary compliance to official standard, all ocean models should progressively switch to the TEOS-10 standard because (1) it is more accurate, being based on an updated database of laboratory measurements, and (2) it uses Conservative Temperature and Absolute Salinity (instead of potential temperature and practical salinity), both variables being more suitable for use as model variables (IOC et al., 2010; Graham and McDougall, 2013). Consequently, all comparisons between model and observation data should now be done using  and SA , and every ocean general circulation model should change to the TEOS10 standard. In ocean models, computation of density, thermal compression and haline contraction is generally performed at each time step, on meshes typically made of several millions of grid points. This requires to have at the same time accurate and numerically efficient algorithms. Analysis of hydrographic data is generally less computing intensive, but requires the highest possible degree of precision. Still, it is generally recognized that the reference expressions based on the Gibbs function described in Feistel (2003) (see also IOC et al., 2010) are too computing intensive for most practical applications, especially when the temperature variable is Conservative Temperature. We have been able to find polynomial expressions for the equation of state whose accuracy is greater than that of the laboratory data to which the TEOS-10 Gibbs function was fitted. Four different approximate polynomial expressions have been proposed1 :

1 Matlab and python codes are available as Supplementary Materials to compute density (or specific volume), thermal expansion and haline contraction.

•

The polyTEOS10-bsq equation of state is particularly well adapted to ocean models making the Boussinesq assumption. In this form, model density is obtained from a 55-term polynomial expression, while polynomials for thermal expansion and haline contraction require 35-terms each. One major advantage of the polynomial form is that it is straightforward to derive analytically its depthintegral, and hence the dynamic enthalpy, opening the way for accurate budgets of energy in the model. The polyTEOS10-stif equation of state is the equivalent of the polyTEOS10-bsq for Boussinesq ocean models using the stiffened density approach (Shchepetkin and McWilliams, 2011). It has the same number of parameters as polyTEOS10-bsq, and a very similar accuracy. The polyTEOS10-55t equation of state follows the same structure as for polyTEOS10-bsq, but provides specific volume instead of seawater density. It is designed for non-Boussinesq compressible models, and features the same level of accuracy and numerical efficiency as for the Boussinesq versions, combined with the ability to easily derive dynamic enthalpy. The polyTEOS10-75t equation of state gives specific volume from a 75-term polynomial, featuring a particularly high accuracy suitable for any applications, including hydrographic data analysis. In particular, the accuracy of the derived sound speed is largely better than any other TEOS-10 approximation, including the TEOS-10 48-term rational approximation. At the same time, its polynomial form makes it highly suitable for deriving other thermodynamic properties, such as specific enthalpy.

While the three former approximations are specifically designed for ocean models, achieving a sensible trade-off between accuracy and numerical efficiency, the latter polyTEOS10-75t features a greater accuracy which also makes it suitable for hydrographic data analysis. Supplementary materials Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ocemod.2015.04.002. Acknowledgments We would like to thank the three anonymous reviewers 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.

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

Appendix A. Algorithm to compute rho using polyTEOS10-bsq (for Boussinesq models) A.1. Vertical reference profile of density Zu = 1e4; zz = −Z / Zu R00 = 4.6494977072e+01; R01 = −5.2099962525e+00; R02 = 2.2601900708e−01; R03 = 6.4326772569e−02; R04 = 1.5616995503e−02; R05 = −1.7243708991e−03; r0 = (((((R05 ∗ zz+R04) ∗ zz+R03 ) ∗ zz+R02 ) ∗ zz+R01) ∗ zz+R00) ∗ zz Check value for Z = −1000 m: r0 = 4.59763035 kg m−3 .

A.2. Density anomaly SAu = 40 ∗ 35.16504/35; CTu = 40; Zu = 1e4 deltaS = 32 R000 = 8.0189615746e+02; R100 = 8.6672408165e+02; R200 = −1.7864682637e+03; R300 = 2.0375295546e+03; R400 = −1.2849161071e+03; R500 = 4.3227585684e+02; R600 = −6.0579916612e+01; R010 = 2.6010145068e+01; R110 = −6.5281885265e+01; R210 = 8.1770425108e+01; R310 = −5.6888046321e+01; R410 = 1.7681814114e+01; R510 = −1.9193502195e+00; R020 = −3.7074170417e+01; R120 = 6.1548258127e+01; R220 = −6.0362551501e+01; R320 = 2.9130021253e+01; R420 = −5.4723692739e+00; R030 = 2.1661789529e+01; R130 = −3.3449108469e+01; R230 = 1.9717078466e+01; R330 = −3.1742946532e+00; R040 = −8.3627885467e+00; R140 = 1.1311538584e+01; R240 = −5.3563304045e+00; R050 = 5.4048723791e−01; R150 = 4.8169980163e−01; R060 = −1.9083568888e−01; R001 = 1.9681925209e+01; R101 = −4.2549998214e+01; R201 = 5.0774768218e+01; R301 = −3.0938076334e+01; R401 = 6.6051753097e+00; R011 = −1.3336301113e+01; R111 = −4.4870114575e+00; R211 = 5.0042598061e+00; R311 = −6.5399043664e−01; R021 = 6.7080479603e+00; R121 = 3.5063081279e+00; R221 = −1.8795372996e+00; R031 = −2.4649669534e+00; R131 = −5.5077101279e−01; R041 = 5.5927935970e−01; R002 = 2.0660924175e+00; R102 = −4.9527603989e+00; R202 = 2.5019633244e+00; R012 = 2.0564311499e+00; R112 = −2.1311365518e−01; R022 = −1.2419983026e+00; R003 = −2.3342758797e−02; R103 = −1.8507636718e−02; R013 = 3.7969820455e−01; ! ss = sqrt ( (SA+deltaS)/SAu ); tt = CT / CTu; zz = − Z / Zu; rz3 = R013 ∗ tt + R103 ∗ ss + R003 rz2 = (R022 ∗ tt+R112 ∗ ss+R012) ∗ tt+(R202 ∗ ss+R102) ∗ ss+R002 rz1 = (((R041 ∗ tt+R131 ∗ ss+R031) ∗ tt & & + (R221 ∗ ss+R121) ∗ ss+R021) ∗ tt & & + ((R311 ∗ ss+R211) ∗ ss+R111) ∗ ss+R011) ∗ tt & & + (((R401 ∗ ss+R301) ∗ ss+R201) ∗ ss+R101) ∗ ss+R001 rz0 = (((((R060 ∗ tt+R150 ∗ ss+R050) ∗ tt & & + (R240 ∗ ss+R140) ∗ ss+R040) ∗ tt & & + ((R330 ∗ ss+R230) ∗ ss+R130) ∗ ss+R030) ∗ tt & & + (((R420 ∗ ss+R320) ∗ ss+R220) ∗ ss+R120) ∗ ss+R020) ∗ tt & & + ((((R510 ∗ ss+R410) ∗ ss+R310) ∗ ss+R210) ∗ ss+R110) ∗ ss+R010) ∗ tt & & +(((((R600 ∗ ss+R500) ∗ ss+R400) ∗ ss+R300) ∗ ss+R200) ∗ ss+R100) ∗ ss+R000 r = ( ( rz3 ∗ zz + rz2 ) ∗ zz + rz1 ) ∗ zz + rz0 Check values for SA = 30 g/kg, CT = 10 ◦ C, Z = −1000 m: r = 1022.85377 kg m−3 , ρ = 1027.45140 kg m−3 a = 0.179646281 kg m−3 K−1 , b = 0.765555368 kg m−3 (g/kg)−1 .

Appendix B. Algorithm to compute rho using polyTEOS10-stif (for Boussinesq models) B.1. Vertical reference profile of density ratiso Zu = 1e4; zz = − Z / Zu R10 = 4.5238001132e−02; R11 = −5.0691457704e−03; R12 = 2.1990865986e−04; R13 = 6.2587720090e−05; R14 = 1.5194795322e−05; R15 = −1.6777531159e−06; r1 = (((((R15 ∗ zz+R14) ∗ zz+R13 ) ∗ zz+R12 ) ∗ zz+R11) ∗ zz+R10) ∗ zz + 1.0 Check value for Z = − 1000 m : r1 = 1.0044733.

39

40

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

B.2. Stiffened density Same as for the density anomaly variable (see Appendix A), but with the following coefficients: R000 = 8.0185969865e+02; R100 = 8.6694400009e+02; R200 = −1.7869886804e+03; R300 = 2.0381548497e+03; R400 = −1.2853207957e+03; R500 = 4.3240996601e+02; R600 = −6.0597694893e+01; R010 = 2.6018938602e+01; R110 = −6.5349778881e+01; R210 = 8.1938301114e+01; R310 = −5.7075044019e+01; R410 = 1.7778972887e+01; R510 = −1.9385277061e+00; R020 = −3.7047586671e+01; R120 = 6.1469677240e+01; R220 = −6.0273564965e+01; R320 = 2.9086148686e+01; R420 = −5.4641152062e+00; R030 = 2.1645370911e+01; R130 = −3.3415216177e+01; R230 = 1.9694119502e+01; R330 = −3.1710487657e+00; R040 = −8.3587252144e+00; R140 = 1.1301873449e+01; R240 = −5.3494910385e+00; R050 = 5.4258439083e−01; R150 = 4.7964142211e−01; R060 = −1.9098973897e−01; R001 = 2.1989266503e+01; R101 = −4.2043785572e+01; R201 = 4.8565183139e+01; R301 = −3.0473875201e+01; R401 = 6.5025798027e+00; R011 = −1.3731591552e+01; R111 = −4.3667272405e+00; R211 = 5.2899251242e+00; R311 = −7.1323430544e−01; R021 = 7.4843328725e+00; R121 = 3.1443007263e+00; R221 = −1.8141790613e+00; R031 = −2.6010169530e+00; R131 = −4.9866824866e−01; R041 = 5.5882379134e−01; R002 = 1.1144133231e+00; R102 = −4.5413503279e+00; R202 = 2.7242118606e+00; R012 = 2.8508453933e+00; R112 = −4.4471541371e−01; R022 = −1.5059281279e+00; R003 = 1.9817129968e−01; R103 = −1.7905418130e−01; R013 = 2.5254195922e−01; Check values for SA = 30 g/kg, CT = 10 ◦ C, Z = − 1000 m: r_dot = 1022.87574 kg m−3 , ρ = r1 x r_dot = 1027.45140 kg m−3 a=0.179649406 kg m−3 K−1 , b = 0.765554495 kg m−3 (g/kg)−1 . Appendix C. Algorithm to compute v using polyTEOS10 (for non-Boussinesq models) C.1. Vertical reference profile of specific volume Pu = 1e4; pp = p / Pu V00 = −4.4015007269e−05; V01 = 6.9232335784e−06; V02 = −7.5004675975e−07; V03 = 1.7009109288e−08; V04 = −1.6884162004e−08; V05 = 1.9613503930e−09; v0 = (((((V05 ∗ pp+V04) ∗ pp+V03 ) ∗ pp+V02 ) ∗ pp+V01) ∗ pp+V00) ∗ pp Check value for p = 1000 dbar : v0 = −4.333016903e−06 m3 kg−1 . C.2. Specific volume anomaly: polyTEOS10-55t Note that this equation of state needs 55 terms if both the vertical profile and the specific anomaly are added up together. SAu = 40 ∗ 35.16504/35; CTu = 40; Pu = 1e4 deltaS = 24 V000 = 1.0772899069e−03; V100 = −3.1263658781e−04; V200 = 6.7615860683e−04; V300 = −8.6127884515e−04; V400 = 5.9010812596e−04; V500 = −2.1503943538e−04; V600 = 3.2678954455e−05; V010 = −1.4949652640e−05; V110 = 3.1866349188e−05; V210 = −3.8070687610e−05; V310 = 2.9818473563e−05; V410 = −1.0011321965e−05; V510 = 1.0751931163e−06; V020 = 2.7546851539e−05; V120 = −3.6597334199e−05; V220 = 3.4489154625e−05; V320 = −1.7663254122e−05; V420 = 3.5965131935e−06; V030 = −1.6506828994e−05; V130 = 2.4412359055e−05; V230 = −1.4606740723e−05; V330 = 2.3293406656e−06; V040 = 6.7896174634e−06; V140 = −8.7951832993e−06; V240 = 4.4249040774e−06; V050 = −7.2535743349e−07; V150 = −3.4680559205e−07; V060 = 1.9041365570e−07; V001 = −1.6889436589e−05; V101 = 2.1106556158e−05; V201 = −2.1322804368e−05; V301 = 1.7347655458e−05; V401 = −4.3209400767e−06; V011 = 1.5355844621e−05; V111 = 2.0914122241e−06; V211 = −5.7751479725e−06; V311 = 1.0767234341e−06; V021 = −9.6659393016e−06; V121 = −7.0686982208e−07; V221 = 1.4488066593e−06; V031 = 3.1134283336e−06; V131 = 7.9562529879e−08; V041 = −5.6590253863e−07; V002 = 1.0500241168e−06; V102 = 1.9600661704e−06; V202 = −2.1666693382e−06; V012 = −3.8541359685e−06; V112 = 1.0157632247e−06; V022 = 1.7178343158e−06; V003 = −4.1503454190e−07; V103 = 3.5627020989e−07; V013 = −1.1293871415e−07; ! ss = sqrt ( (SA+deltaS)/SAu ); tt = CT/CTu; pp = p/Pu; vp3 = V013 ∗ tt+V103 ∗ ss+V003 vp2 = (V022 ∗ tt+V112 ∗ ss+V012) ∗ tt+(V202 ∗ ss+V102) ∗ ss+V002 vp1 = (((V041 ∗ tt+V131 ∗ ss+V031) ∗ tt & & + (V221 ∗ ss+V121) ∗ ss+V021) ∗ tt &

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

& + ((V311 ∗ ss+V211) ∗ ss+V111) ∗ ss+V011) ∗ tt & & + (((V401 ∗ ss+V301) ∗ ss+V201) ∗ ss+V101) ∗ ss+V001 vp0 = (((((V060 ∗ tt+V150 ∗ ss+V050) ∗ tt & & + (V240 ∗ ss+V140) ∗ ss+V040) ∗ tt & & + ((V330 ∗ ss+V230) ∗ ss+V130) ∗ ss+V030) ∗ tt & & + (((V420 ∗ ss+V320) ∗ ss+V220) ∗ ss+V120) ∗ ss+V020) ∗ tt & & + ((((V510 ∗ ss+V410) ∗ ss+V310) ∗ ss+V210) ∗ ss+V110) ∗ ss+V010) ∗ tt & & +(((((V600 ∗ ss+V500) ∗ ss+V400) ∗ ss+V300) ∗ ss+V200) ∗ ss+V100) ∗ ss+V000 delta = ( ( vp3 ∗ pp + vp2 ) ∗ pp + vp1 ) ∗ pp + vp0 Check values for SA = 30 g/kg, CT = 10 °C, p =1000 dbar : −1 −1 delta = 9.776150635e−04 m3 kg , v = 9.732820466e−04 m3 kg −1 −1 α = 1.748553121e−04 K , β = 7.450974025e−04 (g/kg) C.3. Specific volume anomaly: polyTEOS10-75t Note that this equation of state needs 75 terms if both the vertical profile and the specific anomaly are added up together. SAu = 40 ∗ 35.16504/35; CTu = 40; Pu = 1e4 deltaS = 24 V000 = 1.0769995862e−03; V100 = −3.1038981976e−04; V200 = 6.6928067038e−04; V300 = −8.5047933937e−04; V400 = 5.8086069943e−04; V500 = −2.1092370507e−04; V600 = 3.1932457305e−05; V010 = −1.5649734675e−05; V110 = 3.5009599764e−05; V210 = −4.3592678561e−05; V310 = 3.4532461828e−05; V410 = −1.1959409788e−05; V510 = 1.3864594581e−06; V020 = 2.7762106484e−05; V120 = −3.7435842344e−05; V220 = 3.5907822760e−05; V320 = −1.8698584187e−05; V420 = 3.8595339244e−06; V030 = −1.6521159259e−05; V130 = 2.4141479483e−05; V230 = −1.4353633048e−05; V330 = 2.2863324556e−06; V040 = 6.9111322702e−06; V140 = −8.7595873154e−06; V240 = 4.3703680598e−06; V050 = −8.0539615540e−07; V150 = −3.3052758900e−07; V060 = 2.0543094268e−07; V001 = −1.6784136540e−05; V101 = 2.4262468747e−05; V201 = −3.4792460974e−05; V301 = 3.7470777305e−05; V401 = −1.7322218612e−05; V501 = 3.0927427253e−06; V011 = 1.8505765429e−05; V111 = −9.5677088156e−06; V211 = 1.1100834765e−05; V311 = −9.8447117844e−06; V411 = 2.5909225260e−06; V021 = −1.1716606853e−05; V121 = −2.3678308361e−07; V221 = 2.9283346295e−06; V321 = −4.8826139200e−07; V031 = 7.9279656173e−06; V131 = −3.4558773655e−06; V231 = 3.1655306078e−07; V041 = −3.4102187482e−06; V141 = 1.2956717783e−06; V051 = 5.0736766814e−07; V002 = 3.0623833435e−06; V102 = −5.8484432984e−07; V202 = −4.8122251597e−06; V302 = 4.9263106998e−06; V402 = −1.7811974727e−06; V012 = −1.1736386731e−06; V112 = −5.5699154557e−06; V212 = 5.4620748834e−06; V312 = −1.3544185627e−06; V022 = 2.1305028740e−06; V122 = 3.9137387080e−07; V222 = −6.5731104067e−07; V032 = −4.6132540037e−07; V132 = 7.7618888092e−09; V042 = −6.3352916514e−08; V003 = −3.8088938393e−07; V103 = 3.6310188515e−07; V203 = 1.6746303780e−08; V013 = −3.6527006553e−07; V113 = −2.7295696237e−07; V023 = 2.8695905159e−07; V004 = 8.8302421514e−08; V104 = −1.1147125423e−07; V014 = 3.1454099902e−07; V005 = 4.2369007180e−09; ! ss = sqrt ( (SA+deltaS)/SAu ); tt = CT/CTu; pp = P/Pu; vp5 = V005 vp4 = V014 ∗ tt + V104 ∗ ss + V004 vp3 = ( V023 ∗ tt + V113 ∗ ss + V013 ) ∗ tt & & + ( V203 ∗ ss + V103 ) ∗ ss + V003 vp2 = ( ( ( V042 ∗ tt + V132 ∗ ss + V032 ) ∗ tt & & + ( V222 ∗ ss + V122 ) ∗ ss + V022 ) ∗ tt & & + ( ( V312 ∗ ss + V212 ) ∗ ss + V112 ) ∗ ss + V012 ) ∗ tt & & + ( ( ( V402 ∗ ss + V302 ) ∗ ss + V202 ) ∗ ss + V102 ) ∗ ss + V002 vp1 = ( ( ( ( V051 ∗ tt + V141 ∗ ss + V041 ) ∗ tt & & + ( V231 ∗ ss + V131 ) ∗ ss + V031 ) ∗ tt & & + ( ( V321 ∗ ss + V221 ) ∗ ss + V121 ) ∗ ss + V021 ) ∗ tt & & + ( ( ( V411 ∗ ss + V311 ) ∗ ss + V211 ) ∗ ss + V111 ) ∗ ss + V011 ) ∗ tt & & + ( ( ( ( V501 ∗ ss + V401 ) ∗ ss + V301 ) ∗ ss + V201 ) ∗ ss + V101 ) ∗ ss + V001 vp0 = ( ( ( ( ( V060 ∗ tt + V150 ∗ ss + V050 ) ∗ tt & & + ( V240 ∗ ss + V140 ) ∗ ss + V040 ) ∗ tt & & + ( ( V330 ∗ ss + V230 ) ∗ ss + V130 ) ∗ ss + V030 ) ∗ tt & & + ( ( ( V420 ∗ ss + V320 ) ∗ ss + V220 ) ∗ ss + V120 ) ∗ ss + V020 ) ∗ tt & & + ( ( ( ( V510 ∗ ss + V410 ) ∗ ss + V310 ) ∗ ss + V210 ) ∗ ss + V110 ) ∗ ss + V010 ) ∗ tt &

41

42

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43

0

0

(b)

1000

1000

2000

2000 pressure [dbar]

depth [m]

(a)

3000

4000

5000

6000 -10

3000

4000

5000

3σ interval 6-order fit r o (z)

-5

0

6000 -10

5

density [kg/m3]

3σ interval 6-order fit vo(z)

-5

0 specific volume [m3/kg]

5 ×10 -6

Fig. D6. Range of (a) in situ density and (b) specific volume as a function of pressure (as determined by the ± three standard deviation interval), computed using the WOA13 climatological data set. Deviation profiles from the global horizontal mean are superimposed for the 6-order polynomial fit as well as for the corresponding reference profile (ro and vo respectively).

& + ((((( V600 ∗ ss + V500 ) ∗ ss + V400 ) ∗ ss + V300 ) ∗ ss + V200 ) ∗ ss + V100 ) ∗ ss + V000 delta = ( ( ( ( vp5 ∗ pp + vp4 ) ∗ pp + vp3 ) ∗ pp + vp2 ) ∗ pp + vp1 ) ∗ pp + vp0 Check values for SA = 30 g/kg, CT = 10 ◦ C, p = 1000 dbar : −1 −1 delta = 9.776149797e−04 m3 kg , v = 9.732819628e−04 m3 kg −1 −1 α = 1.748439401e−04 K , β = 7.451213159e−04 (g/kg) Appendix D. Global-mean vertical profiles based on the WOA13 climatology D.1. Mean in situ density A 6-order polynomial fit of the global-mean in situ density is proposed here (Fig. D.6a), valid between the surface and 10,000 m depth. Below 5,500 m, the values of in situ density that are used correspond to those of a water mass with mean bottom properties (i.e., SA = 34.9 g/kg,  = 0.94 ◦ C). Coefficients of the polynomial fit are, Rm0 = 1025.33940359618; Rm1 = 95.684723022965; Rm2 = −361.123729646476; Rm3 = 1224.7453087754; Rm4 = −2105.34088973636; Rm5 = 1749.80669778349; Rm6 = −559.162341562398 zz = − z ∗ 1e-4 r\_mean = (((((Rm6 ∗ zz+Rm5) ∗ zz+Rm4 ) ∗ zz+Rm3 ) ∗ zz+Rm2) ∗ zz+Rm1) ∗ zz + Rm0 Check value for Z = −1000 m : r_mean = 1032.32778 kg m−3 . D.2. Mean specific volume A 6-order polynomial fit of the global-mean specific volume is proposed here (Fig. D.6b), valid between the surface and 10,000 m depth, Vm0 = 9.7528903794e−04; Vm1 = −9.0862908805e−05; Vm2 = 3.4629096173e−04; Vm3 = −1.1681685446e−03; Vm4 = 2.0069833711e−03; Vm5 = −1.6677920142e−03; Vm6 = 5.3289991116e−04 pp = p ∗ 1e−4 v\_mean = (((((Vm6 ∗ pp+Vm5) ∗ pp+Vm4 ) ∗ pp+Vm3 ) ∗ pp+Vm2) ∗ pp+Vm1) ∗ pp + Vm0 Check value for p=1000 dbar : v_mean = 9.686845698e−04 m3 kg

References Adcroft, A., Hallberg, R., Harrison, M., 2008. A finite volume discretization of the pressure gradient force using analytic integration. Ocean Model. 22 (3–4), 106–113. Brodeau, L., Barnier, B., Treguier, A.-M., Penduff, T., Gulev, S., 2010. An ERA40-based atmospheric forcing for global ocean circulation models. Ocean Model. 31 (3–4), 88–104.

−1

.

Dukowicz, J.K., 2001. Reduction of density and pressure gradient errors in ocean simulations. J. Phys. Oceanogr. 31 (7), 1915–1921. Feistel, R., 2008. A Gibbs function for seawater thermodynamics for −6 to 80°C and salinity up to 120 g kg−1 . Deep Sea Res. Part I: Oceanogr. Res. Pap. 55 (12), 1639– 1671. Feistel, R., 2003. A new extended Gibbs thermodynamic potential of seawater. Prog. Oceanogr. 58 (1), 43–114.

F. Roquet et al. / Ocean Modelling 90 (2015) 29–43 Graham, F.S., McDougall, T.J., 2013. Quantifying the nonconservative production of conservative temperature, potential temperature, and entropy. J. Phys. Oceanogr. 43 (5), 838–862. Griffies, S., 2004. Fundamentals of Ocean Climate Models, p. 518, Princeton, NJ, Princeton University Press. Hallberg, R., 2005. A thermobaric instability of lagrangian vertical coordinate ocean models. Ocean Model. 8 (3), 279–300. IOC, SCOR, IAPSO, 2010. The international thermodynamic equation of seawater - 2010: Calculation and use of thermodynamic properties. Intergovernmental Oceanographic Commission, Manuals and Guides No. 56, UNESCO (English), 196 pp. Locarnini, R. A., Mishonov, A. V., Antonov, J. I., Boyer, T. P., Garcia, H. E., Baranova, O. K., Zweng, M. M., Paver, C. R., Reagan, J. R., Johnson, D. R., Hamilton, M., Seidov, D., 2013. World Ocean Atlas, In: S. Levitus, A. Mishonov, Temperature Technical Report, vol. 1, NOAA Atlas NESDIS 73, p. 40, McDougall, T.J., 2003. Potential enthalpy: a conservative oceanic variable for evaluating heat content and heat fluxes. J. Phys. Oceanogr. 33, 945–963. McDougall, T.J., Barker, P.M., 2011. Getting Started with TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox. In: Technical Report, SCOR/IAPSO WG127, p. 28.

43

McDougall, T.J., Jackett, D.R., Wright, D.G., Feistel, R., 2003. Accurate and computationally efficient algorithms for potential temperature and density of seawater. J. Atmos. Oceanic Technol. 20 (5), 730–741. McDougall, T.J., Wotherspoon, S.J., 2014. A simple modification of Newton’s method to achieve convergence of order 1 + 2. Appl. Math. Lett. 29, 20– 25. Nycander, J., 2011. Energy conversion, mixing energy and neutral surfaces with a nonlinear equation of state. J. Phys. Oceanogr. 41, 28–41. Roquet, F., 2013. Dynamical potential energy: a new approach to ocean energetics. J. Phys. Oceanogr. 43 (3), 457–476. Shchepetkin, A.F., McWilliams, J.C., 2003. A method for computing horizontal pressuregradient force in an oceanic model with a nonaligned vertical coordinate. J. Geophys. Res.: Oceans 108 (C3), 3090. Shchepetkin, A.F., McWilliams, J.C., 2011. Accurate Boussinesq oceanic modeling with a practical, “stiffened”, equation of state. Ocean Model. 38 (1–2), 41–70. Vallis, G. K., 2006. Atmospheric and Oceanic Fluid Dynamics, Cambridge University Press, Cambridge. Young, W.R., 2010. Dynamic enthalpy, conservative temperature, and the seawater Boussinesq approximation. J. Phys. Oceanogr. 40, 394–400.