Numerical implementation and oceanographic application of the thermodynamic potentials of liquid water, water vapour, ice, seawater and humid air – Part 1: Background and equations

A new seawater standard referred to as the International Thermodynamic Equation of Seawater 2010 (TEOS10) was adopted in June 2009 by UNESCO/IOC on its 25th General Assembly in Paris, as recommended by the SCOR/IAPSO Working Group 127 (WG127) on Thermodynamics and Equation of State of Seawater. To support the adoption process, WG127 has developed a comprehensive source code library for the thermodynamic properties of liquid water, water vapour, ice, seawater and humid air, referred to as the Sea-Ice-Air (SIA) library. Here we present the background information and equations required for the determination of the properties of single phases and components as well as of phase transitions and composite systems as implemented in the library. All results are based on rigorous mathematical methods applied to the Primary Standards of the constituents, formulated as empirical thermodynamic potential functions and, except for humid air, endorsed as Releases of the International Association for the Properties of Water and Steam (IAPWS). Details of the implementation in the TEOS-10 SIA library are given in a companion paper. Correspondence to: R. Feistel (rainer.feistel@io-warnemuende.de)


Introduction
The recent availability of highly accurate mathematical formulations of thermodynamic potentials for fluid water (Wagner and Pruß, 2002;IAPWS, 2009a), air (Lemmon et al., 2000;Feistel et al., 2010a), ice (Feistel and Wagner, 2006;IAPWS, 2009b), and of seawater (Feistel, 2003(Feistel, , 2008;;Millero et al., 2008;IAPWS, 2008a;IOC et al., 2010) permit the description of thermodynamic properties of the ocean and the atmosphere in an unprecedented comprehensive and consistent manner (Feistel et al., 2008;IOC et al., 2010).In June 2009 these formulations, jointly referred to as the International Thermodynamic Equation of Seawater 2010 (TEOS-10), were adopted by UNESCO/IOC as the successor of the International Equation of State of Seawater 1980 (EOS-80).To assist potential users in implementing the new equations of state, the SCOR 1 /IAPSO 2 Working Group 127 (WG127) on Thermodynamics and Equation of State of Seawater has developed a source code library in Fortran and Visual Basic (also for use with Excel) which provides an extended set of functions for the computation of numerous properties of the geophysical fluids, their composites and phase transitions.Equivalent versions of the SIA library implemented in MatLab and C/C++ are planned.634 R. Feistel et al.: Oceanographic application and numerical implementation of TEOS-10: Part 1 Additional library versions are available for specific application purposes such as the Gibbs SeaWater (GSW) Library for oceanographic models (IOC et al., 2010), consistent with the one described in this paper.This SIA library replaces (i.e., updates and extends) the code versions published previously (Feistel, 2005;Feistel et al., 2005).
In this paper we provide formulas, derivations and explanations for the quantities implemented in the library and their thermodynamic relations to the potential functions.The latter are basic relations which remain valid independent of any details of the fits used to approximate the potentials or the numerics used to evaluate the functional relations; they will not change when new approximations to the thermodynamic potentials are determined in the future.
The oceanographic applications of the quantities discussed here are explained in detail in IOC et al. (2010).Formulas for seawater properties at the ocean-atmosphere interface such as the vapour pressure and the latent heat in equilibrium with humid air consistent with the library functions described here are developed in a separate article (Feistel et al., 2010a).Additional details on the library structure and the numerical implementation are available from a companion paper (Wright et al., 2010a).
The code is organised in vertical columns representing the four constituents: fluid water, ice, seawater and air, and in horizontal levels with increasing complexity at higher levels.The code is hierarchically organised; higher level routines make use of routines from lower levels but code at each level neither needs nor has access to procedures from higher levels.Level 0 provides fundamental constants, mathematical methods and formulae for conversions between Absolute Salinity and Practical Salinity given the measurement location plus some general relations.Level 1 contains the Primary Standards, i.e. the thermodynamic potentials, their basic constants and coefficients, and further necessary empirical or theoretical functions.Levels 2-4 contain properties derived directly or indirectly from level 1 by mathematical and numerical manipulations without additional empirical formulas or constants.Level 2 contains explicit relations.Levels 3 and 4 require iteration procedures to numerically solve implicit equations.Levels 1-3 describe only singlephase properties, while level 4 procedures calculate phase equilibria and properties of composite, multi-phase systems.Level 5 is an additional application layer within which input and output units may be more convenient for the user than the basic SI units that are used rigorously at the lower levels.Level 5 also contains additional empirical coefficients in speed-optimized code, derived as "approximate" correlation functions representing reduced data sets computed from the lower, "exact" levels.
The following sections explain the thermodynamic relations and auxiliary equations implemented in the library level by level, with the exception of level 0, which provides general information required by the library, and level 5, which is only numerically different from the lower ones.Within each level, the sections consider the four columns or their combinations.In the appendix, the formulas are described which are implemented in the library to iteratively solve implicit equations.When appropriate, function names defined in the library are mentioned in the text of this paper along with the related equations; they are emphasized by means of a distinct font type.The complete list of modules and functions available from the library is described in the companion paper, Part 2 (Wright et al., 2010a).
The core levels 1-4 of the SIA library take as input parameters absolute pressure in Pa, absolute ITS-90 temperature in K, and Absolute Salinity in kg/kg.Absolute Salinity of Standard Seawater is most accurately estimated by Reference Salinity (Millero et al., 2008) which is proportional to Practical Salinity in the valid range of the 1978 Practical Salinity Scale (PSS-78).For seawater with composition anomalies, estimates for the difference between Absolute and Reference Salinity are available for the global ocean and the Baltic Sea (McDougall et al., 2009;IOC et al., 2010;Feistel et al., 2010b).A routine for conversions between Practical Salinity and Absolute Salinity is included in the SIA library at level 0 (functions asal from psal and psal from asal) using an algorithm adapted from the GSW library (McDougall et al., 2009).A function that permits the computation of Density Salinity as an estimate of Absolute Salinity from measured in situ density, e.g., in the laboratory, is also available in the SIA library (function sea sa si).Further details are available from the companion paper, Part 2 (Wright et al., 2010a).Details on the use and the uncertainties of the different salinity scales are described by Wright et al. (2010b).Conversion routines between different pressure units (function cnv pressure), temperature scales and units (function cnv temperature), as well as between various salinity measures (function cnv salinity) are available at level 5 of the library for convenience.
In order to shorten the main article for technical reasons, 29 tables with groups of thermodynamic properties and their equations were moved to the Digital Supplement of this paper.They are referred to here as Table S1 to S29; their equations are numbered as S1.1 etc.Because of their key role in the SIA library, the various potential functions themselves are summarised in Table 1 rather than in the supplement.
The information provided in this paper is a reference for mathematical and thermodynamic details of the algorithms implemented in the SIA library.It is intended to ease the readability of the open source code and its numerous comment lines, and to encourage users to add new required properties, correlations or applications to the levels 2-5, guided by the examples and equations expained here.For high-speed requirements, tailored look-up tables for arbitrary property combinations can easily be compiled from the SIA code which is rather comprehensive but not speed-optimized.
While this paper was under review, the authors of the dry-air formulation (Lemmon et al., 2000) used in the SIA (5.78) + humid air g W F03 (T , P ) 5 fit liq g f03 si Fit of f F Liquid water IAPWS (2009c) g W IF97 (T , P ) 5 fit liq g if97 si Fit of f F Liquid water IAPWS (2007) g V IF97 (T , P ) 5 fit vap g if97 si Fit of f F Water vapour IAPWS (2007) library decided that the published molar equation can be converted to the mass-based form used here and in the planned IAPWS document (IAPWS, 2010), and that this should be implemented using the latest value for the molar mass of dry air (Picard et al., 2008) rather than the originally published one.For consistency with the IAPWS formulation, the molar mass of dry air of the SIA library is updated in the SIA version 1.1, attached as a supplement to the companion paper (Wright et al., 2010a), in contrast to the obsolete value used in SIA version 1.0 which is consistent with the formulation of Feistel et al. (2010a).
2 Level 1: Thermodynamic potentials -the primary standard As described in the related background articles, the vast amount of quantitative information available from extensive sets of experimental thermodynamic data for water, ice, seawater and air is represented in a compact way by the empirical coefficients of only four independent functions, a Helmholtz function f F (T ,ρ) of fluid water referred to as IAPWS-95 (IAPWS, 2009a;Wagner and Pruß, 2002), a Gibbs function g Ih (T ,P ) of ice, referred to as IAPWS-06 (IAPWS, 2009b;Feistel and Wagner, 2006), a Gibbs function g S (S A ,T ,P ) of sea salt dissolved in water which is referred to as IAPWS-08 (IAPWS, 2008a;Feistel, 2003Feistel, , 2008)), and a Helmholtz function f A (T ,ρ) for dry air    Magnified view of the small region corresponding to the standard oceanographic ("Neptunian") range.TP: triple point gas-liquid-solid, CP: critical point.The deviation of the vapourpressure line from the 101 325 Pa isobar in the liquid region is below the graphical resolution of panel (b).Freezing-point lowering occurs with the addition of sea salt.To deal with this effect in the case of seawater, the extension of the pure water properties into the metastable liquid region just above the line marked "Freezing Point Lowering" is required.al., 2000) in combination with air-water cross-virial coefficients (Hyland and Wexler, 1983;Harvey and Huang, 2007;Feistel et al., 2010a).These potential functions are used as the Primary Standard for pure water (liquid, vapour and solid), seawater and humid air from which all other properties are derived by mathematical operations, i.e. without the need for additional empirical functions.

Sea Salt Dissolved in Water
The Gibbs function (IAPWS, 2008a;Feistel, 2008) is expressed as the sum of a Gibbs function for pure water, ( ) , numerically available from the IAPWS-95 formulation, and a saline part, ( ) Here, salinity is expressed as Absolute Salinity SA, the mass fraction of dissolved salt in seawater, which for standard seawater equals the Reference-Composition Salinity within experimental uncertainty (Millero et al., 2008;Wright et al., 2010b).
In representing the properties of Standard Seawater, the range of validity of the Gibbs function for seawater is shown in Fig. 3.For temperatures in the oceanographic standard range, salinities up to 40 g/kg are properly described up to 100 MPa.For higher salinities up to 120 g/kg and temperatures up to 80 °C, the application is restricted to atmospheric pressure (101 325 Pa).Up to saturation, the salinity of cold concentrated brines agrees well with Antarctic sea-ice data (Fig. 6).For hot concentrates (region F in Fig. 3) the partial derivatives of density are not reliable.New density measurements (Millero and Huang, 2009;Safarov et al., 2009) have led to some recent improvements and an option to use an extension introduced

Fluid water
The validity range of the IAPWS-95 Helmholtz potential f F (T ,ρ) for fluid water (IAPWS, 2009a;Wagner and Pruß, 2002) as a function of temperature T and density ρ is shown in Fig. 1a in a density-temperature diagram.It is confined to the pressure interval between the isobars of 10 nPa and 1 GPa, below the upper temperature bound of 1000 • C and by the phase transition lines with ice and the liquid-vapour 2-phase region.Below the critical temperature, this region separates the stable vapour phase at low density from the stable liquid phase at high density.Only a small fraction of this region (a subset of the sliver to the right of the high density side of the phase transition boundary) belongs to the "Neptunian" oceanographic standard range (Fig. 1b).In the presence of dissolved sea salt, the freezing point is lowered so that the liquid phase of water is extended into the ice and vapour regions indicated in Fig. 1b.The Helmholtz function f F (T ,ρ) together with its first and second partial derivatives is implemented as the library function flu f si.

Ice
The IAPWS-06 Gibbs function g Ih (T ,P ) of hexagonal ice Ih (Feistel and Wagner, 2006;IAPWS, 2009b) covers the entire region of its stable existence (Fig. 2).In the region of low temperature and high pressure the function behaves reasonably although no experimental data were available when the function was constructed.Below 100 K, there are still open scientific questions regarding the possible phase transition to a proton-ordered ice XI or the existence of a density minimum.The Gibbs function is valid to even lower pressures (Feistel and Wagner, 2007) not shown here because the sublimation curve is restricted by the validity of the IAPWS-95 equation for vapour, Fig. 1a.In the library, an extension of by Feistel (2010) is included in the library (Wright et al., 2010a).Nevertheless, reliability of results for this region remains limited by the sparseness of data and the possibility of precipitation of calcium minerals (Marion et al., 2009)  The saline part ( ) of the Gibbs function together with its first and second partial derivatives is implemented as the library function sal_g_si.
Note that the arguments of the Gibbs function are temperature and pressure rather than temperature and density as in the Helmholtz function.Since the Gibbs function of pure water is expressed in terms of the corresponding Helmholtz function, sea_g_si is only available at library level 3 where implicitly determined quantities, such as density in terms of temperature and pressure, are considered.the vapour equation down to 50 K is implemented that permits the computation of sublimation properties to this limit (IAPWS, 2008c;Feistel et al., 2010a).Vapour cannot reasonably be expected to exist below 50 K (Feistel and Wagner, 2007).No ice forms other than Ih occur naturally under oceanographic conditions.
The Gibbs function g Ih (T ,P ) together with its first and second partial derivatives is implemented as the library function ice g si.

Sea salt dissolved in water
The Gibbs function g SW (S A ,T ,P ) of seawater (IAPWS, 2008a;Feistel, 2008) is expressed as the sum of a Gibbs function for pure water, g W (T ,P ), numerically available from the IAPWS-95 formulation, and a saline part, g S (S A ,T ,P ): g SW (S A ,T ,P ) = g W (T ,P ) + g S (S A ,T ,P ). (2.1) Here, salinity is expressed as Absolute Salinity S A , the mass fraction of dissolved salt in seawater, which for standard seawater equals the Reference-Composition Salinity within experimental uncertainty (Millero et al., 2008;Wright et al., 2010b).
In representing the properties of Standard Seawater, the range of validity of the Gibbs function for seawater is shown in Fig. 3.For temperatures in the oceanographic standard range, salinities up to 40 g/kg are properly described up to 100 MPa.For higher salinities up to 120 g/kg and temperatures up to 80 • C, the application is restricted to atmospheric pressure (101 325 Pa).Up to saturation, the salinity of cold concentrated brines agrees well with Antarctic sea-ice data (Fig. 6).For hot concentrates (region F in Fig. 3) the partial derivatives of density are not reliable.New density measurements (Millero and Huang, 2009;Safarov et al., 2009) have led to some recent improvements and an option to use an extension introduced by Feistel ( 2010) is included in the library (Wright et al., 2010a).Nevertheless, reliability of results for this region remains limited by the sparseness of data and the possibility of precipitation of calcium minerals (Marion et al., 2009) which would degrade the accuracy of the Reference Composition approximation in this region.
The saline part g S (S A ,T ,P ) of the Gibbs function together with its first and second partial derivatives is implemented as the library function sal g si.
Note that the arguments of the Gibbs function are temperature and pressure rather than temperature and density as in the Helmholtz function.Since the Gibbs function of pure water is expressed in terms of the corresponding Helmholtz function, sea g si is only available at library level 3 where implicitly determined quantities, such as density in terms of temperature and pressure, are considered.
The function g S (S A ,T ,P ) is constructed as a series expansion with respect to salinity.Based on the theory of ideal and electrolytic solutions (Planck, 1888;Landau and Lifschitz, 1964;Falkenhagen et al., 1971), this expansion consists of salinity-root and logarithmic terms and takes the form Here, the expansion coefficients are defined as with g u =1 J kg −1 , S u =35.16504 g kg −1 ×40/35, and the coefficients g ij k are given in the IAPWS Release 2008.The reduced temperature is τ =(T −T 0 )/T * , T 0 =273.15K, T * =40 K, the reduced pressure is π=(P −P 0 )/P * , P 0 =101 325 Pa, P * =10 8 Pa.The explicit separation of the expansion coefficients of Eq. (2.2) is required for the accurate determination of certain properties which in the zero-salinity limit possess a numerical singularity that can be analytically resolved.In some equations such as for the computation of potential temperature, one or more terms of the expansion (Eq.2.2) cancel analytically.Implementing the numerical solution of such an www.ocean-sci.net/6/633/2010/Ocean Sci., 6, 633-677, 2010 equation without the cancelling terms increases speed and accuracy.In such cases a function is more naturally (and easily) implemented by calling the separate functions (Eqs.2.3-2.5)rather than their combination in g S (S A ,T ,P ).
The expansion terms g i (T ,P ), Eqs.(2.3-2.5),together with their partial derivatives are available from the library function sal g term si.

