Chemical speciation models based upon the pitzer activity coefficient equations, including the propagation of uncertainties: Artificial seawater from 0 to 45 °C

Accurate chemical speciation models of solutions containing the ions of seawater have applications in the calculation of carbonate system equilibria and trace metal speciation in natural waters, and the determination of pH. Existing models, based on the Pitzer formalism for the calculation of activity coefficients, do not yet agree with key experimental data (potentiometric determinations of H + and Cl activity products in acidified artificial seawaters) and, critically, do not include uncertainty estimates. This hampers applications of the models, and also their further development (for which the uncertainty contributions of individual ion interactions and equilibrium constants need to be known). We have therefore implemented the models of Waters and Millero (Mar. Chem. 149, 8-22, 2013) and Clegg and Whitfield (Geochim. et Cosmochim. Acta 59, 2403-2421, 1995) for artificial seawater, within a generalised treatment of uncertainties, as a first step towards a more complete model of standard seawater and pH buffers. This addition to the model enables both the total uncertainty of any model-calculated quantity (e.g., pH, speciation) to be estimated, and also the contributions of all interaction parameters and equilibrium constants. Both models have been fully documented (and some corrections made). Estimates of the variances and covariances of the interaction parameters were obtained by Monte Carlo simulation, with simplifying assumptions. The models were tested Jo ur na l P re -p ro of Journal Pre-proof

These differences compare with an overall uncertainty in the measured EMFs of about 0.04 mV.
Total uncertainties for calculated EMFs of the solutions were dominated by just a few contributions: mainly H + -Cl -, Na + -Cl -, and H + -Na + -Clionic interactions, and the thermodynamic dissociation constant of HSO 4 -

