Potentiometric measurements

A free concentration [A] or activity {A} of a species A is measured by means of an ion selective electrode such as the glass electrode. If the electrode is calibrated using activity standards it is assumed that the Nernst equation applies in the form

${\displaystyle E=E^{0}+{\frac {RT}{nF}}\ln \mathrm {\{A\}} }$

where E⁰ is the standard electrode potential. When buffer solutions of known pH are used for calibration the meter reading will be a pH.

${\displaystyle \mathrm {pH} ={\frac {nF}{RT}}\left(E^{0}-E\right)}$

At 298 K, 1 pH unit is approximately equal to 59 mV.^[2]

When the electrode is calibrated with solutions of known concentration, by means of a strong acid–strong base titration, for example, a modified Nernst equation is assumed.

${\displaystyle E=E^{0}+s\log _{10}\mathrm {[A]} }$

where s is an empirical slope factor. A solution of known hydrogen ion concentration may be prepared by standardization of a strong acid against borax. Constant-boiling hydrochloric acid may also be used as a primary standard for hydrogen ion concentration.

Spectrophotometric measurements

NMR chemical shift measurements

Chemical exchange is assumed to be rapid on the NMR time-scale. An individual chemical shift δ is the mole-fraction-weighted average of the shifts δ of nuclei in contributing species.

${\displaystyle {\bar {\delta }}={\frac {\sum x_{i}\delta _{i}}{\sum x_{i}}}}$

Example: the pK_a of the hydroxyl group in citric acid has been determined from ¹³C chemical shift data to be 14.4. Neither potentiometry nor ultraviolet–visible spectroscopy could be used for this determination.^[3]

Calorimetric measurements

Simultaneous measurement of K and ΔH for 1:1 adducts is routinely carried out using isothermal titration calorimetry. Extension to more complex systems is limited by the availability of suitable software.

The competition method

The competition method may be used when a stability constant value is too large to be determined by a direct method. It was first used by Schwarzenbach in the determination of the stability constants of complexes of EDTA with metal ions.

For simplicity consider the determination of the stability constant ${\displaystyle K_{AB}}$ of a binary complex, AB, of a reagent A with another reagent B.

${\displaystyle K_{AB}={\frac {[AB]}{[A][B]}}}$

where the [X] represents the concentration, at equilibrium, of a species X in a solution of given composition.

A ligand C is chosen which forms a weaker complex with A The stability constant, K_AC, is small enough to be determined by a direct method. For example, in the case of EDTA complexes A is a metal ion and C may be a polyamine such as diethylenetriamine.

${\displaystyle K_{AC}={\frac {[AC]}{[A][C]}}}$

The stability constant, K for the competition reaction

${\displaystyle AC+B\leftrightharpoons AB+C}$

can be expressed as

${\displaystyle K={\frac {[AB][C]}{[AC][B]}}}$

It follows that

${\displaystyle K_{AB}=K\times K_{AC}}$

where K is the stability constant for the competition reaction. Thus, the value of the stability constant ${\displaystyle K_{AB}}$ may be derived from the experimentally determined values of K and ${\displaystyle K_{AC}}$ .

The chemical model

The chemical model consists of a set of chemical species present in solution, both the reactants added to the reaction mixture and the complex species formed from them. Denoting the reactants by A, B..., each complex species is specified by the stoichiometric coefficients that relate the particular combination of reactants forming them.

${\displaystyle {\ce {{{\mathit {p}}A}+{\mathit {q}}B\cdots <=>A_{\mathit {p}}B_{\mathit {q}}\cdots }}}$: ${\displaystyle \beta _{pq\cdots }={\frac {[{\ce {A}}_{p}{\ce {B}}_{q}\cdots ]}{[{\ce {A}}]^{p}[{\ce {B}}]^{q}\cdots }}}$

When using general-purpose computer programs, it is usual to use cumulative association constants, as shown above. Electrical charges are not shown in general expressions such as this and are often omitted from specific expressions, for simplicity of notation. In fact, electrical charges have no bearing on the equilibrium processes other that there being a requirement for overall electrical neutrality in all systems.