Humid air
For a correct description of the thermodynamic properties at the ocean-atmosphere interface a thermodynamic potential of humid air is required and available from the literature (Feistel et al., 2010a).A related document is in preparation by IAPWS (2010).The Helmholtz function for dry air of Lemmon et al. (2000) has the form of the molar Helmholtz energy, f A,mol T ,ρ mol , depending on absolute temperature T (ITS-90) and molar air density, ρ mol .For its conversion to the specific Helmholtz energy, f A , depending on the mass density, ρ, the molar mass of air, M A =28.965 46 g mol −1 , is computed from the recent highly accurate air composition model of Picard et al. (2008).The dry-air part (Eq.2.6) can be combined with the vapour part, f V ≡f F (IAPWS-95, Sect.2.1), involving the second virial coefficient B AW (T ) of Harvey and Huang (2007) and the third virial coefficients C AAW (T ) and C AWW (T ) of air-vapour interaction reported by Hyland and Wexler (1983), to obtain the Helmholtz function of humid air, f AV , as Here, ρ is the density of humid air, A is the mass fraction of dry air in humid air, q=1−A is the specific humidity, (1−A)ρ is the absolute humidity, and r=(1−A)/A the humidity ratio or mixing ratio (van Wylen and Sonntag, 1965;Gill, 1982;Emanuel, 1994).R=8.314 51 J mol −1 K −1 is the molar (or universal) gas constant used by Lemmon et al. (2000), rather than the most recent value of R=8.314 472 J mol −1 K −1 (Mohr et al., 2008), and M W =0.018015268 kg mol −1 is the molar mass of pure water (IAPWS, 2008b).The effective molar mass of humid air M AV depends on the mass fraction A in the form The mass fraction A of air is computed from the mole fraction x A of dry air as and it follows that the mass fraction of vapour is given by .
The inverse function of Eq. (2.9), i.e., the mole fraction of air as a function of the mass fraction of air, is and the related mole fraction of vapour is .
The Helmholtz potential (Eq.2.7) is formally symmetric in the fractions of air and of water vapour.We note that the Helmholtz functions f V and f A that we have chosen to use in Eq. (2.7) are complete expressions rather than truncated expansions in terms of powers of density.Consequently, they include contributions corresponding to higher powers of density than included in the cross-virial terms represented by the third term in Eq. (2.13), ) is thus an inhomogeneous approximation formula with respect to the powers of density and the related correlation clusters.However, its validity is not restricted to small specific humidity, q=(1−A), such as some 1-3% often assumed for empirical equations used in meteorology.It can even be applied to physical situations in which air is the minor fraction, such as condensers of desalination plants or headspaces over subglacial lakes.The mass fraction A rather than the specific humidity q is chosen as a composition variable of humid air for its analogy to Absolute Salinity; the two describe the amount of natural mixtures, gases or salts, contained in ambient water in either gaseous or liquid form.This leads to thermodynamic equations that are formally similar in A and S A (Feistel et al., 2010a).
The range of validity is bounded by the simultaneous validity of the vapour formula (IAPWS-95), of the dry-air formula (Lemmon et al., 2000) and of the cross-virial expansion.The dry-air function correctly describes reliable experimental data for pressures up to 70 MPa and for temperatures from 60 to 873 K; the maximum air density in this region is 1035.8kg/m 3 .The temperature range where all three   (Hyland and Wexler 1983), not shown.
The air fraction is bound between 0 and 1 but is additionally limited by the vapour saturation condition, Fig. 5.At high total pressures, the restriction to vapour pressures below the saturation value represents a significant limitation on the upper limit of 1 -A that can be achieved in thermodynamic equilibrium.For total pressures below the vapour pressure of liquid water or the sublimation pressure of ice at the given temperature, the value of A may take any value between 0 and 1. (2.7).For oceanographic and meteorological applications it is unnecessary to consider liquid or solid air.Thus, we restrict consideration of Eq. (4.37) as follows: (i) for temperatures above the critical temperature of dry air, T >T c =132.5306 K, all density values occurring between the pressure bounds are permitted; and (ii) for subcritical temperatures T <T c , only densities below the dewpoint curve of dry air, indicated by "Dewpoint" are permitted.The resulting validity boundary for dry air is shown in bold."CP" is the critical point of dry air.To consider humid air, virial coefficients are required.The validity range in temperature of the third virial coefficients is shown by horizontal lines.Additionally, the pressure on saturated humid air is restricted to 5 MPa (Hyland and Wexler 1983), not shown.
virial coefficients are valid is from −80 to +200 • C, Fig. 4, (Hyland and Wexler, 1983).Consequently, the most limiting conditions for the validity of Eq. (2.7) are the temperature restrictions on the viral coefficients and the requirement for validity of the truncated virial expansion, i.e. the omitted terms of f AV proportional to A 3 (1−A)ρ 3 , A 2 (1−A) 2 ρ 3 and A(1−A) 3 ρ 3 must be negligibly small in comparison to the retained terms.A rough estimate for a maximum valid density is 100 kg m −3 as concluded from a comparison with experimental data for saturated air in which substantial fractions of both vapour and air are present (Feistel et al., 2010a;Fig. 8).When significant amounts of both air and water vapour are present, the valid temperature range is determined by the validity range for the virial coefficients.As the density of either the air or vapour component is decreased, the contribution from the virial coefficients decreases and the validity range in temperature extends to higher values, reaching 873 K when water vapour is eliminated and 1273 K when air is eliminated.The air fraction is bound between 0 and 1 but is additionally limited by the vapour saturation condition, Fig. 5.At high total pressures, the restriction to vapour pressures below the saturation value represents a significant limitation on the upper limit of 1−A that can be achieved in thermodynamic Saturation curves A sat (T ,P ) of humid air at the pressures 101.325, 50, 20 and 10 kPa, as indicated.Panel (a) shows results for temperatures above the freezing point, computed by solving Eq. (5.48) using the library function liq air massfraction air si, Eq. (S21.9), and panel (b) shows results for temperatures below the freezing point, computed by solving Eq. (5.70) using the function ice air massfraction air si, Eq. (S25.10).Valid air fraction values A are located above the particular saturation curve, A≥A sat (T ,P ), in the region indicated by "HU-MID AIR".In the presence of ice-free seawater, the validity range for A is more restricted, A≥A cond (S A ,T ,P )≥A sat (T ,P ), by the condensation value A cond , computed from the function sea air massfraction air si, Eq. (S29.1).equilibrium.For total pressures below the vapour pressure of liquid water or the sublimation pressure of ice at the given temperature, the value of A may take any value between 0 and 1.
The Helmholtz function f A (T ,ρ) for dry air together with its partial derivatives is implemented as the library function www.ocean-sci.net/6/633/2010/Ocean Sci., 6, 633-677, 2010 dry f si.The Helmholtz function for air-vapour interaction, together with its partial derivatives is implemented as the library function air f mix si.The cross-virial coefficients B AW , C AAW and C AWW are implemented as the library functions air baw m3mol, air caaw m6mol2 and air caww m6mol2.The Helmholtz function of humid air, f AV (A,T ,ρ), Eq. (2.7), together with its partial derivatives is implemented as the library function air f si.
For convenience of use, some auxiliary conversion functions, Eqs.(2.9-2.12), are also implemented at level 0, Table S1.Deviating from the original formulation given by Lemmon et al. (2000), in the library the adjustable constants of dry air are specified such that the entropy and the enthalpy of dry air are zero at the standard ocean state, T =273.15K and P =101325 Pa (Feistel et al., 2010a).This choice does not affect any measurable thermodynamic properties.

Level 2: Directly derived properties
From the level-one functions described in Sect.2, various thermodynamic properties can be computed directly if the corresponding independent variables are known.If some of the input variables need to be derived first from other known ones, based on thermodynamic relations, then the function will be found at level 3 (Sect.4) or higher.
The required input variables for level 2 functions are temperature and density of fluid pure water, either liquid or vapour (Sect.3.1), temperature and pressure for ice (Sect.3.2), and Absolute Salinity, temperature and pressure for dissolved sea salt (Sect.3.3).For moist air, level 2 routines require inputs of temperature, density and the mass fraction of (dry) air in the mixture.Specifying the air mass fraction as 1 gives the dry air limit.
The Jacobi method developed by Shaw (1935) is the mathematically most elegant way of transforming the various partial derivatives of different potential functions into each other, exploiting the convenient formal calculus of functional determinants (Margenau and Murphy, 1943;Landau and Lifschitz, 1964).Conversion tables (Feistel, 2008) between the potentials f (T ,ρ), g(S A ,T ,P ) and h(S A ,η,P ) are given in Sect. 5.

Fluid water
The total differential of the Helmholtz function f F (T ,ρ) of fluid water has the form where v=1/ρ is the specific volume.Therefore, the first derivatives of f F give the specific entropy, η, and the absolute pressure, P , A list of properties derived from f F (T ,ρ) by means of Eqs.(3.2) and (3.3) is given in Table S2.Partial derivatives with respect to these two independent variables are written as subscripts.Whether the property belongs to liquid water or vapour depends on the density used, i.e. on the location in the diagram in Fig. 1.

Ice
The total differential of the Gibbs function g Ih (T ,P ) of ice Ih has the form dg Ih = − ηdT + vdP .(3.6) A list of properties derived from g Ih (T ,P ) is given in Table S3.Partial derivatives with respect to the two independent variables are written as subscripts.