Introduction
Artificial seawater contains the major ions Na + , Mg 2+ , Ca 2+ , K + , Cland SO 4 2-(with H + , HSO 4 -, OHand MgOH + included as minor species for the purposes of acid-base calculations), and chemical speciation models using the Pitzer equations for activity coefficients have been developed for solutions containing these ions over a range of temperatures (Campbell et al., 1993;Clegg and Whitfield, 1995;. Such models form the main element of larger models of standard seawater (which also contains additional solutes including fluoride, and various borate and carbonate species), and can also be extended to include the buffers used to define the 'total' pH scale (Hansson, 1973;DelValls and Dickson, 1998). These buffers are solutions of artificial seawater containing equimolal Tris (or THAM, 2-amino-2-hydroxymethyl-1,3-propanediol, CAS 77-86-1) and its conjugate acid TrisH + (DelValls and Dickson, 1998 ) and H + concentration and activity; and in calculations of trace metal speciation in natural waters (e.g., Pierrot and Millero, 2016). For pH, the models can assist in the extension of the total scale to lower salinities, and the development of pH standards for saline waters containing the ions of seawater but in different ratios to those in standard seawater. The use of Pitzer-based speciation models has also been suggested in order to address metrological issues concerning the total pH scale (Dickson et al., 2016).
However, the speciation models referred to above do not yet agree to within experimental uncertainty with available thermodynamic measurements that can be used for validation (chiefly electromotive forces (EMFs) of acidified artificial seawater), nor do they include uncertainty estimates. Indeed, no comprehensive Pitzer model does so.
Quantification of the uncertainty of the model outputs (speciation and solute and solvent activities, and quantities calculated from them) is essential if the speciation models are to be used in the applications noted above. The estimation of uncertainties ideally requires that the variances and covariances of all Pitzer interaction parameters are known, based upon the statistics of the fits to the original thermodynamic measurements and taking account of their experimental uncertainties. The uncertainties associated with each thermodynamic equilibrium constant should also be known. These quantities, together with standard error propagation techniques, allow both the total uncertainty in speciation model outputs to be calculated and also the individual contributions of all parameters and equilibrium constants. The latter are needed to identify the chemical systems for which new measurements are required, or existing data (and the applicable Pitzer interaction parameters) need to determining the influences on the predicted activity coefficients of H + , NH 4 + , and NH 3 in artificial seawater media caused by fractional changes in the molalities of the major ions in order to determine the principal controls of ammonia speciation. However, such approaches are indirect, do not include a consideration of uncertainties, and results are difficult to interpret in complex solutions involving multiple equilibria.
In order to improve current speciation models, and make them more useful for important applications such as pH, it is necessary to be able to estimate both the total uncertainty of any calculated quantity (such as pH, activity, or stoichiometric dissociation constant), and the individual contributions of thermodynamic equilibrium constants and Pitzer model interaction parameters. The fact that these interaction parameters are not independent of one another must also be taken into account. In the next section we describe the inclusion of these capabilities, for the first time, in a Pitzer-based speciation model.

Methods
The two models of artificial seawater involve three chemical equilibria: the dissociation of HSO 4 -(HSO 4 -⇌ H + + SO 4 2-), that of water (H 2 O ⇌ H + + OH -), and the formation of the MgOH + ion pair (Mg 2+ + OH -⇌ MgOH + ). They are solved iteratively to obtain the equilibrium speciation and ion and solvent activities, using techniques that are described in the Supporting Information.
Uncertainties of model-calculated properties can be determined by standard methods of error propagation such as used by Orr et al. (2018) with models of the marine carbonate system if values of the variances and covariances of the Pitzer interaction coefficients, and uncertainties of equilibrium constants in the model are known. There are none available for the interaction parameters in either the  model or that of Clegg and Whitfield (2013) (they are rarely specified in any study), and it would be impractical to attempt to re-derive the models to determine them. It is therefore necessary to adopt a simplified approach to estimating parameter variances and covariances, with the expectation that those which contribute the most to the uncertainty of model-calculated properties can later be evaluated individually. Such a simplified approach should be sufficient to identify those interactions and parameters that need revision (perhaps from new experimental measurements) to improve the accuracy of the models, and to provide estimated or indicative uncertainties in model outputs.
Interaction parameters for pure aqueous electrolyte solutions (β ca activities to a single temperature); (ii) solute activities from EMF measurements of various electrochemical cells, of which the most important for the present work is the Harned cell (see section 4) which yields the activity product aH + ·aCl -. The values of 'mixture' parameters (ζ cc' , ψ cc'a , ζ aa' , and ψ aa'c ) can be determined from the same types of measurement, and also from solubilities of salts in solution mixtures (e.g., Harvie et al., 1984). In this work we assume that, for all single solute solutions and mixtures for which there are interaction parameters, their values have been determined from single datasets of osmotic coefficients (ϕ) subject to the random and systematic errors that are typical of isopiestic measurements. Osmotic coefficients were adopted because a simple approach was needed for the comprehensive uncertainty analysis undertaken here, and because isopiestic measurements are a principal method of activity measurement for solutions of non-volatile electrolytes and other solutes at room temperature and above (Rard and Platford, 1991). These include all of those salts present in artificial seawater. Thus, variances and covariances based on the assumption of osmotic coefficient datasets can be expected to be representative of many of the cation-anion interactions that most affect the solute activity coefficients and speciation in seawater. It is our assessment that osmotic coefficients constitute more than 75% of the total datasets for single solute solutions at 25 o C, but less than 50% for mixtures (although they have been measured for mixtures of some of the solutes that have the highest concentrations in seawater, or those that are important for the calculation of acid-base equilibria. We have estimated variances and covariances of the Pitzer interaction parameters in the  model at the single temperature of 25 o C. We expect the values for the Clegg and Whitfield (1995) (Harvie et al., 1984), and its standard uncertainty was estimated to be ±0.022 in ln(K). Variances of both the equilibrium constants, that of the dissociation of water, and some further details are listed in the Supporting Information.
The approach to the estimation of variances and covariances outlined above corresponds to how Pitzer-based chemical speciation models of solution mixtures are typically developed, but it does not capture all types of model uncertainty. First, there are limits to the maximum molalities to which osmotic and activity can be represented accurately (about 6 mol kg -1 for many electrolytes containing singly charged ions). While this is not directly a limitation at seawater ionic strengths, fits over extended ranges of concentration can result in reduced accuracy at all molalities. Second, mixture parameters ζ ii' and ψ ii'j (where i and i' are dissimilar ions of one sign, and j is an ion of the opposite sign) are often determined from measurements at high molalities, at which their influence is greatest, and errors in the representation of the cation-anion interactions from the pure solution parameters will be reflected in the fitted values of ζ ii' and ψ ii'j . The same problem is apparent when the mixture parameters are used with sets of pure solution parameters different from the ones they were originally determined with. It is important to note that systematic errors will be introduced into model calculations at temperatures other than 25 o C if there are strong interactions for which the variation of the Pitzer parameters (or equilibrium constants) with temperature is not known. Such errors are not accounted for in the parameter variances that we have determined here, and we have confined our analysis of uncertainties to25 o C.
Next we describe the determination of the functions characterising the errors typical of isopiestic experiments, their application within a Monte Carlo approach to determine the variances and covariances of the interaction parameters, and the error propagation method used to calculate the corresponding uncertainty in model-calculated outputs.

Errors and uncertainties of osmotic coefficients
The isopiestic method, described in detail by Rard and Platford (1991)

J o u r n a l P r e -p r o o f
Osmotic coefficients are subject to both random measurement errors that affect each data point differently, and systematic errors that affect the entire dataset consistently. The latter might be caused, for example, by incomplete equilibrium in the isopiestic chamber (although we recognise that the degree of disequilibrium might vary between sample cups). For each measurement i, of molality m(i), in an artificial dataset we simulate both a random and a systematic measurement error.
These are then added to the osmotic coefficient calculated by the model (ϕ mod (i)) to obtain a simulated measurement (ϕ sim (i)): where ϕ mod (i) is the osmotic coefficient calculated by the model for molality m(i). The quantities Δϕ sys (i) and Δϕ rdm (i) are, respectively, measures of the systematic and random errors for the same molality. These quantities are calculated using three empirical parameters (δ sys , ε rdm , and ε rdm ): Δϕ sys (i) = δ sys (1 + 1/(m(i) + 0.03)) Δϕ rdm (i) = (ε rdm + ε rdm /(m(i) + 0.03)) The forms of these equations (notably the use of the denominator (m(i) + 0.03)) reflect the fact that errors in osmotic coefficients from isopiestic experiments increase at low molalities, because of the relationship between water activity and osmotic coefficient. The D and H(i) coefficients in equation (1) are simulated by a pseudorandom number generator (Mersenne twister). They have the following distributions: These expressions have the following meanings: the value of D is drawn, once for each artificial dataset, from a normal distribution (N) with a mean of 0 and variance of 1. A value of H(i) is drawn, once for each measurement i, from a normal distribution with a mean of 0 and variance of π/2. The use of the variance of π/2, rather than unity, reflects that fact that parameter ε rdm is determined as a mean absolute deviation, which must be multiplied by the factor (π/2) 0.5 to obtain the equivalent standard deviation (and (π/2) 0.5 N(0,1) is equal to N(0, π/2)).
In this study, we used constant values of δ sys = 0.00116, ε rdm = 0.003 and ε rdm = 0 in the simulations. They represent the average random and systematic offsets in our analysis of real osmotic coefficients of multiple electrolyte pairs. We illustrate this analysis here using the J o u r n a l P r e -p r o o f Journal Pre-proof electrolytes KCl and NaCl, which have been extensively studied and for which comprehensive evaluations of their thermodynamic properties exist (Archer, 1992;1999). Figure 1a shows the difference between osmotic coefficients of aqueous NaCl determined from measured isopiestic ratios against aqueous KCl (with aqueous KCl osmotic coefficients calculated from Archer, 1999), and those calculated directly from the thermodynamic model of Archer (1992) Figure   1a show Δϕ sys (i) calculated using this result and equation (2). The solid lines show the result using the larger value of δ sys (0.00116) determined using the complete dataset of many electrolyte pairs. Figure 1b shows absolute values of Δϕ from part (a), but with systematic offsets (calculated using individual fits of equation (3) to each dataset) subtracted. The mean value for these datasets is indicated by the dashed line, and the solid line shows the value (ε rdm equal to 0.003) adopted on the basis our study of data the full range of electrolyte pairs, and the assessment of Rard and Platford (1991).
The systematic component of error δ sys was estimated from 24 isopiestic datasets involving combinations of NaCl, KCl and CaCl 2 . Values of δ sys from the majority of the datasets are well represented by a normal distribution, for δ sys within the range ±4×10 −3 . The presence of a few outliers (only 3 out of 24 separate studies) indicates that sometimes larger systematic errors occur.
However, there is at present insufficient evidence to support using any other distribution, and the Gaussian distribution remains is the best option due to the central limit theorem. The random deviations, whose absolute values are shown in Figure 1(b), conform to a normal distribution if 8 values out of a total of 205 are omitted. These values of |Δϕ -δ sys (1 + 1/(m + 0.03))| are all greater than about 3 times the mean shown in the figure (dashed line), and are mostly at low molalities for which very small deviations from equilibrium in the isopiestic apparatus produce large errors in the osmotic coefficient. Data such as these would typically be omittedas obviously erroneous -from fits of osmotic coefficients to obtain Pitzer model parameters.

Variances and covariances of the Pitzer interaction parameters
We first consider pure aqueous solutions containing hypothetical electrolytes consisting of every cation-anion pair in the model for which there are non-zero (β ca (0-2) as follows: (i) for each solute a realistic measurable range of molality m min to m max is first defined.
The value of m min is set to the practical lower limit of isopiestic measurements (0.2 mol kg -1 ), and m max generally to the solubility of the solute (but never greater than 6 mol kg -1 ). (ii) We assign a set of N meas measurements per dataset, evenly spaced from m min to m max with respect to the √m because this yields values that are typical of real experiments. In this work N meas is equal to 50, and we assume a single dataset for simplicity. (The use of multiple datasets significantly reduces simulated uncertainties.) (iii) Next, the osmotic coefficient at each molality is calculated, and sets of systematic and random measurement errors determined using equations (1) to (5) are added to the modelled osmotic coefficients, which are then fitted to obtain a new set interaction parameters. These simulation and fitting steps are repeated many times (typically 10 4 ) and, finally, variances and covariances of the interaction parameters across all steps are evaluated and then saved as a matrix of values. This matrix is later used to propagate the uncertainty associated with the corresponding set of interaction parameters through the Pitzer model to any of its output quantities.
Second, there are ternary interactions for mixtures of two electrolytes with a common ion (parameters ζ cc' and ψ cc'a , or ζ aa' and ψaa'c ). The procedure for determining variances and covariances is similar to that described above. There is now a pair of (m min , m max ) values, and the N meas solution compositions (25) are assigned in five groups each of which has a nearly constant molality ratio of the two solutes (resembling typical experiments). Variances and covariances of (ζ cc' and ψ cc',Cl ) and (ζ aa' and ψ aa',Na ) are determined first, from simulations of electrolyte mixtures containing Cland Na + counterions, respectively. The variances of each ζ cc' determined in this way was then used in further simulations to obtain all ψ cc'a , where anion a is any anion except Cl -. The same approach was applied to obtain ψ aa'c (where c is any cation except Na + ). It is a consequence of this approach, which is analogous to how parameter values themselves are determined, that there are no covariances between different ψ cc'a , or ψ aa'c .
The assumption of a single simulated dataset for pure aqueous solutions represents a worst case (there are more sets of measurements in most cases Details of the simulations, the methods used, and the arrangements of values in the uncertainty matrices, are described in the Supporting Information.

Determination of uncertainties in model outputs
The propagation of uncertainties, to calculate those in model outputs such as EMF, species concentration, or activity, is carried out using a matrix approach as described by Orr et al. (2018).
See their equation 2 and equations A.1 to A.5. The additional quantities that are required, in addition to the variances and covariances of the parameters and logarithms of the equilibrium constants, are the partial derivatives of the model output quantity with respect to each individual interaction parameter and ln(K). These are determined numerically, with a centred finite difference formula using two to six points depending upon the accuracy required. The calculation of both the total variance, and the individual parameter and equilibrium constant contributions to it, are summarised in Appendix B of this work.

Examples: aqueous NaCl and HCl
In the Supporting Information we describe the results of applying the above methods to aqueous NaCl and HCl and compare the calculated uncertainties to several evaluations (including the models presented in this work) of osmotic and activity coefficients at 25 o C. These two solutes were chosen, first, because NaCl is the major dissolved constituent of seawater, and interactions of dissolved Clwith H + (calculated using parameters obtained from thermodynamic properties of aqueous HCl) are central to the calculation of pH. Second, both solutions have been extensively studied and their thermodynamic properties are the subject of a number of critical evaluations. For aqueous NaCl, osmotic and mean activity coefficients from the model of  and from the following studies were compared with values from the critical evaluation of Archer (1992) (which is used in the model of Clegg and Whitfield (1995)): Robinson and Stokes (1970), J o u r n a l P r e -p r o o f Clarke and Glew (1985), Gibbard et al. (1974), and Hamer and Wu (1972). For aqueous HCl, calculated Harned cell EMFs and mean activity coefficients from the models of , Clegg and Whitfield (1995) and the following evaluations were compared with values from Hamer and Wu (1972): Robinson and Stokes (1970), Partanen et al. (2007), and Holmes et al. (1987). The main findings of the above comparisons are as follows: (i) Osmotic and activity coefficients of NaCl lie well within the calculated ranges of uncertainty for ionic strengths up to those equivalent to salinity 45 seawater. This suggests that the uncertainty contributions of NaCl, and similarly well-studied electrolytes, to calculated activities in seawater solutions are likely to be conservatively estimated (i.e., somewhat overestimated) by the model.
(ii) The maximum difference between values of γ NaCl from the  model, and the equation of Archer (1992), is about 0.001 at seawater ionic strengths. A change of this magnitude in the activity coefficient of the Clion would correspond to a change of only about 0.04 mV in the EMF of a Harned cell, thus both models agree very closely with each other for Na + -Clinteractions at seawater ionic strengths.
(iii) For ionic strengths up to that of seawater, the model of  probably represents the interactions of H + with Clin pure aqueous solution to better than 0.1 mV (equivalent to about 0.002 in γ HCl ). The model of Clegg and Whitfield (1995) predicts EMFs that are about 0.05 mV lower (hence γ HCl greater by about 0.001).
(iv) The estimated uncertainty envelope for predicted EMFs and mean activity coefficients of aqueous HCl appears to be too large at most molalities. There are two main reasons for this: first, mean activity coefficients γ HCl are related to the integral of (ϕ -1) with respect to m 1/2 , which means that model-calculated uncertainties (based upon Pitzer parameters determined only from simulated osmotic coefficient data) increase monotonically with molality, as can be seen in Figure S4 Figure S4 is greater than these values by about a factor of 6 for an ionic strength equivalent to seawater of salinity 35 (0.7225 mol kg -1 ), and the greatest difference in γ HCl between the models and critical evaluations and those of the reference is about a factor of 1.5. This is likely to lead to overestimates of the contribution of H + -Clinteractions to the uncertainties of calculated HCl activities and Harned cell EMFs used for comparisons in this work, and must be borne in mind when analysing the results. The method of uncertainty estimation described above should in the future be extended to include EMFs of electrochemical cells as a further data type.
J o u r n a l P r e -p r o o f

Data Used to Assess the Models
Calculations of equilibrium chemical speciation require values of the ion activity coefficients in the seawater or other natural water medium. Electromotive force data for acidified artificial seawater are used to evaluate the accuracy of the models, and the uncertainty propagation methods described above are used for two purposes. First, they are applied to compare the uncertainties in model predictions to those of the experimental data. Second, they are used to determine the solute interactions and equilibrium constants most likely to cause the differences between measured and modelled EMFs, and the solute activities derived from them. The sources of data, summarised in where the aqueous solution contains acidified artificial seawater. Note that the study of Khoo et al. (1977) includes results for artificial seawater not including the SO 4 2ion. The EMF, E (V), of cell A is given by the following expression: where E 0 (V) is the standard EMF of the cell at the temperature T (K) of interest, R (8.31446 J mol -1 K -1 ) is the gas constant, F (96485.332 C mol -1 ) is Faraday's constant, and prefix a denotes activity.
The activity product of the H + and Clions can also be written mH + ·mCl -·γ H ·γ Cl or mH + ·mCl -·γ HCl 2 , where γ i is the activity coefficient of solute species i, and γ HCl is the mean activity coefficient of H + and Clin the aqueous solution (γ HCl is equal to (γ H ·γ Cl ) 0.5 J o u r n a l P r e -p r o o f The value of the dissociation constant, K(HSO 4 -), is known only to within about 5% at 25 o C , and probably less well at other temperatures. However, the accurate representation of this equilibrium in seawater media is central to the interconversion of pH on the

free ([H + ]) and total ([H + ] + [HSO 4
-]) scales, and to an assessment of the empirically determined standard electromotive force E* (Dickson, 1990) used in the definition of the "total hydrogen ion" pH scale (DelValls and Dickson, 1998). This is discussed further in section 7. Values and uncertainties of K(HSO 4 -) from various sources are listed in Table 1.
The published EMF data used in this work for artificial seawater and other solutions are normalized in the original studies to a hydrogen partial pressure of 101.325 kPa according to: where E (V) is the normalized EMF, E meas is the experimentally determined EMF, p atm is the ambient pressure, and pH 2 O is the equilibrium partial pressure of water of the test solution. The contributions to the total uncertainty in E arise from those in p atm , pH 2 O, T, and E meas . The uncertainty in the standard EMF E 0 is also important because it is (E -E 0 ) from which the activity product of H + and Clions is determined (see equation (6) Khoo et al., 1977;Dickson, 1990;and Campbell et al., 1990). This value is consistent with the observation of Dickson (1990) that measurements made at 298.15 K during individual runs tended to agree within 0.05 mV.
In this work we have found that the agreement between measurements of acidified artificial seawater from the different studies is best at 298.15 K, but tends to decrease towards the temperature extremes (and the scatter increases). Also, other factors related to the preparation and treatment of electrodes, and the operation of the cells, can have effects on both uncertainty and accuracy that are not accounted for in this analysis. For example, although the results of Khoo et al. (1977) and Dickson (1990) generally agree well and exhibit behavior that is consistent with the standard uncertainties noted above, the measurements of Campbell et al. (1993) are more scattered and data for two solutions were excluded from our analysis because of consistent deviations (of the order of 0.2 mV and above) at all T.

Assessment of the Models at 25 o C
In this section we compare the models with available EMF data, and use uncertainty information generated using the methods described earlier to identify causes of the differences found.

Calculations of uncertainty contributions to modelled quantities
We first carried out model simulations to determine the relative contributions of the equilibrium constants and interaction parameters in the model to uncertainties of calculated EMFs, and H + and Clactivities and activity coefficients. These simulations were for acidified artificial seawater, both with and without SO 4 2-, at 25 o C. The salinity is 35 in both cases, and the compositions of the solutions are listed in Table 3. In the calculations for acidified artificial seawater we use the recipe of Dickson (1990), rather than that of Khoo et al. (1977), because it is the former that has been used in experiments to define pH on the total scale (DelValls and Dickson, 1998). The model of  was used in all uncertainty simulations.
Two sets of simulations were carried out. In the first set the variances and covariances of parameters whose values are set to zero in the model are also set to zero. These parameters are listed in Tables S2 to S4 in the Supplementary Information and include, for example, ζ HSO4,SO4 and those for interactions between pairs of reacting species such as H + and MgOH + . There are also parameters for interactions that are unknown because of a lack of data from which to determine them. These are J o u r n a l P r e -p r o o f assigned values of zero by default, but may be non-zero. Their variances and covariances are simulated in the same way as described for the other parameters.
The second set of uncertainty simulations is intended to explore the influence of model parameters whose values are unknown, but may not be zero (identified by 'U' in Tables S2 to S4). In this case we have substituted mean parameter values (for charge types corresponding to those of the interacting ions) from Tables A1 and A2 in the Appendix, and set their variances equal to the squares of the standard deviations. We have not attempted to estimate covariances of, for example, ζ ii' and ψ ii'j parameters whose values are generally determined simultaneously. This will tend to increase their contributions to the total calculated uncertainty. This substitution of non-zero parameter values into the model means that the calculated quantitiesboth speciation and activity coefficientswill be different from the base model. However, the differences have been found to be very small. Khoo et al. (1977) replaced SO 4 2with charge-equivalent amounts of Clwhen preparing these solutions. The composition is given in Table 3. Figure 2 shows the difference between measured and modelled EMFs, and the same data as mean activity coefficients of HCl. The principal uncertainty contributions to the calculated (E -E 0 ), at ionic strength 0.7228 mol kg -1 , are displayed in Figure 3. The partial derivatives of the calculated (E -E 0 ) with respect to each parameter are shown in Figure S1 of the Supporting Information. Both these results show that H + -Cland Na + -Clinteractions, and the parameters ζ H,Na and ψ H,Na,Cl , dominate the uncertainty in the modelled EMF, which is to be expected given that NaCl constitutes 86 % (on a molar basis) of the dissolved solutes.

Acidified artificial seawater without sulphate
As is also the case for H + -Clinteractions, the parameters ζ H,Na and ψ H,Na,Cl are determined from EMF data, and so their uncertainty contribution may be overestimated for reasons given in section 3.4.
However, this does not change the conclusion that the three sets of interactions dominate the total calculated uncertainty. Separate calculations of EMFs and uncertainty contributions using the averaged parameter values (see the Appendix) in place of those that are unknown yielded identical results to those shown. This is because, for this system, they involve only the ions MgOH + and OHwhich have no influence on speciation or activities in these acidic solutions.
The measured EMFs and γ HCl in Figure 2 lie within the uncertainty envelope of the modelled values. As described above, this is likely to be too large given that H + -Cland H + -Na + -Clinteractions are calculated to contribute the most (90% of the total variance, Figure 3 The comparisons in Figure 2 show that the models of  and Clegg and Whitfield (1995) yield very similar results for these artificial seawater solutions. This is despite that fact that Figure S4 of the Supporting Information indicates that EMFs of pure aqueous HCl solutions, calculated by the two models, differ by as much as 0.05 mV to 0.1 mV over the ionic strength range of the measurements of Khoo et al. (1977). The contributions of Na + -Clinteractions to ln(γ HCl ) in the artificial seawater were very similar for both models, which is also consistent with the result for pure aqueous NaCl ( Figure

Acidified artificial seawater
The compositions of the artificial seawaters used in the measurements of Khoo et al. (1977), Dickson (1993), and Campbell et al. (1993), are listed in Table 3 Table 1), and it is the uncertainty for this value that has been used in our model calculations.
As a test, we set K(HSO 4 -) equal to 0.01086 mol kg -1 in the model, and recalculated both Δ(E -E o ) and γ HCl(stoic.) , see Figure 7 and the dashed lines in Figure 6. The changes in Δ(E -E o ) seen in  (Sippola and Taskinen, 2014; and references therein) tends to support a value of the dissociation constant greater than 0.0105 mol kg -1 , and we conclude that the Pitzer model of aqueous H 2 SO 4 should be revised.

Comparisons with Data from 0 o C to 45 o C
We have evaluated the accuracy of the two models at temperatures other than 25 o C by comparing calculated (E -E 0 ) with all the measurements of Khoo et al. (1977), Dickson (1990) and Campbell et al. (1993). The results, in Figures 8 to 11, are plotted as averaged deviations against temperature and various measures of concentration (ionic strength, salinity, and total molality of H + ).
No estimates of model uncertainties are shown, because their current treatment is restricted to 25 o C as described earlier.
The results for solutions not containing SO 4 2-(from Khoo et al., 1977) are shown in Figure 8.
where γ HCl is the mean activity coefficient of HCl in the solution. Note that, if equation (12) (12) is equivalent to a standard potential of the cell (E* / V) for the temperature and salinity of interest, and is obtained experimentally by the extrapolation of measured EMFs to zero HCl molality (Dickson, 1990). The standard potential is given by: where the superscript (tr) means that the value is for trace (zero) HCl in artificial seawater. In such a solution the molality of HSO 4 is so small relative to that of SO 4 2that the total SO 4 2molality, denoted by superscript (T), can be used without loss of accuracy. This definition is equation (13) of Dickson (1990).
The uncertainty profile for E*, calculated using the model of  for salinity 35, is shown in Figure 12. It is very similar to that of the EMF of acidified artificial seawater  (1977), Dickson (1990), and Campbell et al. (1993) yield values of E* that are consistent to within about 0.1 mV. This is much smaller than the estimated standard deviation in E* calculated directly by the model. The contribution of H + -Clinteractions to this uncertainty (about 26% of the variance) may be too high by a factor of 2 or more, so that the true standard deviation is smaller, but the contribution of ln(K(HSO 4 -)) dominates. Obtaining a more accurate value than is currently used in the modelsand reducing its uncertaintyis essential and appears to be the key to a more accurate prediction of E*, as it is for the EMFs of the acidified seawaters generally. This is examined further below.
As noted above, Dickson (1990) obtained values of E* empirically, by fitting the following quadratic equation: where E (V) is the measured EMF, mHCl is the molality of added HCl in the solution, mClis the total chloride molality, and b 0-2 are the fitted coefficients. The coefficient b 0 is equivalent to E*, for which values are listed in Table 3 of Dickson (1990). The fits of the data, for salinity 35 and temperatures 25 o C and 5 o C, are reproduced in Figure 13 where the plotted quantity E' is the left hand side of equation (14)  We conclude that, first, the results in general confirm that the empirical fit of equation (14) yields satisfactory estimates of E*, and that the uncertainty is dominated by that of the measurements themselves and not artifacts of the fitting procedure. Second, the very close agreement of the model with the data, although not yet attained at other salinities, gives encouragement that an accurate model of this chemical system is achievable. Third, such a model should include a new Pitzer model parameterisation of the aqueous H 2 SO 4 (including revised values of K(HSO 4 -)).

Recommendations for Future Work
Here we summarise new measurements, and reassessments of existing data, that are needed to improve the models of artificial seawater at all temperatures, with the aim of achieving agreement J o u r n a l P r e -p r o o f with measured EMFs to within, or close to, experimental uncertainty. This is needed both for applications in the definition and determination of pH, and as the basis for a larger model of speciation in natural waters. The profiles in Figures 3, 5, and 13 identify the major contributors to the uncertainty of model predictions of EMF and E* at 25 o C, and therefore those interactions and equilibrium constants that must be calculated most accurately in a successful model. These are identified in Table 5, together with the predicted quantities for which they are most important. There are two general points to bear in mind. First, improved agreement between measured and modelled EMFs requires the calculation of activities and activity coefficients to better than 1% accuracy. This is likely to require that the key ternary mixture parameters -notably ζ H,Na and ψ H,Na,Cl for H + -Na + -Clinteractions, but also interactions for other ionsshould always be determined using the same parameters for the cation-anion binary interactions (β ca (1-2) , C ca (0,1) ) as are used in the models.
Second, the analysis in this study is primarily for the single temperature of 25 o C, and there are a number of solute interactions that have only been quantified at this temperature. The use of these parameters in calculations for other temperatures will result in larger errors and uncertainties than would be expected from their positions in the uncertainty profiles. This is apparent in the larger standard deviations, and poorer agreement between models and data shown at the temperature extremes (see panels (b) of Figures 9 to 11)). The parameters can be identified in the tables in the Supplementary Information, and those likely to be most significant are included in Table 5. The contents of the table are discussed below.

Aqueous HCl, and H + -Na + -Clsolutions
The parameter ζ H,Na in the Waters and Millero (2013) model has been revised in this work as noted earlier, by fitting to EMFs of H + -Na + -Claqueous solutions measured by Macaskill et al. (1977). These data were also used by Clegg and Whitfield (1995 Whitfield are -0.05 mV to +0.02 mV. The pattern of deviations in (E -E 0 ), with respect to mHCl in the mixtures, differs between the two models. We recommend that both ζ H,Na and ψ H,Na,Cl parameters be revised, based on a larger set of available measurements. There are numerous sources of early data cited by Harned (1959), the work of Hawkins (1932) at 25 o C for 4 mol kg -1 to 6 mol kg -1 ionic strength solutions, and the study of Macaskill et al. (1977) already mentioned. The work of White et al. (1980), who have measured EMFs of H + -Na + -Mg 2+ -Clsolutions, is also relevant. If the accuracy with which EMFs of H + -Na + -Clsolutions can be reproduced is limited by the model-calculated interactions for aqueous HCl, as appears possible from our calculations, then revisions could be J o u r n a l P r e -p r o o f carried out using sources of data cited by Hamer and Wu (1972), Carslaw et al. (1995), and Partanen et al. (2007) for aqueous HCl. Interactions for Na + -Clare also noted in the lower section of Table 5. The uncertainty profile for EMFs of acidified artificial seawater ( Figure 5) shows that Na + -Clinteractions have about the same contribution as those for H + -Na + -Cl -(and about a half that of H + -Cl -). A revision of the model of Waters and Millero to include parameters from the critical review of Archer (1992) should be investigated. We note that there are small differences between the osmotic and activity coefficients of Archer (1992) and those from Clarke and Glew (1985) ( Figure S5 of the Supporting Information), but it is unclear what significance they might have for the calculation of natural water properties. Clarke and Glew (1985) used a equation equivalent to that of Pitzer (with parameters equivalent to β (0) Na,Cl , β (1) Na,Cl , and C (0) Na,Cl ) but extended with two further terms.  (Pitzer et al., 1977;Sippola, 2013). The model of Clegg et al. (1994) for aqueous H 2 SO 4 should be rederived, based upon higher values of K(HSO 4 -), and including both new data that have become available since their study and some older work that was overlooked. These data include osmotic coefficients, EMFs of two electrochemical cells, and heats of dilution (J. A. Rard, pers. comm.). It is recommended that the values of K(HSO 4 -) obtained by Dickson et al. (1990) are used. We note that the model of Sippola and Taskinen (2014) used K(HSO 4 -) equal to 0.0119 mol kg -1 at 25 o C.

The dissociation constant of HSO
However, the values of the dissociation constant were determined as a part of the fit of the overall model in their study, rather than being determined separately, and their model did not include some enthalpy measurements or any heat capacity measurements as a part of the fitted dataset. Both types of data can be used to quantify the variation of solvent and solute activities with temperature (e.g., Pitzer, 1991).

Other Interactions
The parameters and solutions listed under this heading in Table 5, with the exception of NaCl, are specified only at 25 o C in the model of . (In the model of Clegg and Whitfield (1995) both ζ Cl,SO4 and the parameters for Na + -HSO 4 interactions vary with T, but newer data have become available.) Consequently the uncertainty profiles for EMFs and E* are not likely to represent accurately the contributions of this group of parameters to the overall uncertainty for other temperatures. The contributions are likely to be larger.
The two mixture parameters ζ Cl,SO4 and ζ Cl,HSO4 are generally determined in combination with the parameters ψ Cl,SO4,c and ψ Cl,HSO4,c respectively (where c is H + or a major seawater cation such as Na + or Mg 2+ ). The value of ζ Cl,SO4 equal to 0.020 at 25 o C in both models (Harvie and Weare, 1980;Pitzer and Kim, 1974) is based upon isopiestic measurements of Wu et al. (1968). Rard et al. (2003), using a larger range of data (and different cation-anion parameters for Na + -Cland Na + -SO 4 The Na + -HSO 4 interaction parameters in the model of  come from the study of Hovey et al. (1993), and are based upon the isopiestic measurements at 25 o C of Na 2 SO 4 -H 2 SO 4 aqueous solutions by Rard (1989Rard ( , 1992. The mixture parameters determined by Hovey et al. for this system were not adopted by Waters and Millero. Values of the Na + -HSO 4 interaction parameters in the model of Clegg and Whitfield (1995) are based upon the same isopiestic data, and also EMF measurements (see their section 4.1), and their treatment is selfconsistent. In their model the variation of the Na + -HSO 4 parameters with temperature were obtain by analogy, and obtained by fitted enthalpies of aqueous NaClO 4 . The EMF measurements of Pierrot et al. (1997), of aqueous HCl-Na 2 SO 4 solutions from 0 to 50 o C, are also relevant.

J o u r n a l P r e -p r o o f
The value of ζ Cl,HSO4 in both models (from Harvie et al., 1984) appears to be based upon the EMFs of aqueous mixtures at 25 o C measured by Storonkin et al. (1967). In contrast, Pierrot et al. (1997) found that ζ Cl,HSO4 could be set to zero at all temperatures. It may be possible to determine whether this parameter is indeed redundant from the EMF measurements of HCl-MgSO 4 aqueous mixtures noted above.
The interactions for H + -Mg 2+ -Cland H + -Ca 2+ -Clsolutions are not mentioned in Table 5. We  Table A7 and the data (EMFs, over a range of temperatures) from which they had been derived. It is therefore recommended that these parameters be redetermined as functions of temperature.
Finally, several of the aqueous mixtures discussed above have ions, and therefore parameters, in common. Self-consistency of the parameter set is essential for both the accuracy and reliability of the model, and the suggested improvements will entail other model revisions (e.g., ζ Cl,SO4 and ψ Cl,SO4,Na determined from data for aqueous NaCl-Na 2 SO 4 also occur in the model for H + -Na + -Cl --HSO 4 --SO 4 2solutions). This will be true of the revised model of aqueous H 2 SO 4 , suggested in section 8.2, which will require redetermination of interaction parameters for all acid sulphate mixtures.

Uncertainty calculations
In this study we have applied a highly simplified approach in which all interaction parameters are assumed to be determined from a single type of thermodynamic data (osmotic coefficients determined by isopiestic equilibrium) at a single temperature. A more realistic approach should for the key interactions that contribute the most uncertainty to model outputs, by refitting to original measurements.

J o u r n a l P r e -p r o o f
The extension of the treatment of Pitzer parameter variances and covariances to other temperatures can partly be addressed by (iv) above, for a small number of individual solutes and mixtures. For the more general case it will be necessary to treat differently parameters that are known only at a single temperature (and for which variances should be expanded if they are used at other temperatures), and those for which there are data from which the variation of the parameters with temperature is known.

Summary and Conclusions
Here we have implemented the  and Clegg and Whitfield (1995) chemical speciation models of artificial seawater. The propagation of uncertainties of the interaction parameters and thermodynamic equilibrium constants has been included in the models, for the first time, and a simplified method of estimating their variances and covariances has been developed. This has enabled the estimation of indicative uncertainties in model outputs such as solute activity, speciation, and EMF, and allows the identification of the interaction parameters that contribute most to the overall uncertainty and which should therefore be the subject of future efforts to increase the accuracy of the models.

Comparisons made in sections 5 and 6 show that model calculated EMFs do not yet agree
with measured values to within experimental uncertainty for acidified artificial seawater of various salinities. The calculated total uncertainties in modelled EMFs are large enough that the difference between measured and modelled EMFs at 25 o C lies within the uncertainty envelope of the modelled EMFs in most cases (see Figures 2 and 4). However, this is probably too large given that the variance contribution of H + -Clinteractions is overestimated because of the assumption of osmotic coefficientsrather than EMFsas the underlying data from which they were determined (section 3.4, and Supporting Information). It is important that the treatment of uncertainties be extended to include EMF measurements as a second data type, as noted above.
Total uncertainties in calculated EMFs, including the standard EMF E*, are dominated by contributions from ln(K(HSO 4 -)) and a small number of interaction parameters, mainly those for H + -Cl -. This will greatly assist in improving the accuracy of the model and we have recommended additional work to be carried out.
Calculations at two temperatures of the standard EMF, E*, used in the definition of the total pH scale suggest both that the empirical extrapolation used to determine its values (Dickson, 1990) is satisfactory and that an accurate model of artificial seawater mediaone which agrees with measured EMFs to within experimental uncertaintyis attainable.

J o u r n a l P r e -p r o o f
An extension of the model to include the species Tris and TrisH + is recommended. This will enable the pH of Tris buffers to be modelled directly, and some of the approximations made in the definition of the total pH scale quantified. These include, for example, the effect of the difference term that arises from the use of E* (determined in pure artificial seawater) for solutions containing the buffer. An extension will also help extend the total pH scale to lower salinity, develop definitions of pH for natural water compositions other than that of seawater, and carry out conversions between total pH, and free H + ion concentrations and activities.

Supporting Information, and Software
There are 10 numbered documents of supporting information. Document zero summarises the contents of numbers one to nine, and lists the tables and charts that appear in each one. The subjects covered are: the simulation of uncertainties; values of variances and covariances for interactions and equilibrium constants for HSO 4 dissociation, and MgOH + formation; values of the Pitzer parameters and equilibrium constants for both models; and calculated equilibrium solute molalities and activity coefficients for program verification. Document number four contains an application of the uncertainty simulation method to osmotic and activity coefficients of pure aqueous NaCl and HCl. It is expected that software tools incorporating the model of , as amended in this work, will be released in late 2022. See website marchemspec.org for future announcements, or contact the corresponding author.

Acknowledgements
The

Pitzer interaction parameters
For interactions between cation c and anion a. Not all of these may be used, e.g.,  (2) ca is usually for 2:2 charge types only (e.g., CaSO 4 ), and is set to zero otherwise.  ca ,  (2) ca , ca Coefficients associated with the ionic strength terms in the functions that use parameters  (1) ca , (2) ca , and C (1) ca , respectively.  cc' ,  aa' For interactions between dissimilar cations c and c', and between dissimilar anions a and a', respectively.  cc'a ,  aa'c For interactions between anion a and dissimilar cations c and c', and between cation c and dissimilar anions a and a', respectively.  nc ,  na  For interactions between neutral solute n and cation c, and between neutral solute n and anion a, respectively.  nn , μ nnn  For the self-interaction of neutral solute n.  nca For the interaction between neutral solute n, cation c and anion a. Other symbols used in the text aX Activity (molality basis) of species X, equivalent to mX·γ X where γ X is the activity coefficient of X. Δ r C p o Heat capacity change for a reaction at constant pressure (J mol -1 K -1 ). D Coefficient by which the systematic error parameters Δϕ sys (i) are multiplied to obtain the contributions to the simulated osmotic coefficients in a dataset. See equation (4), and the explanation that follows.

E
Electrode potential (V) in a Harned cell. E 0 Standard electrode potential (V) of a Harned cell.
The difference between model-calculated and measured values of (E -E 0 ), see definitions above and equation (6).

E*
The standard potential (V), on a total H + basis, defined by equation (13). It is obtained empirically by extrapolating Harned cell potentials to zero HCl molality in an artificial seawater of a specified composition and temperature (equation 14), and can also be calculated directly using speciation models.

H(i)
Coefficient by which the random error parameter Δϕ rdm (i) is multiplied to obtain the contribution to each simulated osmotic coefficient in a dataset. See equation (5), and the explanation that follows.

I
Ionic strength, on a molality basis (0.5Σ i m i |z i | 2 , where z i is the charge on ion i and the summation is over all ions). 'Trace' value of the stoichiometric bisulphate dissociation constant in artificial seawater, which is its value in the limit of zero added HCl (i.e., for a composition of pure artificial seawater at some specified salinity and temperature). The pH on the total scale (expressed on a moles per kg of solution basis), defined by Dickson (1990) and DelValls and Dickson (1998). pX The vapour pressure of species X. p atm Atmospheric pressure. R The gas constant (8.31446 J mol -1 K -1 ) S Salinity. Strictly, this definition refers to seawater only and for the artificial seawaters considered in this work S is a nominal salinity. T Temperature (K).

u[p]
Standard uncertainty of a measured or predicted property 'p'. The analogous property U [p], which occurs in the Supporting Information, is the expanded uncertainty (equal to 2u[p]).  X Activity coefficient of species X, on a molality basis.  HCl (tr)


Trace value of the mean activity coefficient of HCl (in pure artificial seawater, containing no added HCl).  HCl(stoich.)  The stoichiometric mean activity coefficient of HCl. This is based upon the total H + and Clmolalities in solution (irrespective of equilibria such as the formation of HSO 4 -). Thus the H + •Clactivity product in solution, which is given by mH + •mCl -• HCl 2 on a free ion basis, can also be expressed as (mH + + mHSO 4 -)•mCl -• HCl(stoich.) 2 in the solutions of interest in this study. ϕ Molal osmotic coefficient of a solution.

Δϕ
The difference between measured osmotic coefficients of aqueous NaCl (from various sources) and values calculated using the thermodynamic model of Archer (1992), and shown in Figure 1 of this work. ϕ sim (i) Simulated osmotic coefficient measurement i, consisting of three components: a model-calculated true value (ϕ mod (i)), plus a systematic error (a term in Δϕ sys (i)) and a random error (a term in Δϕ rdm (i)), see equation (1). δ sys , ε rdm , and ε rdm The empirically determined parameters used in the calculation of systematic and random error contributions to a simulated osmotic coefficient measurement ϕ sim (i) above (see equations (2) and (3)).

Appendix A. Pitzer Interaction Parameters That Are Equal to Zero
An aqueous solution containing six cations and four anions can potentially require twenty four sets of pure solution cation-anion parameters (β ca (0-2) and C ca (0-1) ), and over one hundred cation-cationanion and anion-anion-cation parameters. In solutions containing uncharged, neutral, solutes there are further interaction parameters for neutral-neutral, neutral-ion, and neutral-cation-anionparameters (Pitzer, 1991). However, a large number of parameters are set to zero, and generally fall into the following categories: (i) one or more of the potentially five pure solution ion-interaction parameters for an individual cation-anion pair (β ca In the treatment of the pure solution parameters we have not included the C ca (0) and C ca (1) parameters, but have recognised the correlation of β ca (0) and β ca (1) values for 1:1, 2:1, and 1:2 electrolytes (Table A1). The ranges of values of the individual mixture parameters in Table A2 are wide, even though we have excluded from consideration extreme values that involve either very large ions, or small ions which have a high charge density (such as Li + ). Consequently many of the mean values listed in Table A2 differ from zero by less than one standard deviation.

Appendix B. Calculation of Variances of Model Outputs
This requires, first, a row vector (J) of partial derivatives of the model-calculated output with respect to each individual interaction parameter, and natural logarithm of each equilibrium constant (ln(K)).
These are determined numerically, as noted in the text. Second, there is a square matrix C containing the variances and covariances of all parameters and ln(K), which is assembled from the values determined as described earlier (and in the same order as the partial derivatives), thus: where y is the calculated model output for the solution of interest (for example an EMF, mH + , or a mean activity coefficient), and p (1-n) are the individual Pitzer model parameters and logarithms of equilibrium constants. The matrix C is: .2) J o u r n a l P r e -p r o o f where ζ 2 (p i ) is the variance of parameter (or ln(K)) i, and ζ(p i , p j ) is the covariance of parameters p i and p j . In the Supporting Information the variances and covariances of groups of interaction parameters are presented as individual small matrices (for example a 5x5 matrix for the parameters β ca (0) , β ca (1) , β ca (2) , C ca (0) , and C ca (1) for cation c and anion a). To create matrix C, all of these small matrices, and the individual variances of each ln(K), are combined such that: (i) the variances (ζ 2 (p i )) form the diagonal of C; (ii) the order of the parameters and ln(K), both along rows and down columns, is the same as in row vector J. All other matrix values are set to zero.
The total variance in the model output y, ζ 2 (y), due to the uncertainties in the interaction parameters and ln(K), is obtained by two matrix multiplications: where D is a matrix with a single row and n columns, and: where J T is the transpose of J. It is also valuable to obtain contributions of individual ln(K), interaction parameters, and groups of parameters to the total variance and to be able to separate these , C ca Tables (1 to 5, followed by A1 and A2 for Appendix A)  Marshall and Jones (1966) 0.01017 0.0002 to 0.0006 expt. Young et al. (1978) 0.0103 -expt. Matsushima and Okuwaki (1988) 0.01039 0.00018 expt. Mussini et al. (1989) 0.01043 0.0002 expt. Mussini et al. (1989) 0.01086 0.0005 expt. Dickson et al. (1990) 0.0105 -model b Pitzer et al. (1977) 0.01036 -model c Hovey and Hepler (1990) 0.0105 -model d Clegg et al. (1994) 0.01058 -model e Knopf et al. (2003) 0.0119 0.0012 model f Sippola and Taskinen (2014) 0.0105 0.0005 model g this work Notes: The 'Type' column indicates whether the value of K(HSO 4 -) is determined primarily from experimental measurements that, directly or indirectly, yield values of the equilibrium constant ('expt.') or are determined by the application of the Pitzer activity coefficient model to activity and thermal data ('model'). a Determined from the application of a model to literature data for solubilities. b Largely determined by the application of the Pitzer model (without unsymmetrical mixing terms) to literature data for emfs. The quoted value of 0.0105 is favoured, but an alternative model treatment yielding a value of 0.0120 was also found to be possible. c These authors adopted the equation of Pitzer et al. (1977) for K(HSO 4 -), and added a Δ r C p o (heat capacity) term. d Clegg et al. (1994) adopted equation (6) of ) (terms in p 1 to p 5 ), but added a fixed increment in log 10 (K(HSO 4 -)) at all T to yield exactly 0.0105 mol kg -1 at 25 o C. e Knopf et al. (2003) fitted a range of activity and thermal data, and their own experimentally determined degrees of HSO 4 dissociation. Osmotic coefficients and vapour pressures were omitted, and some of the emfs (yielding H 2 SO 4 activities) are known to be inaccurate (Rard and Clegg, 1995). f Sippola and Taskinen (2014) fitted a range of activity and thermal data (but omitting heat capacities), and their expression for K(HSO 4 -) as a function of temperature lacks a Δ r C p o (heat capacity) term. g The model of Clegg et al. (1994) is used, with the uncertainty taken from the study of Dickson et al. (1990).  Khoo et al. (1977) 20.31, 34.99, 44.55 0.413 -0.929 5 -40 HCl + ASW Khoo et al. (1977) 5 -45 0.100 -0.939 0 -55

J o u r n a l P r e -p r o o f
HCl + ASW Campbell et al. (1993) 5 -45 0.100 -0.939 0 -45 HCl + ASW Dickson (1990) Notes: a These are formal ionic strengths that do not take into account any ion-pairing (see Khoo et al., 1977) or the formation of HSO 4 in the solutions of artificial seawater with added HCl. b Artificial seawater is denoted by ASW.
J o u r n a l P r e -p r o o f Notes: The compositions in the first two columns of molalities are from Khoo et al. (1977), and the third is from Dickson (1990 J o u r n a l P r e -p r o o f Notes: The data of Dickson (1990), Khoo et al. (1977), and Campbell et al. (1993) Dickson et al., 1990). c The calculations were carried out using the Clegg and Whitfield (1995) model. J o u r n a l P r e -p r o o f Table 5 Interactions and equilibrium constants that need reassessment.

Solutions
Parameters or equilibrium constants  Archer (1992), and later work by Archer and Carter (2000), is a likely alternative to the parameters used by . e These parameters are from isopiestic measurements of water activity (in the model of Clegg and Whitfield they vary with temperature). f The existing values in the model of Waters and Millero are from isopiestic (ζ Cl,SO4 , Wu et al., (1968)) and EMF (ζ Cl,HSO4 , Storonkin et al. (1967) measurements). In the model of Clegg and Whitfield parameter ζ Cl,SO4 varies with temperature. g The values of these parameters in the Waters and Millero model are from isopiestic measurements (Rard and Clegg, 1999) and those in the Clegg and Whitfield model (obtained by Harvie et al., 1984) are from EMF measurements yielding stoichiometric mean activity coefficients of H 2 SO 4 (Harned and Sturgis, 1925). In addition, Clegg and Whitfield assumed that the parameters had the same temperature coefficients as those of Mg(ClO 4 ) 2 . J o u r n a l P r e -p r o o f Notes: Values of the parameters were taken from Kim and Frederic (1988); Pitzer and Kim (1974); Pitzer (1991); Clegg and Whitfield (1995); Clegg and Brimblecombe (1988a,b); Harvie et al. (1984); ; Millero and Pierrot (1998) (2000); Lodeiro et al. (2021). Duplicate values were removed from the dataset. a Where there are two or more different values for the same interaction parameter (for the same ions) each was given a reduced weight (1/2 in the case of two values, 1/3 in the case of three, etc.). Values of the mean that differ from zero by more than one standard deviation are marked by '* '. b The standard deviation.
c The range of parameter values that were used to determine the mean.  Symbols: dot - Robinson (1945); cross - Gordon (1943); circle - Robinson and Sinclair (1934); square - Scatchard et al. (1938); star - Kirgintsev and Luk'yanov (1967); triangle - Janis and Ferguson (1939). Lines: solid -systematic offset calculated from equation (2) using the mean |δ sys | for many datasets (not shown); dashed -the same as the solid line but calculated using the mean |δ sys | for the datasets for NaCl/KCl shown in this figure. (b) The absolute value of Δϕ less the expression for systematic offset (equation (2), fitted to each dataset individually). Symbols: the same sources as in (a). Lines: solid -average random offset (ε rdm ) calculated from results for many datasets (not shown); dashedthe same as the solid line but calculated only for the datasets for NaCl/KCl shown in this figure.   Millero (2013)), plotted against salinity (bottom axis) and ionic strength (I) (top axis). Symbols: circlemeasurements of Dickson (1990); dotmeasurements of Khoo et al. (1977); trianglemeasurements of Campbell et al. (1993). The lighter shaded area shows the total uncertainty in the calculated values of (E -E o ), and is centered on the zero line. The inner, darker, shaded area is the uncertainty attributed to that of the thermodynamic dissociation constant of HSO 4 -. The estimated uncertainty in the measured values of (E -E 0 ) (i.e., +/-one standard deviation) is indicated on the plot.   ) with values calculated using the models of  and Clegg and Whitfield (1995). The data are from Khoo et al. (1977). ( The meanings of the symbols and bars are the same as in (a). Note that the results for the Clegg and Whitfield (1995) model are offset to the right of those of  in plots (a) Clegg and Whitfield (1995) model are offset to the right of those of  in both plots to aid legibility. The estimated uncertainty in the measured (E -E 0 ) (i.e., +/-one standard deviation) is indicated on the plots. Figure 10. Comparison of measured EMFs (and γ HCl derived from them) of acidified artificial seawater with values calculated using the models of  and Clegg and Whitfield (1995). The data are from Dickson (1990).  Clegg and Whitfield (1995) model are offset to the right of those of  in these plots to aid legibility. (d) Measured and calculated γ HCl at two temperatures, and salinity 35, plotted against total H + ion molality (mH + (total)). Symbols: solid square -0 o C (0.06 has been subtracted from all values); open square -40 o C. Lines: solidmodel of ; dashedmodel of Clegg and Whitfield (1995). The estimated uncertainties in the measured y variable (i.e., +/-one standard deviation) are indicated on the plots. Figure 11. Comparison of measured EMFs (and γ HCl derived from them) of acidified artificial seawater with values calculated using the models of  and Clegg and Whitfield (1995). The data are from Campbell et al. (1993). The meanings of the symbols and lines in all plots (a) to (d) are the same as those in Figure 10.