With aqueous solutions the concentrations of proton (hydronium ion) and hydroxide ion are constrained by the self-dissociation of water.

${\displaystyle {\ce {H2O <=> H+ + OH-}}}$: ${\displaystyle K_{\mathrm {W} }^{'}={\frac {[H^{+}][OH^{-}]}{[H_{2}O]}}}$

With dilute solutions the concentration of water is assumed constant, so the equilibrium expression is written in the form of the ionic product of water.

${\displaystyle K_{\mathrm {W} }={\ce {[H+]}}[{\ce {OH-}}]\,}$

When both H⁺ and OH⁻ must be considered as reactants, one of them is eliminated from the model by specifying that its concentration be derived from the concentration of the other. Usually the concentration of the hydroxide ion is given by

${\displaystyle [{\ce {OH-}}]={\frac {K_{{\ce {W}}}}{[{\ce {H+}}]}}\,}$

In this case the equilibrium constant for the formation of hydroxide has the stoichiometric coefficients −1 in regard to the proton and zero for the other reactants. This has important implications for all protonation equilibria in aqueous solution and for hydrolysis constants in particular.

It is quite usual to omit from the model those species whose concentrations are considered negligible. For example, it is usually assumed then there is no interaction between the reactants and/or complexes and the electrolyte used to maintain constant ionic strength or the buffer used to maintain constant pH. These assumptions may or may not be justified. Also, it is implicitly assumed that there are no other complex species present. When complexes are wrongly ignored a systematic error is introduced into the calculations.

Equilibrium constant values are usually estimated initially by reference to data sources.

Speciation calculations

A speciation calculation is one in which concentrations of all the species in an equilibrium system are calculated, knowing the analytical concentrations, T_A, T_B etc. of the reactants A, B etc. This means solving a set of nonlinear equations of mass-balance

${\displaystyle {\begin{aligned}{\ce {T_{A}}}&=[{\ce {A}}]+\sum _{1,nk}p\beta _{pq\cdots }[{\ce {A}}]^{p}[{\ce {B}}]^{q}\cdots \\{\ce {T_{B}}}&=[{\ce {B}}]+\sum _{1,nk}q\beta _{pq\cdots }[{\ce {A}}]^{p}[{\ce {B}}]^{q}\cdots \\etc.\end{aligned}}}$

for the free concentrations [A], [B] etc. When the pH (or equivalent e.m.f., E).is measured, the free concentration of hydrogen ions, [H], is obtained from the measured value as

${\displaystyle [\mathrm {H} ]=10^{-\mathrm {pH} }}$ or ${\displaystyle [\mathrm {H} ]=10^{\mathrm {\frac {E-E^{0}}{2.303RT/nF}} }}$

and only the free concentrations of the other reactants are calculated. The concentrations of the complexes are derived from the free concentrations via the chemical model.

Some authors^[5][6] include the free reactant terms in the sums by declaring identity (unit) β constants for which the stoichiometric coefficients are 1 for the reactant concerned and zero for all other reactants. For example, with 2 reagents, the mass-balance equations assume the simpler form.

${\displaystyle {\begin{aligned}T_{\ce {A}}&=\sum _{0,nk}p\beta _{pq}[{\ce {A}}]^{p}[{\ce {B}}]^{q}\\[4pt]T_{\ce {B}}&=\sum _{0,nk}q\beta _{pq}[{\ce {A}}]^{p}[{\ce {B}}]^{q}\\\end{aligned}}}$

${\displaystyle \beta _{10}=\beta _{01}=1}$

In this manner, all chemical species, including the free reactants, are treated in the same way, having been formed from the combination of reactants that is specified by the stoichiometric coefficients.

In a titration system the analytical concentrations of the reactants at each titration point are obtained from the initial conditions, the burette concentrations and volumes. The analytical (total) concentration of a reactant R at the ith titration point is given by

${\displaystyle T_{{\ce {R}}}={\frac {{\ce {R}}_{0}+v_{i}{\ce {[R]}}}{v_{0}+v_{i}}}}$

where R₀ is the initial amount of R in the titration vessel, v₀ is the initial volume, [R] is the concentration of R in the burette and v_i is the volume added. The burette concentration of a reactant not present in the burette is taken to be zero.

In general, solving these nonlinear equations presents a formidable challenge because of the huge range over which the free concentrations may vary. At the beginning, values for the free concentrations must be estimated. Then, these values are refined, usually by means of Newton–Raphson iterations. The logarithms of the free concentrations may be refined rather than the free concentrations themselves. Refinement of the logarithms of the free concentrations has the added advantage of automatically imposing a non-negativity constraint on the free concentrations. Once the free reactant concentrations have been calculated, the concentrations of the complexes are derived from them and the equilibrium constants.

Note that the free reactant concentrations can be regarded as implicit parameters in the equilibrium constant refinement process. In that context the values of the free concentrations are constrained by forcing the conditions of mass-balance to apply at all stages of the process.

Equilibrium constant refinement

The objective of the refinement process is to find equilibrium constant values that give the best fit to the experimental data. This is usually achieved by minimising an objective function, U, by the method of non-linear least-squares. First the residuals are defined as

${\displaystyle r_{i}=y_{i}^{\text{obs}}-y_{i}^{\text{calc}}}$

Then the most general objective function is given by

${\displaystyle U=\sum _{i}\sum _{j}r_{i}W_{ij}r_{j}\,}$

The matrix of weights, W, should be, ideally, the inverse of the variance-covariance matrix of the observations. It is rare for this to be known. However, when it is, the expectation value of U is one, which means that the data are fitted within experimental error. Most often only the diagonal elements are known, in which case the objective function simplifies to

${\displaystyle U=\sum _{i}W_{ii}r_{i}^{2}}$

with W_ij = 0 when j ≠ i. Unit weights, W_ii = 1, are often used but, in that case, the expectation value of U is the root mean square of the experimental errors.

The minimization may be performed using the Gauss–Newton method. Firstly the objective function is linearised by approximating it as a first-order Taylor series expansion about an initial parameter set, p.

${\displaystyle U=U^{0}+\sum _{i}{\frac {\partial U}{\partial p_{i}}}\delta p_{i}}$

The increments δp_i are added to the corresponding initial parameters such that U is less than U⁰. At the minimum the derivatives ∂U/∂p_i, which are simply related to the elements of the Jacobian matrix, J

${\displaystyle J_{jk}={\frac {\partial y_{j}^{\mathrm {calc} }}{\partial p_{k}}}}$

where p_k is the kth parameter of the refinement, are equal to zero. One or more equilibrium constants may be parameters of the refinement. However, the measured quantities (see above) represented by y are not expressed in terms of the equilibrium constants, but in terms of the species concentrations, which are implicit functions of these parameters. Therefore, the Jacobian elements must be obtained using implicit differentiation.

The parameter increments δp are calculated by solving the normal equations, derived from the conditions that ∂U/∂p = 0 at the minimum.

${\displaystyle {\left(J^{\mathrm {T} }WJ\right)\delta p=J^{\mathrm {T} }Wr}}$

The increments δp are added iteratively to the parameters

${\displaystyle \mathbf {p} ^{n+1}=\mathbf {p} ^{n}+\delta \mathbf {p} }$

where n is an iteration number. The species concentrations and y^calc values are recalculated at every data point. The iterations are continued until no significant reduction in U is achieved, that is, until a convergence criterion is satisfied. If, however, the updated parameters do not result in a decrease of the objective function, that is, if divergence occurs, the increment calculation must be modified. The simplest modification is to use a fraction, f, of calculated increment, so-called shift-cutting.

${\displaystyle \mathbf {p} ^{n+1}=\mathbf {p} ^{n}+f\mathbf {\delta p} }$

In this case, the direction of the shift vector, δp, is unchanged. With the more powerful Levenberg–Marquardt algorithm, on the other hand, the shift vector is rotated towards the direction of steepest descent, by modifying the normal equations,

${\displaystyle \mathbf {\left(J^{\mathrm {T} }WJ+\lambda I\right)\delta p=J^{\mathrm {T} }Wr} }$

where λ is the Marquardt parameter and I is an identity matrix. Other methods of handling divergence have been proposed.^[6]

A particular issue arises with NMR and spectrophotometric data. For the latter, the observed quantity is absorbance, A, and the Beer–Lambert law can be written as

${\displaystyle A_{\lambda }=l\sum (\varepsilon _{pq..})_{\lambda }c_{pq..}}$

It can be seen that, assuming that the concentrations, c, are known, that absorbance, A, at a given wavelength, ${\displaystyle \lambda }$, and path length ${\displaystyle l}$, is a linear function of the molar absorbptivities, ε. With 1 cm path-length, in matrix notation

${\displaystyle \mathbf {A} ={\boldsymbol {\varepsilon }}\mathbf {C} \,}$

There are two approaches to the calculation of the unknown molar absorptivities

(1) The ε values are considered parameters of the minimization and the Jacobian is constructed on that basis. However, the ε values themselves are calculated at each step of the refinement by linear least-squares:

${\displaystyle {\boldsymbol {\varepsilon }}=\mathbf {\left(C^{\mathrm {T} }C\right)^{-1}C^{\mathrm {T} }A} }$

using the refined values of the equilibrium constants to obtain the speciation. The matrix

${\displaystyle \mathbf {\left(C^{T}C\right)^{-1}C^{T}} }$

is an example of a pseudo-inverse.

Golub and Pereyra^[7] showed how the pseudo-inverse can be differentiated so that parameter increments for both molar absorptivities and equilibrium constants can be calculated by solving the normal equations.

(2) The Beer–Lambert law is written as

${\displaystyle \mathbf {\boldsymbol {\varepsilon }} _{\lambda }=\mathbf {A} _{\lambda }^{-1}\mathbf {C} \,}$

The unknown molar absorbances of all "coloured" species are found by using the non-iterative method of linear least-squares, one wavelength at a time. The calculations are performed once every refinement cycle, using the stability constant values obtaining at that refinement cycle to calculate species' concentration values in the matrix ${\displaystyle \mathbf {C} }$.

Model selection

Once a refinement has been completed the results should be checked to verify that the chosen model is acceptable. generally speaking, a model is acceptable when the data are fitted within experimental error, but there is no single criterion to use to make the judgement. The following should be considered.

Distribution of residuals

At the minimum in U the system can be approximated to a linear one, the residuals in the case of unit weights are related to the observations by

${\displaystyle \mathbf {r=y^{\mathrm {obs} }-J\left(J^{\mathrm {T} }T\right)^{-1}J^{\mathrm {T} }y^{\mathrm {obs} }} }$

The symmetric, idempotent matrix J(J^TT)⁻¹J is known in the statistics literature as the hat matrix, H. Thus,

${\displaystyle \mathbf {r=\left(I-H\right)y^{\mathrm {obs} }} }$

and

${\displaystyle \mathbf {M^{r}=\left(I-H\right)M^{y}\left(I-H\right)} }$

where I is an identity matrix and M^r and M^y are the variance-covariance matrices of the residuals and observations, respectively. This shows that even though the observations may be uncorrelated, the residuals are always correlated.

The diagram at the right shows the result of a refinement of the stability constants of Ni(Gly)⁺, Ni(Gly)₂ and Ni(Gly)⁻
₃ (where GlyH = glycine). The observed values are shown a blue diamonds and the species concentrations, as a percentage of the total nickel, are superimposed. The residuals are shown in the lower box. The residuals are not distributed as randomly as would be expected. This is due to the variation of liquid junction potentials and other effects at the glass/liquid interfaces. Those effects are very slow compared to the rate at which equilibrium is established.