Dissolved sea salt
The total differential of the saline part g S (S A ,T ,P ) of the Gibbs function of seawater has the form (3.10) The list of properties derived from g S (S A ,T ,P ) is given in Table S4.Partial derivatives with respect to the three independent variables are written as subscripts where the subscript of S A is omitted for simplicity.Details on the definition of osmotic and activity coefficients are given by Falkenhagen et al. (1971), Millero and Leung (1976), Ewing et al. (1994), Lehmann et al. (1996), IUPAC (1997), Feistel and Marion (2007) and Feistel (2008).
The mean practical activity coefficient ln γ of sea salt (S4.1) can be computed from the activity potential ψ (S4.2) as (Feistel and Marion, 2007) Here, m=S A /[(1 − S A )×M S ] is the molality (moles of salt per kg of water) implemented in the library as sal molality si, and γ id =1 kg mol −1 is the asymptotic value of γ at infinite dilution.M S =31.4038218 g mol −1 is the mean molar mass of sea salt with Reference Composition (Millero et al., 2008), R=8.314 472 J mol −1 K −1 is the molar gas constant and (1−S A ) is the mass fraction of water in seawater.The zero-salinity limit of Eq. (3.11) is lim ln γ /γ id =0.The activity potential ψ (S A ,T ,P ), Eq. (S4.2), describes the ion-ion interactions and consists of higher salinity powers O S 3/2 A of the saline part of the Gibbs function (Eq.2.2) in the form (Feistel and Marion, 2007) g S (S A ,T ,P ) = S A g 2 (T ,P ) (3.12) Here, R S =R/M S =264.7599J kg −1 K −1 is the specific gas constant of sea salt.The activity potential is related to the osmotic coefficient φ and the activity coefficient ln γ by (3.13) The zero-salinity limit is lim ψ=0.The activity potential vanishes for ideal solutions.The osmotic coefficient φ , Eq. (S4.11), expresses the activity coefficient of water and can be computed from the activity potential ψ, Eq. (S4.2), as It is related to the chemical potential of pure water, g W (Sect. 4), and the chemical potential of water in seawater, µ W , by (Feistel and Marion, 2007) µ W (S A ,T ,P ) = g W (T ,P ) − mRT φ (S A ,T ,P ). (3.15) The zero-salinity limit is lim φ=1.
The saline excess chemical potential µ WS , Eq. (S4.3), is the difference between the chemical potentials of water in seawater and of pure water, µ WS (S A ,T ,P ) = µ W (S A ,T ,P ) − µ W (0,T ,P ) = − mRT φ. (3.16) The zero-salinity limit is lim The activity of water a w , Eq. (S4.3), is related to the osmotic coefficient by (3.17) Here, M W =18.015268 g mol −1 is the molar mass of water (IAPWS, 2008b) and R W =R/M W =461.523 64 J kg −1 K −1 is the specific gas constant of water.The zero-salinity limit is lim a W =1. At low vapour pressures, a W equals the relative humidity of sea air (Feistel et al., 2010a).
The relative chemical potential µ, Eq. (S4.5), describes the change of the Gibbs energy of a seawater parcel if at constant temperature and pressure a small mass fraction of water is replaced by salt.Its zero-salinity limit possesses a logarithmic singularity, lim The dilution coefficient D, Eq. (S4.6), describes the change of salinity in relation to freezing or evaporation processes, (Feistel and Hagen, 1998;Feistel et al., 2010a), as e.g. in Eqs.(A28), (4.44) or (A38).The zero-salinity limit (Raoult's law) is lim D=R S T .The chemical coefficient (S4.6), D S =S A D, is used for the description of sea air (Feistel et al., 2010a).
The specific enthalpy, entropy and volume of sea salt, Eqs.(S4.12)-(S4.14),provide the enthalpy, entropy and volume per mass of sea-salt particles dissolved in water.The zero-salinity limits are lim v S =(∂g 2 /∂P ) T .The logarithmic singularity of entropy reflects the empirical fact that rigorous purification of a mixture, i.e., complete desalination, is impossible by thermodynamic processes.
Mixing enthalpy, entropy and volume, Eqs.(S4.12)-(S4.14),provide the change of enthalpy, entropy or specific volume if two seawater samples with absolute salinities S 1 , S 2 and mass fractions w 1 , w 2 are mixed at constant temperature and pressure.If the mixing occurs adiabatically at constant pressure, the enthalpy remains constant while entropy is produced and the temperature changes.Since such effects do not occur in ideal solutions, the related quantities can be computed from the activity potential ψ (S A ,T ,P ) alone (Feistel and Marion, 2007).They disappear at infinite dilution.

Humid air
The Helmholtz function f AV (A,T ,ρ) of humid air, Eq. (2.7), permits the direct computation of all thermodynamic properties if temperature T , density ρ and air fraction A are either given or can be obtained from similar quantities such as the specific humidity, q=1−A or the mixing ratio r=(1−A)/A.This does not include properties at given relative humidity which requires the knowledge of vapour saturation, i.e. of the phase equilibrium between vapour and liquid water which is a composite system considered later in Sect.5.8.A list of equations for the computation of humidair properties from Eq. (2.7) is given in Table S5.

Level 3: Functions involving numerical solution of implicit equations
If quantities other than the natural independent variables of the three potential functions of Sect. 2 are given, in particular, if the pressure is known rather that the density of pure water, or the entropy rather than the temperature of seawater, the relevant thermodynamic equations must be inverted analytically or numerically.These steps inevitably add larger numerical uncertainties to all properties that depend on these inversions, and hence on the settings chosen for the associated iteration algorithms.Default values for iteration number or tolerance are specified in the SIA library routines that should be appropriate for most purposes; if necessary, they can be modified by related "set " procedures of the library (Wright et al., 2010a).Quantities that require such inversions appear in the libraries as level-3 procedures.To ensure the stability and uniqueness of the numerical solutions, initial conditions must be chosen appropriately.Various empirical functions are used to provide suitable initial values as discussed in the appendices referenced in Sect.4.1-4.3.While the algorithmic success and speed are sensitive to these choices, the final quantitative results are, within their numerical uncertainty, independent of the details of the initial "guess" functions.Therefore, if desired for certain applications, these auxiliary functions implemented in the library and described in this paper may be replaced by more suitable or effective customised ones without affecting the correctness of the final results.

Gibbs functions for liquid water and water vapour
To compute properties of fluid water at given T and P from its Helmholtz potential, f F (T ,ρ), it is necessary to solve Eq. (S2.11), for the density.Except for spurious or unstable numerical solutions outside the validity range, Fig. 1a, there is exactly one physically meaningful solution at supercritical temperatures.
Depending on the pressure, there may be one or two stable solutions below the critical temperature, given by the intersection points of isobars with isotherms illustrated in Fig. 1, providing the density of liquid water, ρ W (T ,P ), and/or of water vapour, ρ V (T ,P ).
Consequently, there cannot exist a single-valued Gibbs function g(T ,P ) that fully represents the properties of the Helmholtz function f F (T ,ρ) of fluid water.Rather, there are two different Gibbs functions, for liquid water and for water vapour, which coincide under supercritical conditions.Interestingly, critical conditions can be encountered at hydrothermal vents in the abyssal ocean (Reed, 2006;Sun et al., 2008).To implement the above expressions for the Gibbs functions we must determine the liquid and vapour densities corresponding to the temperature and pressure inputs.This requires iterative solution of Eq. (4.1), with considerable care required to select the appropriate root for each case.Details on the iterative numerical method and the conditions used to initialize the iteration procedure are provided in Appendix A1.
Once the liquid or vapour density of water is computed from the Helmholtz function f F at given temperature and pressure, the numerical values of the Gibbs function of water and its partial derivatives can be computed from the formulas of Table S6.The equations given there for water, g W , Eq. (4.2), apply in an analogous manner to vapour, g V , Eq. ( 4.3), if only the density ρ of liquid water is replaced by that of vapour.
Note that the above procedure is required to ensure arbitrarily precise consistency between the Gibbs function of pure water and the corresponding Helmholtz function.As long as this consistency is demanded, determination of the Gibbs function and its derivatives requires an iterative numerical procedure to determine the density argument of the Helmholtz function, so no explicit algebraic expression is possible.Thus, the pure water component of the Gibbs function must be determined at level 3 and it is only at this level that the Gibbs function for seawater can be completely determined.However, once the liquid pure water density is determined, the corresponding Gibbs potential is fully determined and it can be used in the seawater functions described in Sect.4.2 and 4.3.
Finally, note that the library functions listed in Table S2 for pure fluid water in terms of temperature and density are available as similar functions of temperature and pressure with the prefix liq for liquid water and vap for water vapour, respectively, rather than with the prefix flu given in Table S2.

Gibbs function of seawater
The Gibbs function of seawater, Eq. (2.1), is reproduced here as, g SW (S A ,T ,P ) = g W (T ,P ) + g S (S A ,T ,P ), (4.4) and is directly available from the sum of the Gibbs function of pure water computed at level 3, Table S6, and the saline part from the Primary Standard, level 1, Eq. (2.2).Properties of seawater can be computed from the partial derivatives of g S and g SW as given in Tables 4 and 7.
The Gibbs function g SW (S A ,T ,P ) of seawater, Eq. ( 4.4), together with its first and second partial derivatives is implemented as the library function sea g si.

Enthalpy of seawater
Besides the Gibbs and the Helmholtz functions, the specific enthalpy h SW (S A ,η,P ) of seawater, expressed in terms of Absolute Salinity S A , specific entropy, η, and absolute pressure P is a third important thermodynamic potential, useful in oceanography in particular for the computation of properties related to adiabatic processes (Feistel and Hagen, 1995;McDougall, 2003;Feistel, 2008;IOC et al., 2010).
To compute this potential and its partial derivatives from the Gibbs function g SW (S A ,T ,P ) of seawater, the independent variable T appearing in the expression for the enthalpy, to provide the implicit relation T =T (S A ,η,P ).Details on the iterative solution method used in the libraries are given in Appendix A2.
The specific enthalpy h SW (S A ,η,P ) of seawater, Eq. (4.5), as a thermodynamic potential is implemented in the library as the function sea h si.
Once the value of T has been determined as described in the appendix, the partial derivatives of h SW (S A ,η,P ) are obtained from those of g SW (S A ,T ,P ) as given in Table S8.
From the enthalpy and its derivatives, all thermodynamic properties can be computed.A selection is given in Table S9 and additional quantities are given in Table S10 after socalled "potential" properties are introduced.
Many oceanic processes like pressure excursions of a seawater parcel conserve salinity and entropy to very good approximation.In particular, if a parcel is moved this way to some reference pressure P =P r , the thermodynamic properties given in Table S9 can be computed at that reference level from the partial derivatives of h SW (S A ,η,P r ).Such properties derived from the potential function h SW at the reference pressure are commonly referred to as "potential" properties in meteorology and oceanography.Originally introduced by von Bezold (1888), potential temperature is defined as the temperature that a fluid parcel takes if it is moved adiabatically from its in situ pressure to a reference pressure level, which is often specified as the ocean surface.Analogous definitions hold for the potential density and potential enthalpy (IOC et al., 2010).
The potential enthalpy, h θ , is obtained from Eq. (S8.2), the absolute potential temperature, θ, in K, is obtained from Eq. (S9.2), and the potential density, ρ θ , is obtained from Eq. (S9.1), Evidently, for any fixed reference pressure, P r , the values of h SW (S A ,η,P r ) and its partial derivatives, as well as any other arbitrary function depending on this triple of variables, remain unchanged during isentropic (η = const) and isohaline (S A = const) processes.Derived from Eqs. (4.7) and (4.8), three kinds of thermal expansion and haline contraction coefficients are important for numerical models (IOC et al., 2010) and these are considered below as the cases (i) to (vi).In these cases, we have omitted the superscripts SW on the seawater potential functions for simplicity of the expressions.As well, we have always regarded the reference pressure P r as a constant value in each derivative considered here, without explicitly indicating this in the formulas.This implies that potential enthalpy, Eq. ( 4.7), and potential temperature, Eq. (4.8), are pressureindependent functions of salinity and entropy, and in particular, that any derivatives taken at constant (S A , η) can equivalently be taken at constant (S A , θ) or constant (S A , h θ ).An example is the isentropic compressibility, The thermal expansion coefficient, α T , is defined as: It is expressed in terms of derivatives of enthalpy by means of the Jacobi method and Eqs.(S9.1), (S9.2), as (∂v/∂η) S,P (∂T /∂η) S,P = h ηP h P h ηη .
Using Table S8, the partial derivatives of h can be substituted by those of g, with the result The thermal expansion coefficient with respect to potential temperature, α θ , is defined as: Similar to Eq. ( 4.12), with the help of Eq. ( 4.8) we compute Using Table S8, the partial derivatives of h can be substituted by those of g, with the result g P g T T . (4.16) Here, g θ is the potential Gibbs energy defined as g θ ≡ g(S A ,θ,P r ).
(iii) The thermal expansion coefficient with respect to potential enthalpy, α h , is defined as: Similar to Eq. (4.12), with the help of Eq. (4.7) we compute Using Table S8, the partial derivatives of h can be substituted by those of g, with the result The thermal expansion coefficient with respect to conservative temperature, α , is related to Eq. (4.19) by a constant conversion factor, c 0 IOC et al., 2010).Conservative temperature, , is potential specific enthalpy, h θ , Eq. (4.7), expressed in terms of an arbitrarily defined temperature unit, =h θ /c 0 P (McDougall, 2003;IOC et al., 2010); as such, it belongs to level 5 of the library where non-basic-SI units and userdefined functions are implemented.In contrast, potential enthalpy itself is defined at the core level 3 of the SIA library.
(iv) The isothermal haline contraction coefficient, β, is defined as: Similar to Eq. ( 4.12) we write Eq. (4.20) in terms of Jacobians (4.21) Expanding the functional determinant in the numerator yields, with the help of Eqs.(S9.1) and (S9.2) = − Using Table S8, the partial derivatives of h can be substituted by those of g, with the result The haline contraction coefficient with respect to potential temperature, β θ , is defined as: Similar to Eq. ( 4.12) we write Eq. (4.24) in terms of Jacobians Expanding the functional determinant in the numerator yields, with the help of Eqs.(S9.1) and (S9.2) Using Table S8, the partial derivatives of h can be substituted by those of g, with the result The haline contraction coefficient with respect to potential enthalpy, β , is defined as: Similar to Eq. ( 4.12) we write Eq. (4.28) in terms of Jacobians Expanding the functional determinant in the numerator yields, with the help of Eq. (4.7) Using Table S8, the partial derivatives of h can be substituted by those of g, with the result g ST g T P −g SP g T T −g θ S g T P /θ g P g T T . (4.31) The latter equalities in Eqs.(4.13), (4.16), (4.19), (4.23), (4.27) and (4.31) are the results given earlier in Table S7.
The potential quantities written in terms of the enthalpy of seawater are listed in Table S10.Entropy as a function of salinity, temperature and pressure is available from Eq. (S7.2).Potential temperature is defined by the relation η(S A ,T ,P )=η(S A ,θ,P r ), therefore the same function (Eq.S7.2) can be used to compute entropy as a function of salinity, potential temperature and reference pressure.Since the cases (i) to (vi) above, Eqs.(4.13), (4.16), (4.19), (4.23), (4.27) and (4.31), specify the different expansion and contraction coefficients as functions of entropy, these coefficients are available as functions of potential temperature, too, by means of Eq. (S7.2).
From the enthalpy definition Eq. (4.5) and the differential Eq. (3.7) of the Gibbs function, the relation can be inferred.Hence, when enthalpy is used as an independent thermal variable in combination with salinity and pressure, the responsible thermodynamic potential function is entropy, η SW (S A ,h,P ).Note that the superscript "SW" on η is included here to indicate its use as a thermodynamic potential function for seawater, consistent with the inclusion of "SW" on both the Gibbs function g SW and the enthalpy h SW when used as a potential function for seawater.To obtain this function value numerically from its arguments, the Eq.(S8.2) which are really the same functions with different arguments.The iterative inversion algorithm is straight forward and is implemented as the library function sea eta entropy si.It provides entropy η in the form of either η SW (S A ,h,P ) or η SW S A ,h θ ,P r , from which in turn all properties listed in Tables S9, S10 can be determined.Note, however, that we have not implemented an explicit routine for entropy, Eq. (4.35), as a potential function in the library.That is, the function sea eta entropy si provides entropy as a function of salinity, enthalpy and pressure, but it does not provide the partial derivatives of entropy with respect to those variables, nor does it take any orders of derivatives as input parameters.As such, the thermodynamic potential "entropy" is not available in the present SIA library version in the same form as the other potential functions that are summarised in Table 1.Nevertheless, various properties (Tables S9, S10) derived from it are implemented at level 3 and evaluated from indirect algorithms, just as if the potential "entropy" were available.The corresponding routines can be identified in the implementation of the library discussed in Part 2 (Wright et al., 2010a) by an eta instead of an h in the function names given in Table S10, which indicates the implicit use of entropy as the potential function.Consequently, these routines take enthalpy or potential enthalpy as the thermal input parameter rather than entropy.

Gibbs function of humid air
In many practical situations the pressure rather than the density of humid air is available from observations.For this purpose, the appropriate thermodynamic potential is the Gibbs function g AV (A,T ,P ) of humid air, computed from the Helmholtz function f AV (A,T ,ρ) of humid air, Eq. (2.7), by the Legendre transform (Alberty, 2001) and the subsequent substitution of the independent variable ρ by P , obtained from solving numerically the equation For oceanographic and meteorological applications it is not necessary to consider liquid or solid air.Therefore, we restrict consideration of Eq. (4.37) to the following regions: (i) for temperatures above the critical temperature of dry air, T >T c =132.5306 K, all pressures in the range shown in Fig. 4 are included; and (ii) at subcritical temperatures T <T c only temperatures higher than the dewpoint temperature T D , i.e. the condensation point of liquid air, are included.The function T D is available from Lemmon et al. (2000) and is shown in Fig. 4. As a starting value for the density iteration of Eq. ( 4.38) at given pressure P , air fraction A and temperature T , the ideal-gas equation is suitable: The molar mass of humid air M AV is given by Eq. (2.8), and R=8.314472 J mol −1 K −1 is the molar gas constant.Inserting the numerical result for ρ into Eq.(4.37) provides the required function value of g AV at given A, T , P .For the numerical computation of partial derivatives of the Gibbs function, algebraic combinations of analytical derivatives of the Helmholtz function are implemented as given in Table S11.Thermodynamic properties of humid air at given A, T , P are computed from the Gibbs function (Eq.4.37) and its partial derivatives, Table S11, as given in Table S12.The deviation of the compressibility factor Z AV from unity describes the non-ideal behaviour.The adiabatic lapse rate is given with respect to pressure rather than altitude and refers to subsaturated humid air, often referred to as "dry-adiabatic" in the meteorological literature.The air contraction coefficient, β=− 1 v ∂v ∂A T ,P , is the relative density increase if a small mass of vapour in a sample is replaced by air.Additional equations for humid-air properties are discussed in Feistel et al. (2010a).

Enthalpy of humid air
When humid air is lifted adiabatically from the surface to a certain pressure level, its air fraction and its entropy can often be considered as conservative during this process.Thus, the entropy η rather than the temperature T is known for a parcel at some given altitude if the initial entropy was computed at the surface.For this application purpose, the appropriate thermodynamic potential is the specific enthalpy h AV (A,η,P ) of humid air, computed from the Gibbs function g AV (A,T ,P ) of humid air, Eq. (4.37), by the Legendre transform (Alberty, 2001) The subsequent substitution of the independent variable T by η is obtained numerically from solving Eq. (S12.2) for T : As a starting value for the iterative solution of Eq. (4.41) for T at given pressure P , air fraction A and entropy η, we use an ideal-gas approximation of Eq. (S12.2) in the vicinity of the IAPWS-95 triple point (T t , P t ) of water: This expression does not depend on the particular choice made for the adjustable coefficients of the entropy.
The constants are η t =η(A,T t ,P t ) computed from Eq. (S12.2) at T t =273.16K, P t =611.654771Pa, The molar mass of air is M A =0.02896546 kg mol −1 , that of water is M W =0.018015268 kg mol −1 , M AV is given by Eq. (2.8), and R=8.314472 J mol −1 K −1 is the molar gas constant.
Once the value of T has been determined from solving Eq. (4.41), the partial derivatives of h AV (A,η,P ) are obtained from those of g AV (A,T ,P ) as given in Table S13.Thermodynamic properties as given in Table S14 can be calculated from algebraic combinations of these derivatives.
If a humid-air parcel is moved adiabatically to some reference pressure P =P r below its isentropic condensation level (ICL, Emanuel, 1994;Feistel et al., 2010a), all its thermodynamic properties given in Table S14 can be computed at that reference level from the partial derivatives of h AV (A,η,P r ) and the in situ entropy η(A,T ,P ), Eq. (S13.1).As discussed for seawater in Sect.4.3, properties derived from the potential function h AV at the reference pressure (frequently specified as the surface pressure) are commonly referred to as "potential" properties in meteorology.Examples are the potential enthalpy, h θ , library function air potenthalpy si, the potential temperature from Eq. (S14.2),θ, in • C, library function air pottemp si, and the potential density from Eq. (S14.1),ρ θ , library function air potdensity si, Evidently, for any fixed reference pressure, P r , the value of h AV (A,η,P r ) and its partial derivatives, as well as any other arbitrary function depending on this triple of variables, remain unchanged during isentropic processes (η = const) at constant specific humidity (A = const).
Physically reasonable values of the entropy to be used as an independent variable of the enthalpy are restricted to ranges depending on humidity and pressure, between the particular limits given by dry and saturated air, see Sect.5.8.

Level 4: Phase equilibria and composite systems
Equilibrium properties at phase transition boundaries or of coexisting phases are often characterized by drastic spatial or temporal changes, and large values of latent heat exchange or volume expansion, e.g. if seawater freezes or evaporates.Such multi-phase and multi-component properties are available from combinations of the thermodynamic potentials if they are consistently adjusted to reference state conditions which fix the absolute energies and entropies of the substances involved (Feistel et al., 2008).Gibbs functions can be constructed for composite systems such as sea ice (Feistel and Hagen, 1998, Sect. 5.4) that contain two stable phases (e.g.ice and seawater).When the temperature, the volume or the pressure of a composite system is changed, mass is transferred from one phase to the other; for example if sea water freezes or evaporates, brine salinity or vapour pressure adjust to the new conditions imposed and the heat capacity or the thermal expansion of the whole system exhibits very large changes resulting from the changes in latent heat contributions.By utilizing mutually consistent potential functions, rigorous mathematical formulae can be determined for the numerical calculation of latent properties depending on the particular conditions such as isobaric, isochoric or isentropic processes.

Equilibrium liquid water-vapour: saturation
The saturation point of pure water is usually computed at a given temperature T , providing the vapour pressure P = P vap (T ), or at a given pressure P providing the boiling temperature T = T boil (P ).The defining condition is equality of the chemical potentials of liquid and vapour, which equal the Gibbs functions in the case of pure phases, g W (T ,P ) = g V (T ,P ). (5.1) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.1) is expressed by the system which exploits the relations (Eqs.S2.6 and S2.11) to avoid stacked numerical iterations.Eq. ( 5.2) is equivalent to Eq. ( 5.1) and is also known as the "Maxwell condition" in the form . Equations (5.2)-(5.4)provide three equations for the four unknowns T , P , ρ V and ρ W .Any one of these quantities can be specified independently to complete the system and permit the numerical solution as discussed in Appendix A3.
Once the values of T , P , ρ V and ρ W are computed from the iteration of Eqs.(5.2)-( 5.4) at the specified saturation condition, various equilibrium properties can be determined from the formulae given in Table S15.

Equilibrium water-ice: melting and freezing
The melting pressure of ice is usually computed at a given temperature T , giving P melt (T ).Similarly, the freezing temperature of water is normally determined at a given pressure P , giving T frz (P ), which also gives the melting temperature of ice.In either case, the defining condition is equality of the chemical potentials of liquid water and ice, g W (T ,P ) = g Ih (T ,P ). (5.5) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.5) is represented as the system which exploits the relations (Eqs.S2.6 and ES2.11) to avoid stacked numerical iterations.Equations (5.6) and (5.7) supply two equations for the three unknowns T , P and ρ W . Specifying any one of these quantities completes the determination of the system which can then be solved as discussed in Appendix A4.
Once the values of T , P and ρ W are computed from the iteration of Eqs.(5.6), (5.7) at the specified melting condition, various equilibrium properties can be determined from the formulae given in Table S16.

Sublimation equilibrium ice-vapour
The sublimation pressure of ice is usually computed at a given temperature T , giving P subl (T ).Similarly, the sublimation temperature of ice is usually computed at a given pressure P , giving T subl (P ), which also gives the ice-point temperature of vapour at which frost is formed.The defining condition is equality of the chemical potentials of water vapour and ice, g V (T ,P ) = g Ih (T ,P ). (5.8) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.8) is represented by the system (5.10) which exploits the relations (Eqs.S2.6 and S2.11) to avoid stacked numerical iterations.Specifying any one of T , P and ρ V , completes the system and allows numerical solution as discussed in Appendix A5.
Once the values of T , P and ρ V are computed from the iteration of Eqs. ( 59) and (5.10) at the given sublimation condition, various equilibrium properties can be determined from the formulae given in Table S17.

Equilibrium seawater-ice: sea ice
The freezing temperature of seawater with absolute salinity S A is usually computed at a specified pressure P giving T frz (S A , P ).Similarly, the brine salinity of sea ice is calculated at a specified temperature T and pressure P giving S brine A (T ,P ), and the melting pressure at which the solid fraction of sea ice disappears is calculated at specified S A and T giving P melt (S A ,T ).The defining condition for each of these is equality of the chemical potentials of ice and of water in seawater, g Ih (T ,P ) = g SW (S A ,T ,P ) − S A ∂g SW ∂S A T ,P . (5.11) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.11) is represented as the system

S
(5.12) (5.13) which exploits the relations Eqs.(S2.6), (S2.11), (4.4) and (S7.12) to avoid stacked numerical iterations.The function f F T ,ρ W is abbreviated here by f W , and similarly for its partial derivatives.Equations (5.12) and (5.13) provide two conditions for the four unknowns S A , T , P and ρ W .To complete the system, two of these variables must be specified, usually out of the triple S A , T or P .Once two of these variables are specified, the system may be solved iteratively as discussed in Appendix A6.
Once the values of S A , T , P and ρ W are computed from the iteration of Eqs.(5.12) and (5.13) for the specified choice of parameters, various single-phase equilibrium properties can be determined from the formulae given in Table S18.
The composite system "sea ice" consisting of seawater and ice can be described by a suitable Gibbs function g SI (S SI , T , P ) which is available from the equilibrium solution of Eq. (5.11) and can be used to compute all thermodynamic properties of this two-phase system, in particular its latent heat (for details see Feistel and Hagen, 1998): g SI (S SI ,T ,P )=(1−b)g Ih (T ,P )+bg SW (S A ,T ,P ).(5.14) Here, b=S SI /S A ≤1 is the mass fraction of brine and S SI is the given "bulk" or sea-ice salinity, i.e. the mass fraction of salt in sea ice, in contrast to the brine salinity S A (T ,P ), the mass fraction of salt in the liquid part, which is a function of temperature and pressure controlled by the equilibrium Eq. (5.11).For a compact writing of the partial derivatives of Eq. (5.14) it is useful to define a formal latency operator of sea ice, (5.15) Here, z is a certain thermodynamic function.For example, the equilibrium condition (Eq.5.11) can be written in the form SI [g] = 0. (5.16) The total differential of Eq. (5.16) is commonly known as the Clausius-Clapeyron differential equation of this phase transition: The first term yields the chemical coefficient, Eq. (S4.6), (5.18) which has a positive sign as can be concluded from the Second Law of Thermodynamics (Landau and Lifschitz, 1964;IOC et al., 2010).From the second and third terms of Eq. ( 5.17) we infer the derivatives of the brine salinity to be (5.20) With the help of these relations we can compute thermodynamic properties of sea ice from the partial derivatives of the Gibbs function (Eq.5.14) as given in Table S7 for seawater if the salinity S A considered there is substituted by the sea-ice salinity S SI and the Gibbs function g SW (S A ,T ,P ) by g SI (S SI ,T ,P ).The Gibbs function of sea ice, Eq. (5.14), is implemented as the function sea ice g si in the library.Related to the intrinsic phase transition, certain properties of interest are very specific for a composite system like sea ice and not listed in Table S7.Using Eq. (5.19), the isobaric melting rate, i.e. the increase of the brine fraction b=S SI /S A upon warming is (5.21) The isobaric heat capacity of sea ice computed from Eqs. (5.14), (S7.6) and (5.19), S SI ,P (5.22) consists of the single-phase contributions of ice, c Ih P , and brine, c SW P , as well as a latent part, and is implemented as sea ice cp seaice si.In the latter term, the coefficient L SI P in front of the melting rate, Eq. (5.21), is the isobaric latent heat of sea ice, and is available from the function sea ice enthalpy melt si in the library.The enthalpy of the brine, h SW , is computed from Eq. (S7.3), and the enthalpy of ice, h Ih , is computed from Eq. (S3.4).
Similarly, the thermal expansion of sea ice is computed from Eqs. and is available from the function sea ice expansion seaice si in the library.The third term on the right hand side of Eq. (5.24) is the melting rate, Eq. (5.21), multiplied by the isobaric melting volume of sea ice, function sea ice volume melt si, (5.25) The specific volume of the brine, v SW , is computed from Eq. (S7.1), and that of ice, v Ih , from Eq. (S3.13).
Since the freezing-point lowering due to pressure always exceeds the adiabatic lapserate of seawater, cold seawater may freeze and decompose into ice and brine during adiabatic uplift but this can never happen to a sinking parcel.This freezing process can destabilize the water column, e.g. off the Antarctic shelf (Foldvik and Kvinge, 1974), since the thermal expansion of sea ice, α SI =g SI T P /g SI P , Eq. ( 5.24), function sea ice expansion seaice si in the library, and consequently also the adiabatic lapserate (McDougall and Feistel, 2003) of sea ice, SI =−g SI T P /g SI T T , Eq. (S18.14),possess large negative values near the freezing point (Feistel and Hagen, 1998).These and related properties can be evaluated directly from the partial derivatives of the Gibbs function of sea ice, Eq. (5.14), implemented as the function sea ice g si in the library, in terms of the in situ temperature.For a seawater parcel, the potential temperature that corresponds to the freezing point under pressure is somewhat ill-defined physically since it is practically impossible to lift a parcel at the freezing point to the surface isentropically without decomposition into ice and brine.Freezing of a seawater parcel cannot occur at any depth as long as its potential temperature referenced to the surface is higher than its freezing point temperature T frz (S A , P SO ) computed from Eq. (5.11) at the surface pressure P SO , as discussed by Jackett et al. (2006).
As an observational example, the brine salinity, Eq. (S18.2), of Antarctic sea ice at normal pressure is shown in Fig. 6 in comparison to measurements of concentrated brines by Gleitz et al. (1995) and Fischer (2009).Note that only freezing point measurements at salinities less than 40 g kg −1 were used for the construction of IAPWS-08 (Feistel, 2008).

Equilibrium seawater-vapour
The vapour pressure of seawater is usually computed as a function of temperature T and Absolute Salinity S A , givingP vap (S A ,T ).Similarly, the boiling temperature of seawater is computed as a function of absolute pressure P and Absolute Salinity S A , giving T boil (S A ,P ), and the equilibrium brine salinity is computed as a function of T and P , giving S brine A (T ,P ).The defining condition for each of these quantities is equality of the chemical potentials of water www.ocean-sci.net/6/633/2010/Ocean Sci., 6, 633-677, 2010 Brine Salinity of Antarctic Sea Ice Brine salinity computed from Eq. (S18.2) at given temperature and normal pressure, compared with measured results for Antarctic sea ice.Symbol "F": data of Fischer (2009), "G": data of Gleitz et al. (1995).

Equilibrium Seawater-Vapour
The vapour pressure of seawater is usually computed as a function of temperature T and Absolute Salinity SA, giving P vap (SA, T).Similarly, the boiling temperature of seawater is computed as a function of absolute pressure P and Absolute Salinity SA, giving T boil (SA, P), and the equilibrium brine salinity is computed as a function of T and P, giving ) , ( brine A P T S .The defining condition for each of these quantities is equality of the chemical potentials of water vapour and of water in seawater, ( ) ( ) (5.26) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq.
vapour and of water in seawater, g V (T ,P ) = g SW (S A ,T ,P ) − S A ∂g SW ∂S A T ,P . (5.26) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.26) is expressed as the system

S
(5.27) which exploits the relations (Eqs.S2.6, S2.11, 4.4 and S7.12) to avoid stacked numerical iterations.The function f F T ,ρ V is abbreviated here by f V , and similarly for f W and their partial derivatives.Equations (5.27), (5.28), (5.29) provide three conditions for the five unknowns S A , T , P , ρ V and ρ W , so two of these parameters must be specified, usually from the set S A , T , P , to complete the system.Once this choice is made, the system can be solved as discussed in Appendix A7.
Once the values of S A , T , P , ρ V and ρ W are computed from the iteration of Eqs.(5.27)-(5.29) at the given evaporation conditions, various equilibrium properties can be determined from the formulae given in Table S19.
We use the name "sea vapour" to refer to a composite system consisting of seawater and vapour in thermodynamic equilibrium.Its Gibbs function g SV (S SV ,T ,P ) depends on absolute temperature T , absolute pressure P and the mass fraction of salt, which is the "bulk" or "sea-vapour" salinity S SV ; the Gibbs function is expressed as g SV (S SV ,T ,P ) (5.30) = (1 − b)g V (T ,P ) + bg SW (S A ,T ,P ).
Here, b=S SV /S A ≤1 is the mass fraction of brine.The brine salinity S A (T ,P ) is a function of temperature and pressure controlled by the equilibrium Eq. (5.26).For a compact writing of the partial derivatives of Eq. (5.30) it is useful to define a formal latency operator of sea vapour, (5.31) Here, z is a certain thermodynamic function.For example, the equilibrium condition Eq. ( 5.26) can be written in the form SV [g] = 0. (5.32) The total differential of Eq. (5.32) is commonly known as the Clausius-Clapeyron differential equation of this phase transition: The first term is the chemical coefficient (Eq.S4.6), (5.34) From the second and third terms of Eq. (5.32) we infer the derivatives of the brine salinity, (5.36) With the help of these relations we can compute thermodynamic properties of sea vapour from the partial derivatives of the Gibbs function (Eq.5.30) as given in Table S7 for seawater if the salinity S A considered there is substituted by the sea-vapour salinity S SV and the Gibbs function g SW (S A ,T ,P ) by g SV (S SV ,T ,P ).The Gibbs function of sea vapour, Eq. (5.30), is implemented as the function sea vap g si in the library.
Certain properties of interest are very specific for a composite system and not listed in Table S7.Using Eq. (5.35), the isobaric evaporation rate, i.e. the decrease of the brine fraction b=S SV /S A upon warming is (5.37) The isobaric heat capacity of sea vapour computed from Eqs. (5.30), (S7.6) and (5.35), S SV ,P (5.38) consists of the single-phase contributions of vapour, c V P , and brine, c SW P , as well as a latent part, and is implemented as sea vap cp seavap si.In the latter term, the coefficient in front of the evaporation rate (Eq.5.37) is the isobaric latent heat L SV P of seawater, which is available from the function sea vap enthalpy evap si in the library.The brine enthalpy, h SW , is computed from Eq. (S7.3), and the vapour enthalpy, h V , from Eq. (S2.3).

Osmotic equilibrium seawater-liquid
If pure water is separated from seawater by a semi-permeable membrane which lets water molecules pass but not salt particles, water will penetrate into the seawater, this way diluting it and possibly increasing its pressure, until the chemical potential of water in both boxes will be the same (or the pure water reservoir is exhausted).In the usual model configuration, the two samples are thermally coupled but may possess different pressures; the resulting pressure difference required to maintain equilibrium is the osmotic pressure of seawater.An example is desalination by reverse osmosis; if the pressure on seawater in a vessel exceeds its osmotic pressure, freshwater can be squeezed out of the solution through suitable membrane walls (Sherwood et al., 1967).The defining condition for the osmotic equilibrium is equality of the chemical potentials of pure water at the pressure P W and of water in seawater at the pressure P S , g W T ,P W =g SW S A ,T ,P S −S A ∂g SW ∂S A T ,P S .(5.40) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.40) is expressed as the system (5.43) which exploits the relations (Eqs.S2.6, S2.11, 4.4 and S7.12) to avoid stacked numerical iterations.
The function f F T ,ρ W is abbreviated here by f W , and similarly for f S computed at the liquid-water density ρ S related to the pressure P S , as well as their partial derivatives.Equations (5.41), (5.42), (5.43) provide three conditions for the six unknownsS A , T , P S , P W , ρ S and ρ W , so three of these parameters must be specified to complete the system.Once this choice is made, the remaining parameters can be determined as in Appendix A8.
Once the solution for S A , T , P S , P W , ρ S and ρ W has been found, the desired properties of the equilibrium can be computed, in particular the osmotic pressure, P osm =P S −P W .The related function sea liq osmoticpressure si is implemented in the library.

Triple point sea ice -vapour
The equilibrium between sea ice and vapour includes three phases, solid, liquid and gas, and two components, water and salt.Air is not involved.This equilibrium state extends the ordinary triple point of pure water to non-zero salinities, i.e. along a one-dimensional manifold.This curve is shown in Fig. 3 by the "Triple Line" which has the same T −P relation as the sublimation line because ice is in sublimation equilibrium with water vapour at any given brine salinity.Note that saturation is defined as the equilibrium state between water vapour and liquid water above the freezing point of pure water, or, below that temperature, between water vapour and ice (IAPWS, 2010).Hence, as soon as ice is present in an equilibrium system, the water vapour in the gas phase is regarded as saturated.
The equilibrium conditions for temperature, pressure and chemical potentials that determine the locus of triple points are expressed in terms of the Primary Standard as

S
(5.44) Any one of the five parameters may be specified to complete the system which may then be solved as discussed in Appendix A.9.
If any one of the three variables SA, T, P is specified, the other two are determined by the above conditions.Fig. 7 shows the displacement of the triple point along the sublimation line as a function of salinity.In the library, the equilibrium properties P, T and SA of sea-ice vapour are available from the functions sea_ice_vap_pressure_si, sea_ice_vap_temperature_si and sea_ice_vap_salinity_si.Note that the equilibrium conditions are actually determined by Fig. 7. Temperature-pressure phase diagram of seawater in the vicinity of the triple point.At different salinities, the triple point (TP), i.e. the equilibrium between liquid seawater, ice and vapour is displaced along the sublimation line (in bold) of the ice-vapour equilibrium.Note that the triple-point pressure can change by a factor of 2 while the vapour-pressure lowering at constant temperature is only of order 10% or less.
(5.47) Equations ( 5.44)-(5.47)provide four conditions for the five unknowns S A , T , P , ρ V and ρ W .Any one of the five parameters may be specified to complete the system which may then be solved as discussed in Appendix A9.
If any one of the three variables S A , T , P is specified, the other two are determined by the above conditions.Figure 7 shows the displacement of the triple point along the sublimation line as a function of salinity.
In the library, the equilibrium properties P , T and S A of sea-ice vapour are available from the functions sea ice vap pressure si, sea ice vap temperature si and sea ice vap salinity si.Note that the equilibrium conditions are actually determined by calling one of set sea ice vap eq at p, set sea ice vap eq at t or set sea ice vap eq at s, depending on which of pressure, temperature or salinity is specified.Thus, one of these "set "-routines must be called before accessing P , T or S A using the above function calls, but all three equilibrium properties corresponding to the specified parameter choice are available once the appropriate "set "-routine is executed.

Equilibrium humid air -liquid water
The state in which humid air is in equilibrium with liquid water is commonly referred to as "saturated air", the "dewpoint" or the "condensation point".The condition for this equilibrium is equal chemical potentials of liquid water, Eq. (4.2), and of water in humid air, Eq. (S12.15), (5.48) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.48) is expressed using the relations g W (T ,P ) = f F T ,ρ W + P /ρ W (5.49) g AV = f AV A,T ,ρ AV + P /ρ AV (5.51) . (5.52) The independent variables in this scheme are the total pressure, P , the liquid density, ρ W , the humid-air density, ρ AV , the temperature, T , and the air fraction, A. Using Eqs.(5.49) and (5.51) to eliminate the Gibbs potential in favour of the Helmholtz potentials results in three equations for these five unknowns.
For the numerical solution, two of the five unknowns as well as starting values for the remaining unknowns must be specified.Four important cases are considered in detail in Appendix A10.
No matter which of the four cases considered in the appendix is applied to compute the equilibrium between liquid water and humid air, the numerical solution of Eqs.(5.48)-(5.52)provides a consistent set of equilibrium values for A, T ,P , ρ W and ρ AV which is then available for the computation of any other property of either saturated humid air or liquid water in this state.
For example, at given temperature T and total pressure P , the partial vapour pressure of saturated air is available in the form from the solution obtained for A (T , P ), using library function liq air massfraction air si and then converting to the mole fraction of vapour, x AV V =1−x AV A (A), Eq. (S1.5), using the library function air molfraction vap si.The comparison with experimental data for the saturated vapour pressure (Feistel et al., 2010a), Fig. 8, permits an estimate of the effect of the cross-virial coefficients B AW (T ), C AAW (T ) and C AWW (T ) on the validity of this formulation at higher densities.From the data scatter at low densities we may estimate the experimental uncertainty to be better than 0.5%.Ambient air has a density of typically 1 kg m −3 or less.At about 10 kg m −3 , an error of about 2% must be expected if the second virial coefficients B AW (T ) is omitted.The same error occurs at 100 kg m −3 with B AW (T ) included but is reduced to 1% if the coefficients C AAW (T ) and C AWW (T ) are considered, too.
A practically important quantity is the relative humidity, RH, which expresses the deviation of the air fraction A of a given sample of humid air from the saturation value A sat (T ,P ) belonging to the same temperature and pressure, computed as the solution of Eq. (5.48) in the scenario of case 3 from Appendix A10.Out of several options available from the literature, two different definitions are implemented and attributed here to the WMO 3 and to the CCT 4 , RH WMO (A,T ,P ) = 1/A − 1 1/A sat (T ,P ) − 1 (5.54) and RH CCT (A,T ,P ) = x AV V (A) x AV V A sat (T ,P ) .
(5.55) with the mole fraction x AV V =1−x AV A (A) from Eq. (2.12).According to Jacobson (2005), the WMO defines the relative humidity as given in Eq. (5.54).This definition is also given by other sources such as Gill (1982).Alternatively, internal discussion documents of BIPM CCT-WG6 (Jeremy Lovell-Smith, private communication, 2010) consider as a suitable option for the definition of relative humidity the commonly used formula (Eq.5.55).This definition is also recommended in a recent document of WMO (2008), in contrast to Eq. (5.54).The definition of relative humidity given by the International Union of Pure and Applied Chemistry (IUPAC, 1997) is similar to Eq. (5.55) but uses the ratio of the partial pressure of water vapour in humid air to the pressure of saturated, air-free vapour, and does not exactly match 100% at saturation.
In the library, the conversion functions from air fraction to relative humidity are implemented as liq air rh wmo from a si and liq air rh cct from a si.Their inverse functions are and, from Eqs. (2.9) and (2.11), . (5.57) 3 WMO: World Meteorological Organisation, www.wmo.int 4 CCT: Consultative Committee for Thermometry, www.bipm.org/en/committees/cc/cct/cross-virial coefficients BAW(T), CAAW(T) and CAWW(T) on the validity of this formulation at higher densities.From the data scatter at low densities we may estimate the experimental uncertainty to be better than 0.5%.Ambient air has a density of typically 1 kg m -3 or less.At about 10 kg m -3 , an error of about 2% must be expected if the second virial coefficients BAW(T) is omitted.The same error occurs at 100 kg m -3 with BAW(T) included but is reduced to 1% if the coefficients CAAW(T) and CAWW(T) are considered, too.
Fig. 8 Experimental data for the saturated vapour pressure exp sat,

P
of humid air at different pressures P and temperatures T as reported in (Feistel et al., 2010a), in comparison to calc sat, P computed from Eqs. (5.53), (S21.9) and (S1.5).Symbol "o": formula (2.7) without cross-virial coefficients, "B": formula with the second crossvirial coefficient BAW(T), "C": formula with the second and third cross-virial coefficients BAW(T), CAAW(T), CAWW(T).The smaller scatter is magnified in panel b).The improvement realized by including the C coefficients is effective mainly at densities higher than 100 kg m -3 .
A practically important quantity is the relative humidity, RH, which expresses the deviation of the air fraction A of a given sample of humid air from the saturation value A sat (T, P) belonging to the same temperature and pressure, computed as the solution of Eq. (5.48) in the scenario of case 3 from Appendix A.10.Out of several options available from the literature, two different definitions are implemented and attributed here to the WMO 3 and to the CCT 4 , (5.54) and 3 WMO: World Meteorological Organisation, www.wmo.int 4 CCT: Consultative Committee for Thermometry, www.bipm.org/en/committees/cc/cct/Fig. 8. Experimental data for the saturated vapour pressure P sat,exp of humid air at different pressures P and temperatures T as reported in (Feistel et al., 2010a), in comparison to P sat,calc computed from Eqs. (5.53), (S21.9) and (S1.5).Symbol "o": formula Eq. (2.7) without cross-virial coefficients, "B": formula with the second cross-virial coefficient B AW (T ), "C": formula with the second and third cross-virial coefficients B AW (T ), C AAW (T ), C AWW (T ).The smaller scatter is magnified in panel (b).The improvement realized by including the C coefficients is effective mainly at densities higher than 100 kg m −3 .
They are implemented in the library as liq air a from rh wmo si and liq air a from rh cct si.
With "wet air" we refer to a composite system of liquid water and humid air with the mass fractions w A of dry air, w V of vapour and w W of liquid water, w A +w V +w W =1. The mutual equilibrium requires A sat (T ,P )=w A /w AV , with w AV =w A +w V =1−w W being the gaseous mass fraction.The Gibbs function of wet air reads (Feistel et al., 2010a) www.ocean-sci.net/6/633/2010/Ocean Sci., 6, 633-677, 2010 g AW w A ,T ,P = w A A sat (T ,P ) g AV A sat ,T ,P (5.58) and is a linear function of the air fraction, w A .Various wet-air properties are available from combinations of partial derivatives of the potential (Eq.5.58) with A sat (T ,P ) computed from Eqs. (5.48)-(5.52)as in case 3 from Appendix A10, see Table S21.For the computation of the partial T −P derivatives of g AW , the first derivatives of A sat (T ,P ) are required.Taking the respective derivatives of Eq. ( 5.48) we get the isobaric drying rate, and the isothermal drying rate, of humid air, i.e. the decrease of its saturated air fraction A sat due to heating or compression.The chemical coefficient D A is defined in Eq. (S12.16).The latency operator AW of wet air used here is defined for the specific entropy, η AW =−g AW T , of the form (5.61) and for the specific volume v AW =g AW P of the form (5.62) The partial derivatives of the Gibbs function g AW w A ,T ,P , Eq. (5.58), of wet air are given in Table S20.Properties of wet air computed from this Gibbs function are given in Table S21.
For the description of isentropic processes such as the uplift of wet air in the atmosphere, enthalpy h AW w A ,η,P computed from the Gibbs function (Eq.5.58) is a useful thermodynamic potential: (5.63) For this purpose, the temperature T corresponding to a given entropy η, must be determined to evaluate the right side of Eq. (5.63).The appropriate value of T must be obtained by numerically solving the equation (5.64)This is not a trivial task because a good analytical estimate for T AW w A ,η,P is not available, and the Newton iteration of Eq. ( 5.64) tends to be unstable so that the range of starting parameters that yields convergent solutions of Eq. (5.64) is rather restricted.An interval method like Brent's algorithm appeared to be the best choice in this case, applied between upper and lower temperature bounds.These limits follow from the physical conditions that wet air can only exist between freezing and complete evaporation of the liquid water part.Thus, the lower temperature bound T min (w A ,P ) for the solution of Eq. (5.64) is the freezing temperature T melt (P ) of water under the pressure P , g W (T min ,P ) = g Ih (T min ,P ), (5.65) computed from Eq. (5.5).The upper temperature bound T max (w A ,P ) is computed from Eq. ( 5.48) in case 2 from Appendix A10 for a vanishing liquid fraction, Eq. (S21.9),i.e., the air fraction A of humid air equal to that of wet air, w A =A: The entropy range corresponding to the interval T min -T max is shown in Fig. 9.The partial derivatives of the enthalpy h AW w A ,η,P are computed from those of the Gibbs function, Table S20, as given in Table S22.
Selected properties of wet air computed from the enthalpy (Eq.5.63) and its partial derivatives are given in Table S23.
Many meteorological processes such as adiabatic uplift of a wet-air parcel conserve specific humidity and entropy to very good approximation.In particular, if a parcel is moved this way to some reference pressure P =P r , all its thermodynamic properties given in Table S23 can be computed at that reference level from the partial derivatives of h AW w A ,η,P r .Such properties derived from the potential function h AW at the reference pressure are commonly referred to as "potential" properties in meteorology (von Bezold 1888, von Helmholtz 1888).Examples are the potential enthalpy, h θ , (5.67) the potential temperature, θ, in • C, obtained from Eq. (S23.2) .68)and the potential density, ρ θ , obtained from Eq. (S23.1) (5.69) The partial derivatives of the enthalpy ( ) A AW η are computed from those of the Gibbs function, Table S20, as given in Table S22.

M P a M P a M P a M P a T ri p le L in e T ri p le L in e T ri p le L in e T ri p le L in e Ice Format ion Ice Format ion Ice Format ion Ice Format ion D e w p o in t D e w p o in t D e w p o in t D e w p o in t WET AIR WET AIR WET AIR WET AIR HUMID AIR HUMID AIR HUMID AIR HUMID AIR
Fig. 9 Valid entropy values of wet air computed from Eq. (5.64) as arguments of enthalpy ( ) A AW η are restricted to triangular regions, depending on the air fraction w A between 0 and 100% for selected pressures P as shown.At the upper entropy bound, the liquid phase is completely evaporated, given by the solution T(A, P) of case 2 in Appendix A.10, indicated as the "Dewpoint" lines in the figure.At the lower bound, the condensate of the wet-air sample starts freezing, Eq. (5.5), indicated here as the lines radiating from (η, A), and referred to as "Ice Formation" lines.The envelop on which the triangles' air-fraction maxima are located is Fig. 9. Valid entropy values of wet air computed from Eq. (5.64) as arguments of enthalpy h AW w A , η, P are restricted to triangular regions, depending on the air fraction w A between 0 and 100% for selected pressures P as shown.At the upper entropy bound, the liquid phase is completely evaporated, given by the solution T (A, P ) of case 2 in Appendix A10, indicated as the "Dewpoint" lines in the figure .At the lower bound, the condensate of the wet-air sample starts freezing, Eq. (5.64), indicated here as the lines radiating from (η, A), and referred to as "Ice Formation" lines.The envelop on which the triangles' air-fraction maxima are located is the triple line, shown dashed, where ice, liquid water and vapour coexist in the presence of air, as described in Sect.5.10, Eq. (S28.8).Note that the vapour-pressure lowering of water due to dissolved air is neglected in the equations.Freezing curves were computed with the library functions ice liq meltingtemperature si and liq air g entropy si, dewpoint curves using liq air dewpoint si and air g entropy si.
For running w A , the triple line is computed by calling the sequence set liq ice air eq at a, liq ice air temperature si, liq ice air pressure si and air g entropy si.
Potential enthalpy is a measure of the "heat content" of wet air in the sense of von Helmholtz' (1888) suggestion and was introduced into oceanography by McDougall (2003).The formula for the computation of the meteorological wet-bulb temperature from the enthalpy of humid air is given on page I.4-28 of WMO (2008).

Equilibrium humid air -ice
When humid air is in equilibrium with ice, its state is referred to as "saturated ice air" or the "frost point".The thermodynamic relations for this state are quite similar to those of the previous Sect.5.8 except that the Gibbs function g W of liquid water is replaced by the Gibbs function g Ih of ice.The condition for this equilibrium is equality of the chemical potentials of ice, Sect.2.2, and of water in humid air, Eq. (S12.15): (5.70) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. ( 5.70) is expressed by the system g AV = f AV A,T ,ρ AV + P /ρ AV (5.71) . (5.72) The independent variables in this scheme are the total pressure, P , the humid-air density, ρ AV , the temperature, T , and the air fraction, A. Expressing the chemical potential in Eq. ( 5.70) by means of Eq. ( 5.71) gives two equations in these four unknowns.Once two of the unknowns are specified, then the system is closed and may be solved numerically.Four important cases are discussed in detail in Appendix A11.
No matter which of the cases 1-4 considered in Appendix A11 is applied to compute the equilibrium between ice and humid air, the numerical solution of Eqs.(5.70)-(5.72)results in a consistent set of equilibrium values for A, T , P and ρ AV which is then available for the computation of any other property of either saturated humid air or ice in this state.
The related library functions are ice air rh wmo from a si, ice air rh cct from a si, ice air a from rh wmo si and ice air a from rh cct si.
With "ice air" we refer to a composite system of ice and humid air (e.g. a cirrus cloud) with the mass fractions w A of dry air, w V of vapour and w Ih of ice satisfying w A +w V +w Ih =1.The mutual equilibrium requires A sat (T ,P )=w A /w AV , with w AV =w A +w V =1−w Ih being the gaseous mass fraction.The Gibbs function of ice air reads (Feistel et al., 2010a), g AV A sat ,T ,P (5.73) and is a linear function of the air fraction, w A .Various ice-air properties are available from combinations of partial derivatives of the potential (Eq.5.73) with A sat (T ,P ) computed www.ocean-sci.net/6/633/2010/Ocean Sci., 6, 633-677, 2010 from Eqs. (5.70)-( 5.72) and case 3 from Appendix A11, see Table S25.For the computation of the partial T −P derivatives of g AI , the first derivatives of A sat (T ,P ) are required.
Taking the respective derivatives of Eq. ( 5.70) we get the isobaric drying rate, and the isothermal drying rate, of humid air, i.e. the decrease of its saturated air fraction A sat due to heating or compression.The chemical coefficient D A is defined in Eq. (S12.16).The latency operator AI of ice air used here is defined for the specific entropy, η AI =−g AI T , of the form (5.76) and for the specific volume v AI =g AI P of the form (5.77) The partial derivatives of the Gibbs function g AI w A ,T ,P , Eq. (5.73), of ice air are given in Table S24.Properties of ice air computed from this Gibbs function are given in Table S25.
For the description of isentropic processes such as the uplift of ice air in the atmosphere, enthalpy h AI w A ,η,P computed from the Gibbs function (Eq.5.73) is a useful thermodynamic potential: (5.78) For this purpose, temperature T in Eq. (5.78) must be determined from entropy η by numerically solving the equation (5.79) The partial derivatives of the enthalpy h AI w A ,η,P are computed from those of the Gibbs function, Table S24, as given in Table S26.
Selected properties of ice air computed from the enthalpy (Eq.5.78) and its partial derivatives are given in Table S27.
Many meteorological processes such as adiabatic uplift of an ice-air parcel conserve specific humidity and entropy to very good approximation.In particular, if a parcel is moved this way to some reference pressure P =P r , all of the thermodynamic properties given in Table S27 can be computed at that reference level from the partial derivatives of h AI w A ,η,P r .Such properties derived from the potential function h AI at the reference pressure are commonly referred to as "potential" properties in meteorology (von Bezold, 1888;von Helmholtz, 1888).Examples are the potential enthalpy, h θ , (5.80) the potential temperature, θ, in • C, obtained from Eq. (S27.2), and the potential density, ρ θ , from Eq. (S27.1), (5.82) The related library functions are ice air potenthalpy si, Eq. (5.80), ice air pottemp si, Eq. (5.81), and ice air potdensity si, Eq. (5.82).Ice air can exist only below an upper bound of entropy as shown in Fig. 10, given by either melting or the complete sublimation of the ice phase.

Equilibrium humid air -liquid water -ice
With the additional presence of air in the gas phase, the common triple point of water is expanded to a triple line in the A−T −P phase space, similar to the triple line of seawater, Fig. 3, in which the amount of salt present adds a new independent degree of freedom.When humid air, liquid water and ice coexist, the given conditions simultaneously satisfy the equilibrium conditions (Eqs.5.48 and 5.70) of equal chemical potentials of water in all three phases: (5.83) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.83) is expressed by the system g W (T ,P ) = f F T ,ρ W + P /ρ W (5.84) (5.82) The related library functions are ice_air_potenthalpy_si, Eq. (5.80), ice_air_pottemp_si, Eq. (5.81), and ice_air_potdensity_si, Eq. (5.82).Ice air can exist only below an upper bound of entropy as shown in Fig. 10, given by either melting or the complete sublimation of the ice phase.

T ri p le L in e T ri p le L in e T ri p le L in e T ri p le L in
A AI η are bounded above by roof-shaped curves, depending on the air fraction w A between 0 and 100% for selected pressures P as shown.At the entropy bound on the right, the ice phase is completely sublimated, given by the solution T = T subl (P vap ) of case 2 in Appendix A.11, and labelled "Frost Point" in the figure.At the left boundary lines radiating from the lower left portion of the figure, the ice phase starts melting, Eq. (5.5), labelled here as "Melting" lines.The locus of the roof tops at various pressures is the triple line, shown dashed, at which ice, liquid water and vapour coexist in the presence of air, as described in Sect.5.10, Eq. (S28.8).Freezing curves were computed with the library functions ice_liq_meltingtemperature_si and ice_air_g_entropy_si, and frost point curves were determined using ice_air_frostpoint_si and air_g_entropy_si.For running w A , the triple line is computed by calling the sequence set_liq_ice_air_eq_at_a, liq_ice_air_temperature_si, liq_ice_air_pressure_si and air_g_entropy_si.Fig. 10.Valid entropy values of ice air computed from Eq. (5.79) as arguments of enthalpy h AI w A , η, P are bounded above by roof-shaped curves, depending on the air fraction w A between 0 and 100% for selected pressures P as shown.At the entropy bound on the right, the ice phase is completely sublimated, given by the solution T =T subl (P vap ) of case 2 in Appendix A11, and labelled "Frost Point" in the figure.At the left boundary lines radiating from the lower left portion of the figure, the ice phase starts melting, Eq. (5.5), labelled here as "Melting" lines.The locus of the roof tops at various pressures is the triple line, shown dashed, at which ice, liquid water and vapour coexist in the presence of air, as described in Sect.5.10, Eq. (S28.8).Freezing curves were computed with the library functions ice liq meltingtemperature si and ice air g entropy si, and frost point curves were determined using ice air frostpoint si and air g entropy si.For running w A , the triple line is computed by calling the sequence set liq ice air eq at a, liq ice air temperature si, liq ice air pressure si and air g entropy si. (5.87) The independent variables in this scheme are the total pressure, P , the liquid density, ρ W , the humid-air density, ρ AV , the temperature, T , and the air fraction, A. Expressing the chemical potentials in Eq. (5.83) by means of Eqs.(5.84) and (5.86), gives four equations in these five unknowns so that one of the independent variables must be specified to complete the system.Once this is done, the remaining variables may be solved for iteratively as discussed in Appendix A12.Three important cases of different initially known properties corresponding to this system are discussed there.If the relative mass fractions of the three phases are required, then an additional condition is required to fix these quantities, since at constant T and P the water-ice mass ratio The independent variables in this scheme are the total pressure, P, the liquid density, ρ W , the humid-air density, ρ AV , the temperature, T, and the air fraction, A. Expressing the chemical potentials in Eq. (5.83) by means of Eqs.(5.84) and (5.86), gives four equations in these five unknowns so that one of the independent variables must be specified to complete the system.Once this is done, the remaining variables may be solved for iteratively as discussed in Appendix A.12. Three important cases of different initially known properties corresponding to this system are discussed there.If the relative mass fractions of the three phases are required, then an additional condition is required to fix these quantities, since at constant T and P the water-ice mass ratio can still change.The two additional cases, 4 and 5, considered in Appendix A.12 address this requirement.can still change.The two additional cases, 4 and 5, considered in Appendix A12 address this requirement.
Figure 11 corresponds to case 1 in the Appendix A12 with fixed dry air fraction, A. It illustrates that the temperature of wet ice air differs only very little from the triple-point temperature of water, almost independent of pressure, causing the adiabatic lapse rate under these conditions to be very small.Note that the curve shown here neglects the solubility of air in water which could result in temperature effects of similar order.
Figure 12 shows results corresponding to case 5 from Appendix A12 in which the dry-air fraction, w A , entropy, η, and the liquid fraction of the condensed part, w=w W /(w W +w Ih ) are specified.If an air parcel is lifted with the first two quantities fixed, then w varies between 0 at the melting level (completely frozen condensate), and 1 at the freezing level (completely molten condensate).Four valid wedge-shaped Wet-Ice-Air (WIA) regions are shown in this figure corresponding to pressures of 1000, 10 000, 101 325 and 1 000 000 Pa.Only points (w A , η) selected from these wedge-shaped regions permit valid solutions in this case.
Selected properties of wet ice air included as library routines are listed in Table S28.

Equilibrium humid air -seawater
Humid air in equilibrium with seawater, referred to as sea air, is subsaturated because the vapour pressure of seawater is lower than that of pure water.
In contrast to wet air, the liquid part of sea air can neither entirely evaporate nor freeze, i.e., as long as there is salt in the system there will always be a liquid fraction.Since there must be a gas fraction, too, whenever air is present, the composite system seawater -humid air can exist under ambient conditions only in two forms, with or without ice.

T ri p le L in e T ri p le L in e T ri p le L in e T ri p le L in e ICE AIR ICE AIR ICE AIR ICE AIR
F re e z in g F re e z in g F re e z in g F re e z in g M e lt in g M e lt in g M e lt in g M e lt in g

HUMID AIR HUMID AIR HUMID AIR HUMID AIR
WIA WIA WIA WIA Fig. 12 Valid entropy values of wet ice air, "WIA", computed from Eq. (S28.8) are restricted to narrow wedge-shaped regions depending on the air fraction w A between 0 and 100% for selected pressures P as shown.Only points (w A , η) selected from these regions permit valid solutions of case 5 discussed in Appendix A.12.At the upper entropy bound of wet ice air, wet air starts freezing, indicated as "Freezing" on the 1000 Pa case in the diagram.At the lower entropy bound of wet ice air, ice air starts melting, indicated as "Melting".The locus of the wedge tips at various pressures is the triple line, shown dashed, at which ice, liquid water and water vapour coexist in the presence of air, Eq. (S28.8).Freezing and melting curves were computed with the library functions ice_liq_meltingtemperature_si in conjunction with ice_air_g_entropy_si and liq_air_g_entropy_si.For running w A , the triple line is computed by calling the sequence set_liq_ice_air_eq_at_a, liq_ice_air_temperature_si, liq_ice_air_pressure_si and air_g_entropy_si.
Fig. 11 corresponds to case 1 in the Appendix A.12 with fixed dry air fraction, A. It illustrates that the temperature of wet ice air differs only very little from the triple-point temperature of water, almost independent of pressure, causing the adiabatic lapse rate under these conditions to be very small.Note that the curve shown here neglects the solubility of air in water which could result in temperature effects of similar order.
Fig. 12 shows results corresponding to case 5 from Appendix A.12 in which the dry-air fraction, w A , entropy, η, and the liquid fraction of the condensed part, w = w W / (w W + w Ih ) are specified.If an air parcel is lifted with the first two quantities fixed, then w varies between 0 at the melting level (completely frozen condensate), and 1 at the freezing level (completely molten condensate).Four valid wedge-shaped Wet-Ice-Air (WIA) regions are shown in this Fig.12. Valid entropy values of wet ice air, "WIA", computed from Eq. (S28.8) are restricted to narrow wedge-shaped regions depending on the air fraction w A between 0 and 100% for selected pressures P as shown.Only points (w A , η) selected from these regions permit valid solutions of case 5 discussed in Appendix A12.At the upper entropy bound of wet ice air, wet air starts freezing, indicated as "Freezing" on the 1000 Pa case in the diagram.At the lower entropy bound of wet ice air, ice air starts melting, indicated as "Melting".The locus of the wedge tips at various pressures is the triple line, shown dashed, at which ice, liquid water and water vapour coexist in the presence of air, Eq. (S28.8).Freezing and melting curves were computed with the library functions ice liq meltingtemperature si in conjunction with ice air g entropy si and liq air g entropy si.For running w A , the triple line is computed by calling the sequence set liq ice air eq at a, liq ice air temperature si, liq ice air pressure si and air g entropy si.
The first case is considered in Sect.5.12.Note that sea air does not contain ice at temperatures above the freezing point of seawater.Nonetheless, air saturation and relative humidity of humid air is defined relative to ice if the temperature is below the freezing point of pure water, even though no stable ice phase is present in the interval between the freezing temperatures of pure water and of the system's seawater component.Similar to Eq. (5.26), the condition for this equilibrium is equal chemical potentials of water in seawater, Eq. (S7.12), and of water in humid air, Eq. (S12.15): ∂g SW ∂S A T ,P . (5.88) In terms of the Primary Standard functions and their independent variables (Sect.2), Eq. (5.88) is expressed by the system g SW (S A ,T ,P )=f F T ,ρ W +P /ρ W +g S (S A ,T ,P ) (5.89) g AV = f AV A,T ,ρ AV + P /ρ AV (5.91) . (5.92) The independent variables in this scheme are the total pressure, P , the pure-water density, ρ W , the humid-air density, ρ AV , the temperature, T , the Absolute Salinity, S A , and the air fraction, A. Note that ρ W is merely a formal property here -the density that liquid pure water has at given T and P .Expressing the chemical potentials in Eq. ( 5.88) by means of Eqs.(5.89) and (5.91), provides three equations in the six unknowns so three of the independent variables must be specified to complete the system.Once this is done, the remaining unknowns may be determined by iterative numerical methods.Two important cases are considered in detail in Appendix A13.
Selected properties of sea air are given in Table S29.The latent heat L SA P of sea air is defined here as the enthalpy required to evaporate a small amount of water from seawater to humid air by heating at constant pressure.A derivation of the latent-heat equation is given in (Feistel et al., 2010a), similar to the latent heat of melting sea ice, Eq. (5.23).

Equilibrium humid air -seawater -ice
In contrast to sea air, Sect.5.11, humid air in equilibrium with sea ice, referred to as sea-ice air here, is saturated because it is in equilibrium with salt-free ice, Sect.5.9.The phases of sea-ice air are simultaneously in pairwise mutual equilibria, seawater with ice (sea ice, Sect.5.4), ice with humid air (ice air, Sect.5.9), and seawater with humid air (sea air, Sect.5.11).Most of the properties of sea-ice air are available from the related library functions described in those sections, therefore we have refrained from implementing a special sea-ice-air module.For completeness, we mention that the equilibrium conditions for sea-ice air consist of two equations between the chemical potentials of water in the three present phases, Eqs.(5.11), (5.70) and (5.88): The latent heat of sea-ice air includes the transfer of water between the phases by melting, evaporation and sublimation.
The resulting isobaric latent heat of sea-ice air is, expressed per kg of molten ice (Feistel et al., 2010a), =w S /S A are the gaseous and the liquid fractions, and w A and w S are the given constant mass fractions of air and of salt in the sea-ice-air sample.

Summary and short discussion
The mutually consistent formulations of thermodynamic potentials for liquid water, water vapour, ice, seawater and humid air are now available and permit the numerical computation of a wealth of thermodynamic properties of the geophysical fluids, their mixtures, composites and phase transitions.The new seawater standard TEOS-10 (IOC et al., 2010) together with its collection of background papers developed by WG127 in cooperation with IAPWS is based on this physically and mathematically rigorous building-block concept (Feistel et al., 2008).To support the practical use and general implementation of TEOS-10, WG127 has developed a source code library that provides easy access to a large selection of properties and may serve as a guide for writing customized application code using the new standard.
The library is hierachically organized; all available properties are computed exclusively from the Primary Standard, i.e., level 1 of the code, by merely mathematical and numerical methods.The concept of the Primary Standard is intentionally similar to axiomatic systems in mathematics which possess the general properties of consistency, independence and completeness.These properties ensure that the Primary Standard contains all necessary but no redundant components, and prevents the computation of contradicting results.The higher levels obey the conditions of a mathematical semi-order structure; code of a given level does not refer to code of higher levels, thus avoiding direct or indirect recursion.
In the case of seawater, it would be most natural to provide access to only the saline component of the Gibbs function (Eq.2.2) at level 1 and not permit access to the individual coefficients (Eqs.2.3-2.5) of the salinity expansion.However, it is necessary to have access to the individual temperature and pressure dependent coefficients in order to rigorously consider numerical limits as S A tends to zero.Thus, these fundamental building blocks are made individually available at level 1.To obey the independence rule for level 1 routines, it is then necessary to place the Gibbs function (Eq.2.2) at level 2, which is not subject to this condition.A similar situation appears in the case of humid air.The Primary Standard provides the Helmholtz function of dry air (Eq.2.6) together with the air-water virial coefficients as the fundamental information from which the properties of humid air can be computed.To ensure independence for level 1 routines, the Helmholtz function of humid air, Eq. (2.7), and the cross-over Helmholtz function (Eq.2.13) are then implemented in level 2 of the library.Note that while the library is constructed to strictly adhere to the development based on axiomatic results at level 1, we have discussed the potentials of seawater and humid air together with the level-1 functions in Sect. 2 of this paper because of their close logical relations.
In addition to the Primary Standard, the library provides easy access to other thermodynamic potential functions derived from the Primary Standard.Available are Helmholtz functions that are computed from temperature and density, Gibbs functions computed from temperature and pressure, enthalpy functions computed from entropy and pressure, and implicitly entropy as a potential computed from enthalpy and pressure.A list of explicitly implemented potential functions is given in Table 1.From each of these potential functions, all thermodynamic properties of the particular system can be computed; the library provides an extensive but still selective set of relevant properties.For additional composite systems such as seawater with humid air, several properties are available from the library even though related potential functions were not implemented explicitly.
Further details on organization, content and access to the library are contained in the companion paper (Wright et al., 2010a).
We begin with a first guess for the density and linearize Eq. (3.3) with respect to changes of density under the assumption that our first guess is sufficiently near the desired root that the linearization is valid.This gives which, by Newton iteration, permits the computation of a density improvement ρ from a given estimate ρ at fixed values of T and P .The iteration will converge to a single fluid value for supercritical conditions, and to one or the other of the distinct vapour and liquid density values, depending on the initial density estimate, for subcritical conditions.Once the solution ρ is known, all thermodynamic properties of fluid, liquid or vapour, can be computed by either of two formally different methods: (i) For the direct access to liquid water or vapour properties, the required function of T and ρ is called from Sect.3.1, Table S2.
(ii) For the indirect use of water properties as a part of e.g.seawater properties, the Gibbs function, Eqs.(A1) or (A2) and its derivatives must be made available to the related calling functions, see Eq. (2.1), Sect.4.2 and Table S7.
To determine liquid or vapour solutions of Eq. (A1) where either one or both of these may exist along with possible spurious numerical solutions, the choice of an initial starting point must be made carefully to lie inside the "convergence radius" of the desired attractor.It is therefore useful to consider liquid and vapour separately in each of the subcritical range, the critical region and the supercritical range, as shown in Fig. A1.Liquid and vapour can be distinguished from each other by their different densities and entropies in the vicinity of the saturation line which is the curve connecting the triple point (TP) with the critical point (CP) in Fig. A1.On the saturation line, both phases can coexist in physical space, separated by an interface (the "water surface") across which the properties change abruptly.The saturation line is defined by equal chemical potentials, i.e. equal specific Gibbs energies (Eqs.A1, A2) of the two phases.Except for this mutual equality, there is no particular distinguishing property within either phase which might separate the saturation state from its surrounding T −P states (Landau and Lifschitz, 1964).In the vicinity of the saturation line, the phase with the lower Gibbs energy (Eqs.A1, A2) is stable, the other state exists, but is metastable.Here, metastable means stable with respect to infinitesimal fluctuations but unstable with respect to certain macroscopic perturbations, namely the emergence of finite volumes of the coexisting phase (nucleation of supercritical bubbles or droplets).At a greater distance from the saturation line, the state with higher Gibbs energy may be unstable if thermodynamic stability criteria such as positive compressibility are violated.(Independent of any simultaneous existence of other phases, a negative compressibility would amplify any pressure-density fluctuation, causing the fluid to collapse.)The boundary between metastable and unstable existence is regarded as the spinodal line beyond which the phase can no longer stably and homogeneously exist.
Since it is practically impossible to measure thermodynamic properties on or beyond the spinodal, its location in the phase diagram is not exactly known from empirical equations of state.Even though the IAPWS-95 formulation extrapolates well into the metastable regions (Wagner and Pruß, 2002;Feistel et al., 2008), with increasing distance from the saturation line the values computed from the Helmholtz function outside its validity range will become unreliable.This turns out not to be a major issue for determination of the seawater Gibbs function since we require only relatively small excursions into the metastable regions to deal with the effects of shifted phase transition boundaries in the presence of sea salt.
Panels a and b of Fig. A1 indicate the initializations used in the library for our iterative solutions for the liquid and vapour phases, respectively.Figure A2 provides additional information regarding the industrial formulation IF-97 (IAPWS, 2007;Wagner and Kretzschmar, 2008) referred to in Fig. A1.
Starting in either fluid state at a point near CP located along the saturation line TP and CP, we may circumscribe the critical point along a closed T −P path.Along any such curve, the properties change only gradually; nothing like a transition point between liquid and vapour is encountered.Since, for numerical purposes, we distinguish between the Gibbs functions of liquid, Eq. (A1), and vapour, Eq. (A2), we need to specify such a transition point for technical rather than for physical reasons.Here we define the Gibbs functions Eqs.(A1) and (A2) to be different at subcritical conditions (T <T C and P <P C ) and to be identical at supercritical conditions (T ≥T C or P ≥P C ).The critical temperature of water is T C =647.096 K and the critical density is ρ C =322 kg m −3 (IAPWS, 2009a); the critical pressure follows from Eq. (4.1) to be P C =22.064 MPa.In order to cover metastable states of liquid water as required in the regions of vapour-pressure lowering or freezing-point lowering caused by the presence of dissolved sea salt, the Gibbs function for liquid water is also available for T and P in regions extending somewhat beyond the saturation curve and beyond the melting curve.
According to our numerical definition of liquid, vapour and fluid states, the initial values required for the iteration of Eq. (A3) can be chosen identically for the fluid density in the supercritical region (T ≥T C or P ≥P C ), as shown in Fig. A1, but must be different for liquid and vapour in the subcritical quarter (T <T C and P <P C ).In the subcritical region, separate Gibbs functions are available from the industrial formulation IF-97 (IAPWS, 2007;Wagner and Kretzschmar, 2008) have been determined by regression to IAPWS-95 data points in the stable liquid, vapour and fluid region with an r.m.s.deviation of 1% or less for each phase; the resulting coefficients are given in Table A1.The cubic polynomial used for Eq.(A4) permits analytical inversion to determine ( ) in regions 1 and 2 (Fig. A2) as defined therein which provide excellent starting values for the liquid and the vapour state.These Gibbs functions can also be used for the fluid, region 1 below 623.15K and 100 MPa, and region 2 between 273.15 K and 1073.15K, and below 16.529 MPa.In the sublimation region and in the supercritical region, the ideal-gas density, ρ=P /(R W T ), provides a sufficient starting estimate below 273.15K and above 650 K, and a constant value of ρ=1000 kg m −3 can be used below 650 K, Fig. A1.These latter choices are sufficient to ensure numerical convergence but do not necessarily optimize the speed of the code.Additional considerations apply to the immediate neighbourhood of the critical point as discussed below.Note that all of these rules are built into the library routines discussed in Part 2 (Wright et al., 2010a) so that the user can make use of the routines without dealing with (or even being fully aware of) the details.
The critical region is defined here as the T −P rectangle 623.15-650K and 16.529-35 MPa, Fig. A1.The coefficients of an auxiliary cubic polynomial equation of state have been determined by regression to IAPWS-95 data points in the stable liquid, vapour and fluid region with an r.m.s.deviation of 1% or less for each phase; the resulting coefficients are given in Table A1.The cubic polynomial used for Eq.(A4) permits analytical inversion to determine ρ (T ,p) for both the liquid and the vapour branches.The critical point of Eq. (A4) was chosen to be identical with the The initial densities for the iteration, Eq. (A3), in the critical region are computed from the intersection points of the horizontal isobars with the isotherms as shown in Fig. A3.In the subcritical range, T < TC and P < PC, there exist three solutions, the vapour density to the left of the isotherm maximum, the liquid density to the right of the minimum, and an extraneous unstable solution in between the extrema.The curve (not shown) connecting the minima and the maxima of adjacent isotherms, which passes smoothly through the critical point, is the spinodal of the auxiliary equation.Beneath the spinodal, the compressibility is negative, , thermodynamic stability is violated and no stable single-phase states can exist.By means of this stability gap, the spinodal separates low-density vapour from high-density liquid on the particular isotherm.At the critical point, maximum, minimum and inflection point coincide, and at supercritical temperatures only one fluid solution exists for any given Fig. A3.Selected isotherms of the auxiliary cubic equation of state, Eq. (A4), in the critical region 623.15<T<650 K, 16.529<P <35 MPa.Shown in bold is the saturation-pressure curve of IAPWS-95, separating the single-phase region above from the two-phase region below.The critical point is at T C =647.096 K, ρ C =322 kg m −3 , P C =22.064 MPa.Given a line of constant subcritical pressure P and an isotherm T , their intersection points with positive slopes provide either the density estimate for ρ V (T , P ) of vapour (on the left branch) or for ρ V (T , P ) of the liquid (on the right), separated from each other by the unstable region of negative slopes.
IAPWS-95 critical point (Fig. A3) through the specifications of T C , P C , and ρ C .
The initial densities for the iteration, Eq. (A3), in the critical region are computed from the intersection points of the horizontal isobars with the isotherms as shown in Fig. A3.In the subcritical range, T <T C and P <P C , there exist three solutions, the vapour density to the left of the isotherm maximum, the liquid density to the right of the minimum, and an extraneous unstable solution in between the extrema.The curve (not shown) connecting the minima and the maxima of adjacent isotherms, which passes smoothly through the critical point, is the spinodal of the auxiliary equation.Beneath the spinodal, the compressibility is negative, (∂ρ/∂P ) T <0, thermodynamic stability is violated and no stable single-phase states can exist.By means of this stability gap, the spinodal separates low-density vapour from high-density liquid on the particular isotherm.At the critical point, maximum, minimum and inflection point coincide, and at supercritical temperatures only one fluid solution exists for any given pressure.Below the critical temperature, a single solution from the liquid branch is computed for P >P C which is considered a supercritical fluid state according to our numerical definition of the liquid and vapour functions (Eqs.A1 and A2).Very close to the critical point, initial densities computed from the auxiliary cubic equation of state may falsely be located inside the spinodal of IAPWS-95 and thus prevent convergent iteration.In this highly specialized case, applications may need better starting values than those from the cubic polynomial, e.g.find exact densities at the gas and liquid spinodal points from the condition (∂P /∂ρ) T =0 and use one of them to confine ρ (T ,P ) for a bisection iteration method such as the secant or Brent algorithms.Details of the universal critical properties are available from Stanley (1971), Anisimov (1991), Kurzeja et al. (1999), Skripov andFaizullin (2006), or Ivanov (2008).
A2 Seawater temperature from salinity, entropy and pressure (Sect.4.3) To compute the specific enthalpy potential and its partial derivatives from the Gibbs function g SW (S A ,T ,P ) of seawater, the value of T appearing in the expression for the enthalpy, Eq. (4.5), must be obtained from knowledge of the entropy, η, along with the salinity and pressure values.The required temperature is obtained by numerically solving the Eq.(4.6) to give T =T (S A ,η,P ).
To solve Eq. (4.6) we first linearize η=−g SW T with respect to small changes of temperature to obtain the equation which can be used to iteratively update the value of T at given values of S A , η, P .Since the heat capacity of water is rather constant under different oceanic conditions and Eq. ( A5) has an unambiguous solution in the region of oceanographic interest, the simple linear estimate provides a sufficiently accurate initial temperature to ensure convergent iteration of Eq. (A5).

A3 Saturated water vapour conditions (Sect. 5.1)
To numerically determine the conditions corresponding to the saturation point (frequently referred to as the boiling point or dewpoint) of pure water, we first linearize the three Eqs.(5.2)-(5.4)with respect to small changes of the four unknowns T , P , ρ V and ρ W , which gives: For brevity, f F T ,ρ W is abbreviated here by f W , and similarly for f V as well as their partial derivatives.To obtain Eq. (A7), Eq. ( 5.2) was first expanded and then simplified by using Eqs.(A8) and (A9).When the equilibrium point is reached, Eq. (A7) takes the form of the Clausius-Clapeyron equation as its right-hand side vanishes.To solve the system (Eqs.A7-A9) for T , P , ρ V and ρ W , a fourth equation must be added which specifies an additionally imposed condition, usually one of T =0 (for specified temperature) or P =0 (for specified pressure).
Auxiliary empirical equations are used to determine initial estimates for T , P , ρ V and ρ W .
For 550 K<T <T C , an estimate of the vapour density ρ V on the saturation curve is available from the cubic correlation polynomial To iteratively determine freezing and melting conditions for pure water, we first linearize the two Eqs.(5.6), (5.7) with respect to small changes of the three unknowns T , P and ρ W to obtain: The function f F T ,ρ W is abbreviated here by f W , and similarly for its partial derivatives.To iteratively solve the system (Eqs.5.9, 5.10) for T , P and ρ W using Eqs.(A15), (A16), a third equation must be added which specifies an additionally imposed condition, usually T =0 (if the temperature is specified) or P =0 (if the pressure is specified).Auxiliary empirical equations are used to determine initial estimates for T , P and ρ W . www.ocean-sci.net/6/633/2010/Ocean Sci., 6, 633-677, 2010 An estimate of the freezing temperature T in K as a function of the absolute pressure P in Pa is obtained from a correlation fit between 252 and 273 K: An estimate of the density ρ W in kg m −3 of the freezing liquid as a function of absolute temperature T in K is obtained from a correlation fit between 252 and 273 K: with rms error equal to 1.2 E-4 in ρ W /ρ W t .The constants are ρ W t =999.792520 031 621 kg m −3 , a 1 =−1.785829 814 921 13, a 2 =−12.232508 430 673 4, a 3 =−52.823693 643 352 9.

A5 Equilibrium conditions for ice and water vapoursublimation (Sect. 5.3)
To determine conditions for sublimation, we first linearize the two Eqs.(5.9, 5.10) with respect to small changes of the three unknowns T , P and ρ V to obtain: The function f F T ,ρ V is abbreviated here by f V , and similarly for its partial derivatives.To iteratively solve the system (Eqs.5.9, 5.10) for T , P and ρ V using Eqs.(A19), (A20), a third equation must be added which specifies an additionally imposed condition, usually T =0 (if the temperature is specified) or P =0 (if the pressure is specified).Auxiliary empirical equations are used to determine initial estimates for T , P and ρ V .
The sublimation temperature T in K as a function of the absolute pressure P in Pa is estimated using the Clausius-Clapeyron equation (for details see Feistel and Wagner, 2007): The constants are P t =611.654771 007 894 Pa, T t =273.16K, the sublimation heat, h=2 834 359.445 433 54 J kg −1 and the specific gas constant of water, R W =461.518 05 J kg −1 K −1 .The density ρ V in kg m −3 of the condensing vapour is estimated as a function of absolute temperature T in K from the Clausius-Clapeyron, Eq. (A21), in combination with the ideal-gas equation: The constants are the same as for Eq.(A21).
A6 Equilibrium conditions for ice in seawater (Sect.5.4) To determine the conditions under which ice exists in equilibrium with seawater, we first linearize the two Eqs.(5.12), (5.13) with respect to small changes of the four unknowns S A , T , P and ρ W to obtain: To iteratively solve the system (Eqs.5.12, 5.13) for S A , T , P and ρ W using Eqs.(A23), (A24), two further equations must be added which specify an additionally imposed pair of conditions, commonly taken to be T =0 and S A =0 (if the temperature and the salinity are specified) or P =0 and S A =0 (if the pressure and the salinity are specified).Auxiliary empirical equations are used to determine initial estimates forS A , T , P and ρ W .
In the oceanographic range, the pure-water part Eq. (A25) of the International Equation of State of Seawater EOS-80 provides a very good estimate of the density ρ W =1/v W as a function of temperature and pressure (Millard, 1987):  A2.Equation (A18) provides an estimate of the density ρ W in kg m −3 of freezing pure water as a function of absolute temperature T in K from a correlation fit between 252 and 273 K.
Equation (A26) provides an estimate of the brine salinity S A in kg kg −1 of sea ice at given absolute temperature T in K and absolute pressure P in Pa from the empirical laws of Clausius-Clapeyron and Raoult: The IAPWS-95 triple point is P t =611.654771 007 894 Pa, T t =273.16K, the Raoult coefficient is α=−0.217(Feistel et al., 2008) and the Clausius-Clapeyron coefficient is χ =−74.3×10−9 K Pa −1 (Feistel and Wagner, 2006).When solved for the freezing temperature T (S A ,P ), Eq. (A26) gives and when solved for the melting pressure P (S A ,T ), it gives A7 Conditions for seawater in equilibrium with water vapour (Sect.5.5) To determined conditions under which water vapour will be in equilibrium with seawater, we first linearize the three Eqs.(5.27)-(5.29)with respect to small changes of the five unknowns S A , T , P , ρ V and ρ W to obtain: To iteratively solve the system (Eqs.5.27, 5.28, 5.29) for S A , T , P , ρ V and ρ W using Eqs.(A29)-(A31), two further equations must be added which specify an additionally imposed pair of conditions, commonly T =0 and S A =0 (if the temperature and the salinity are specified) or P =0 and S A =0 (if the pressure and the salinity are specified).
Auxiliary empirical equations are used to determine initial estimates for S A , T , P , ρ V and ρ W .
The function (Eq.A32) is the inverse of Eq. (A10) and estimates the boiling temperature T in K of the seawatervapour equilibrium at given brine salinity S A in kg kg −1 and absolute pressure P in Pa from the Clausius-Clapeyron and Raoult laws for T <640 K:  (Feistel, 2008).For the estimates required here, the raw conversion T 48 ≈T and S P ≈1000 S A is sufficiently precise.Formula (Eq.A34), obtained from Eq. (A33) and Raoult's law, computes a brine salinity estimate S A in kg kg −1 for seawater-vapour equilibrium at given absolute temperature T in K and absolute pressure P in Pa: and expand the resulting three Eqs.(5.48), (5.50) and (5.52) with respect to small changes of the five variables: For brevity, f F T ,ρ W is abbreviated here by f W as well as its partial derivatives.For the numerical solution, two additional conditions are required, such as specifying the temperature and pressure, so that T =0 and P =0.Appropriate starting values of the remaining unknowns must be specified for their iterative determination.Four important special cases are considered in the following.
Case 1: Equilibrium at given air fraction, A, and temperature, T At given A and T , humid air can approximately be considered as an ideal mixture of air and vapour.The partial pressure P vap of vapour is computed from the vapour pressure of liquid water at given T by solving Eq. (5.1).The vapour density follows from Eq. (4.3) as ρ V =1/g V P (T ,P vap ).For the dry-air density we have ρ A =ρ V ×A/(1−A).The partial pressure of dry air is computed from Eq. (S5.11) as P air = ρ A 2 f AV ρ 1,T ,ρ A .This provides an estimate for the total pressure, P =P vap +P air .With A, T and P available, the required density estimate of liquid water, ρ W =1/g W P (T ,P ), and of humid air, ρ AV =1/g AV P (A,T ,P ), are easily calculated from the Gibbs functions, Eqs.(4.2) and (4.37).Using A=0 and T =0, the linear system (Eqs.A45-A47) can now be solved iteratively for P , ρ W and ρ AV .
In particular, this solution provides the pressure P (A,T ) of saturated humid air as a function of the air fraction and the temperature.
The equilibrium is computed this way by the library call set liq air eq at a t or by the function liq air condensationpressure si.
Case 2: Equilibrium at given air fraction, A, and pressure, P At given A and P , humid air can approximately be considered as an ideal mixture of air and vapour.The partial pressure P vap =x AV V P of vapour is computed from the total pressure P and the mole fraction x AV V (A), Eq. (2.11).In turn, the boiling temperature T =T boil (P vap ) of water is computed from Eq. (5.1).With A, T and P available, the required density estimate of liquid water, ρ W =1/g W P (T ,P ), and of humid air, ρ AV = 1/g AV P (A,T ,P ), are easily calculated from the related Gibbs functions, Eqs.(4.2) and (4.37).Using A=0 and P =0, the linear system (Eqs.A45-A47) can now be solved iteratively for T , ρ W and ρ AV .
In particular, this solution provides the dewpoint temperature T (A,P ) of humid air as a function of the air fraction and the pressure.
This approach is used to compute the equilibrium with the library call set liq air eq at a p or using the function liq air dewpoint si.
Case 3: Equilibrium at given temperature, T , and pressure, P At given T and P , humid air can approximately be considered as an ideal mixture of air and vapour.The partial pressure P vap of vapour is computed from the vapour pressure of liquid water at given T from solving Eq. (5.1).The vapour density follows from Eq. (4.3) as ρ V =1/g V P (T ,P vap ) and the air density from ρ A =1/g AV P (1,T ,P −P vap ).Now the air fraction is available from A=ρ A / ρ A +ρ V .With A, T and P available, the required density estimate of liquid water, ρ W =1/g W P (T ,P ), and of humid air, ρ AV =1/g AV P (A,T ,P ), are easily calculated from the related Gibbs functions, Eqs.(4.2) and (4.37).Using T =0 and P =0, the linear system (Eqs.A45-A47) can now be solved iteratively for A, ρ W and ρ AV .
In particular, this solution provides the specific humidity q=1 − A(T ,P ) of saturated humid air as a function of the temperature and the pressure.
The equilibrium is computed using this approach with the library call set liq air eq at t p or using the function liq air massfraction air si.

Fig. 1 .
Fig. 1.Panel (a) Validity region (bounded by bold lines) of the IAPWS-95 Helmholtz potential for fluid water with isobars as indicated.Panel (b) Magnified view of the small region corresponding to the standard oceanographic ("Neptunian") range.TP: triple point gas-liquid-solid, CP: critical point.The deviation of the vapourpressure line from the 101 325 Pa isobar in the liquid region is below the graphical resolution of panel (b).Freezing-point lowering occurs with the addition of sea salt.To deal with this effect in the case of seawater, the extension of the pure water properties into the metastable liquid region just above the line marked "Freezing Point Lowering" is required.

Fig. 2 :
Fig. 2: Range of validity (bold curves) of the Gibbs function of ice Ih and uncertainty of density.

Fig. 2 .
Fig. 2. Range of validity (bold curves) of the Gibbs function of ice Ih and uncertainty of density.

Fig. 3
Fig.3Range of validity of the IAPWS-08 Gibbs function of seawater and uncertainty of density estimates calculated from this function.Region A: oceanographic standard range, B: extension to higher salinities, C: hot concentrates, D: zero-salinity limit, E: extrapolation into the metastable region below 0 °C.

Fig. 3 .
Fig. 3. Range of validity of the IAPWS-08 Gibbs function of seawater and uncertainty of density estimates calculated from this function.Region A: oceanographic standard range, B: extension to higher salinities, C: hot concentrates, D: zero-salinity limit, E: extrapolation into the metastable region below 0 • C.

Fig. 4 .
Fig. 4. Validity range of the Helmholtz function for humid air, Eq.(2.7).For oceanographic and meteorological applications it is unnecessary to consider liquid or solid air.Thus, we restrict consideration of Eq. (4.37) as follows: (i) for temperatures above the critical temperature of dry air, T >T c =132.5306 K, all density values occurring between the pressure bounds are permitted; and (ii) for subcritical temperatures T <T c , only densities below the dewpoint curve of dry air, indicated by "Dewpoint" are permitted.The resulting validity boundary for dry air is shown in bold."CP" is the critical point of dry air.To consider humid air, virial coefficients are required.The validity range in temperature of the third virial coefficients is shown by horizontal lines.Additionally, the pressure on saturated humid air is restricted to 5 MPa(Hyland and Wexler 1983), not shown.
Fig. 5.Saturation curves A sat (T ,P ) of humid air at the pressures 101.325, 50, 20 and 10 kPa, as indicated.Panel (a) shows results for temperatures above the freezing point, computed by solving Eq. (5.48) using the library function liq air massfraction air si, Eq. (S21.9), and panel (b) shows results for temperatures below the freezing point, computed by solving Eq. (5.70) using the function ice air massfraction air si, Eq. (S25.10).Valid air fraction values A are located above the particular saturation curve, A≥A sat (T ,P ), in the region indicated by "HU-MID AIR".In the presence of ice-free seawater, the validity range for A is more restricted, A≥A cond (S A ,T ,P )≥A sat (T ,P ), by the condensation value A cond , computed from the function sea air massfraction air si, Eq. (S29.1).
from knowledge of salinity, entropy and pressure.Given values of S A , η and P , the corresponding value of T is obtained by numerically solving the equation

Fig. 7 :
Fig.7: Temperature-pressure phase diagram of seawater in the vicinity of the triple point.At different salinities, the triple point (TP), i.e. the equilibrium between liquid seawater, ice and vapour is displaced along the sublimation line (in bold) of the ice-vapour equilibrium.Note that the triple-point pressure can change by a factor of 2 while the vapour-pressure lowering at constant temperature is only of order 10% or less.
Pressure Data of Saturated Humid Air, magnified

Fig. 11 :
Fig.11: Temperature of wet ice air as a function of the air fraction, T(A), computed as described under case 1, Appendix A.12.

Fig. 11 .
Fig. 11.Temperature of wet ice air as a function of the air fraction, T (A), computed as described under case 1, Appendix A12.
Fig. A1.Initial estimates used in the library for the numerical computation of density from pressure by iteratively solving Eq. (A3), either for liquid water, panel (a), or water vapour, panel (b).The region surrounding the critical point (CP) is treated separately, as shown in Fig.A3.Here we make use of the Gibbs functions "g(p, T )" in two of the five regions defined in IF-97, region 1 (liquid/fluid) and region 2 (vapour/fluid) as shown in Fig.A2.Panels (a) and (b) differ only for subcritical conditions T <T c and P <P c ; otherwise there is only one solution corresponding to the unique Gibbs function for fluid water.TP is the ice-liquid-vapour triple point.The saturation curve connects TP with CP and separates liquid above from vapour below.To the far left, ice Ih (ICE) is separated from the liquid by the melting curve, above TP, and from the vapour by the sublimation curve, below TP.Note that the melting curve above 200 MPa belongs to forms of ice other than Ih, the ambient hexagonal phase.For vapour below 273.15 K, the ideal-gas equation is used.

60Fig.
Fig. A2: Thermodynamic relations available from the Industrial Formulation IF-97 (IAPWS, 2007; Wagner and Kretzschmar, 2008) defined in different temperature-pressure regions, derived with reduced accuracy from the IAPWS-95 Helmholtz function for different independent variables.Here we make use of the Gibbs functions "g(p, T)" available in region 1 (liquid/fluid) and in region 2 (vapour/fluid) as shown in Fig. A1.In region 3, separate equations for the specific volume are available for various subregions which are not used here.Region 4 is the saturation curve.Graphics reproduced from IAPWS (2007), with permission of IAPWS.
Fig. A2.Thermodynamic relations available from the Industrial Formulation IF-97 (IAPWS, 2007; Wagner and Kretzschmar, 2008) defined in different temperature-pressure regions, derived with reduced accuracy from the IAPWS-95 Helmholtz function for different independent variables.Here we make use of the Gibbs functions "g(p, T )" available in region 1 (liquid/fluid) and in region 2 (vapour/fluid) as shown in Fig. A1.In region 3, separate equations for the specific volume are available for various subregions which are not used here.Region 4 is the saturation curve.Graphics reproduced from IAPWS (2007), with permission of IAPWS.

Table 1 .
Hierarchy of thermodynamic potentials and their elementary constituents implemented in the library.

Table A1 :
Coefficients of the auxiliary critical equation of state, Eq. (A4).The critical region is defined here as the T-P rectangle 623.15 -650 K and 16.529 -35 MPa, Fig. A1.The coefficients of an auxiliary cubic polynomial equation of state chosen to be identical with the IAPWS-95 critical point (Fig.A3) through the specifications of Tc, Pc, and C Selected isotherms of the auxiliary cubic equation of state, Eq. (A4), in the critical region 623.15 < T < 650 K, 16.529 < P < 35 MPa.Shown in bold is the saturation-pressure curve of IAPWS-95, separating the single-phase region above from the two-phase region below.The critical point is at TC = 647.096K, ρC = 322 kg m -3 , PC = 22.064 MPa.Given a line of constant subcritical pressure P and an isotherm T, their intersection points with positive slopes provide either the density estimate for ρ V (T, P) of vapour (on the left branch) or for ρ W (T, P) of the liquid (on the right), separated from each other by the unstable region of negative slopes.

Table A1 .
Coefficients of the auxiliary critical equation of state, Eq. (A4).
a 1