# A tutorial on potentiometric data processing. Analysis of software for optimization of protonation constants

Silvia Berto, Salvador Blasco, Lorenzo Castellino, Aleksandar Cvetkovski, Concetta de Stefano, Sofia Gama, Enrique García-España, Petr Hermann, Gabriele Lando, Matteo Marafante, et al.

► To cite this version:

Silvia Berto, Salvador Blasco, Lorenzo Castellino, Aleksandar Cvetkovski, Concetta de Stefano, et al. A tutorial on potentiometric data processing. Analysis of software for optimization of protonation constants. Analytica Chimica Acta, 2024, 1303, pp.342476. 10.1016/j.aca.2024.342476. hal-04544012

HAL Id: hal-04544012

https://hal.science/hal-04544012

Submitted on 12 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A tutorial on potentiometric data processing. Analysis of software for optimization of protonation constants

Silvia Berto, Salvador Blasco, Lorenzo Castellino, Aleksandar Cvetkovski, Concetta de Stefano, Sofia Gama, Enrique García-España, Petr Hermann, Gabriele Lando, Matteo Marafante, et al.

► To cite this version:

Silvia Berto, Salvador Blasco, Lorenzo Castellino, Aleksandar Cvetkovski, Concetta de Stefano, et al. A tutorial on potentiometric data processing. Analysis of software for optimization of protonation constants. Analytica Chimica Acta, 2024, 1303, pp.342476. 10.1016/j.aca.2024.342476. hal-04544012

HAL Id: hal-04544012

https://hal.science/hal-04544012v1

Submitted on 12 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A tutorial on potentiometric data processing. Analysis of software for optimization of protonation constants

Silvia Berto<sup>a, *</sup>, Salvador Blasco<sup>b</sup>, Lorenzo Castellino<sup>a</sup>, Aleksandar Cvetkovski<sup>c</sup>, Concetta De Stefano<sup>d</sup>, Sofia Gama<sup>e</sup>, Enrique García-España<sup>b</sup>, Petr Hermann<sup>f</sup>, Gabriele Lando<sup>d</sup>, Matteo Marafante<sup>a</sup>, Michel Meyer<sup>g</sup>, Winfried Plass<sup>h</sup>, Lauryn Quinodoz<sup>g,i</sup>, Demetrio Milea<sup>d</sup>

## HIGHLIGHTS

* Five software used to optimize logK from potentiometric data are critically evaluated.

* A simulated titration dataset of a hexaprotic acid was processed by five software.

* The impact on logK of systematic errors occurring during titrations is evaluated.

* The impact on logK of ionic strength variations occurring during titrations is evaluated.

* Guidelines for data acquisition and treatment are given.

## GRAPHICAL ABSTRACT

<figure><img src="image_2.png" /><figcaption>Analysis of tools for optimization of equilibrium constants from potentiometric data</figcaption></figure>

## ARTICLE INFO

**Keywords:**

Equilibrium constants

Chemical thermodynamics

Potentiometry

Data analysis

## ABSTRACT

Defining the distribution of the chemical species in a multicomponent system is a task of great importance with applications in many fields. To clarify the identity and the abundance of the species that can be formed by the interaction of the components of a solution, it is fundamental to know the formation constants of those species. The determination of equilibrium constants is mainly performed through the analysis of experimental data obtained by different instrumental techniques. Among them, potentiometry is the elective technique for this

<footnote><sup>a</sup> Università di Torino, Dipartimento di Chimica, Via P. Giuria 7, 10125, Torino, Italy</footnote>

<footnote><sup>b</sup> Institute of Molecular Sciences, Departamento de Química Inorgánica, University of Valencia, C/Catedrático José Beltrán Martínez 2, 46980, Paterna, Valencia, Spain</footnote>

<footnote><sup>c</sup> Faculty of Medical Sciences, Goce Delcev Universit, Krste Misirkov Bb, PO 201, 2000, Štip, Macedonia</footnote>

<footnote><sup>d</sup> Università degli Studi di Messina, Dipartimento di Scienze Chimiche, Biologiche, Farmaceutiche ed Ambientali, CHIBIOFARAM, V.le F. Stagno d'Alcontres, 31, 98166, Messina, Italy</footnote>

<footnote><sup>e</sup> Centro de Ciências e Tecnologias Nucleares, Instituto Superior Técnico, Universidade de Lisboa, Estrada Nacional 10 (km 139.7), 2695-066, Bobadela LRS, Portugal</footnote>

<footnote><sup>f</sup> Department of Inorganic Chemistry, Faculty of Science, Charles University, Hlavova 8, 128 40 Prague 2, Czech Republic</footnote>

<footnote><sup>g</sup> Institut de Chimie Moléculaire de l'Université de Bourgogne (ICMUB), UMR 6302, CNRS, Université de Bourgogne, 9 avenue Alain Savary, BP 47870, 21078, Dijon, Cedex, France</footnote>

<footnote><sup>h</sup> Institut für Anorganische und Analytische Chemie, Friedrich Schiller University Jena, Humboldtstr 8, D-07743, Jena, Germany</footnote>

<footnote><sup>i</sup> University of Applied Sciences Western Switzerland HES-SO, Haute Ecole d'Ingénierie et d'Architecture, Pérolles 80, 1700, Fribourg, Switzerland</footnote>

<footnote>* Corresponding author.
E-mail address: silvia.berto@unito.it (S. Berto).</footnote>
# A tutorial on potentiometric data processing. Analysis of software for optimization of protonation constants

Silvia Berto<sup>a, *</sup>, Salvador Blasco<sup>b</sup>, Lorenzo Castellino<sup>a</sup>, Aleksandar Cvetkovski<sup>c</sup>, Concetta De Stefano<sup>d</sup>, Sofia Gama<sup>e</sup>, Enrique García-España<sup>b</sup>, Petr Hermann<sup>f</sup>, Gabriele Lando<sup>d</sup>, Matteo Marafante<sup>a</sup>, Michel Meyer<sup>g</sup>, Winfried Plass<sup>h</sup>, Lauryn Quinodoz<sup>g,i</sup>, Demetrio Milea<sup>d</sup>

## HIGHLIGHTS

* Five software used to optimize logK from potentiometric data are critically evaluated.

* A simulated titration dataset of a hexaprotic acid was processed by five software.

* The impact on logK of systematic errors occurring during titrations is evaluated.

* The impact on logK of ionic strength variations occurring during titrations is evaluated.

* Guidelines for data acquisition and treatment are given.

## GRAPHICAL ABSTRACT

### Analysis of tools for optimization of equilibrium constants from potentiometric data

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

**Keywords:**

Equilibrium constants

Chemical thermodynamics

Potentiometry

Data analysis

## ABSTRACT

Defining the distribution of the chemical species in a multicomponent system is a task of great importance with applications in many fields. To clarify the identity and the abundance of the species that can be formed by the interaction of the components of a solution, it is fundamental to know the formation constants of those species. The determination of equilibrium constants is mainly performed through the analysis of experimental data obtained by different instrumental techniques. Among them, potentiometry is the elective technique for this

<footnote><sup>a</sup> Università di Torino, Dipartimento di Chimica, Via P. Giuria 7, 10125, Torino, Italy</footnote>

<footnote><sup>b</sup> Institute of Molecular Sciences, Departamento de Química Inorgánica, University of Valencia, C/Catedrático José Beltrán Martínez 2, 46980, Paterna, Valencia, Spain</footnote>

<footnote><sup>c</sup> Faculty of Medical Sciences, Goce Delcev Universit, Krste Misirkov Bb, PO 201, 2000, Štip, Macedonia</footnote>

<footnote><sup>d</sup> Università degli Studi di Messina, Dipartimento di Scienze Chimiche, Biologiche, Farmaceutiche ed Ambientali, CHIBIOFARAM, V.le F. Stagno d'Alcontres, 31, 98166, Messina, Italy</footnote>

<footnote><sup>e</sup> Centro de Ciências e Tecnologias Nucleares, Instituto Superior Técnico, Universidade de Lisboa, Estrada Nacional 10 (km 139.7), 2695-066, Bobadela LRS, Portugal</footnote>

<footnote><sup>f</sup> Department of Inorganic Chemistry, Faculty of Science, Charles University, Hlavova 8, 128 40 Prague 2, Czech Republic</footnote>

<footnote><sup>g</sup> Institut de Chimie Moléculaire de l'Université de Bourgogne (ICMUB), UMR 6302, CNRS, Université de Bourgogne, 9 avenue Alain Savary, BP 47870, 21078, Dijon, Cedex, France</footnote>

<footnote><sup>h</sup> Institut für Anorganische und Analytische Chemie, Friedrich Schiller University Jena, Humboldtstr 8, D-07743, Jena, Germany</footnote>

<footnote><sup>i</sup> University of Applied Sciences Western Switzerland HES-SO, Haute Ecole d'Ingénierie et d'Architecture, Pérolles 80, 1700, Fribourg, Switzerland</footnote>

<footnote>* Corresponding author.
E-mail address: silvia.berto@unito.it (S. Berto).</footnote>
purpose. As such, a survey was run within the NECTAR COST Action – Network for Equilibria and Chemical Thermodynamics Advanced Research, to identify the most used software for the analysis of potentiometric data and to highlight their strengths and weaknesses. The features and the calculation processes of each software were analyzed and rationalized, and a simulated titration dataset of a hypothetic hexaprotic acid was processed by each software to compare and discuss the optimized protonation constants. Moreover, further data analysis was also carried out on the original dataset including some systematic errors from different sources, as some calibration parameters, the total analytical concentration of reagents and ionic strength variations during titrations, to evaluate their impact on the refined parameters. Results showed that differences on the protonation constants estimated by the tested software are not significant, while some of the considered systematic errors affect results. Overall, it emerged that software commonly used suffer from many limitations, highlighting the urgency of new dedicated and modern tools. In this context, some guidelines for data generation and treatment are also given.

# 1. Introduction

Chemical speciation studies [1] are of great importance in many fields: (i) geochemists and environmental chemists are interested in studying the real composition of all types of natural waters; (ii) biochemists are interested in the speciation in fluids such as blood, plasma, urine or saliva; (iii) for industrial chemists, it is important to define the species formed in formulations (detergents, personal care products, drugs, lubricants, etc.); (iv) chemists, in general, are often keen in the knowledge of the composition of multiphasic/multicomponent systems for kinetic, synthetic, or other kinds of studies. To estimate the abundance of complexes or adducts formed at equilibrium in solutions under well-defined conditions (e.g. solvent, pH, concentrations, temperature, pressure, ionic strength, salinity, etc.), it is fundamental to know the identity of the various species and the associated thermodynamic parameters of the chemical equilibria describing their formation. If this information is not available or the data are unreliable, the determination of the stoichiometric formation constant for each species prevailing in a given medium at constant temperature and pressure is a prerequisite. In thermodynamic terms, a chemical equilibrium analysis can be carried out by solving a set of nonlinear equations consisting of mass balance equations in which the unknowns to determine are commonly the concentrations of the formed species and/or their stability constants.

In this light, a significant boost to the study of chemical equilibria has been given by the development of Information Technology (IT) products. In 1961, the first general computer program for the evaluation of equilibrium constants from equilibrium data was announced [2]. From that date onwards, many research groups have developed their own software independently, with some being sold as a commercial product. In the 1980s, several reviews were written about computer programs thought for the evaluation of formation constants [3–5]. The main features for judging the advantages of applicable methods relate to their robustness in convergence for many different systems, their flexibility and effortlessness in designing a model to describe a certain chemical process, and the fast achievement of quadratic convergence [6].

To solve the system of mass balance equations, most of the software, such as LETAGROP [7], MICROQL [8], SUPERQUAD [9], ESTA [10], IMPACT [11], CHESS [12], BSTAC [13], HYPERQUAD [14], or PHREEQC [15], use the iterative Newton-Raphson method or some of its modifications. The Newton-Raphson method is subject to non-convergence when inaccurate initial conditions are selected and it is often implemented by techniques that apply a preliminary selection of initial solutions [12,16,17]. Other computer programs have been developed in the successive decades and, for those that are still in use, the code of the current versions has been updated or rewritten with different programming languages. In any case, the mathematical approach has not changed much in the last 20 years. Unfortunately, most of these software applications are no longer in an active development state and several of them have already lost compatibility with current operating systems.

Recently, under the frame of the aims of the COST Action CA18202 NECTAR – Network for Equilibria and Chemical Thermodynamics Advanced Research, a survey among several CA18202 NECTAR

participants was conducted regarding the most used software for potentiometric data analysis and the problems that users are facing. The outcomes show a high fragmentation in the software used, associated with a general dissatisfaction concerning the user experience. This highlights the need of a critical evaluation of the most used software for the analysis of potentiometric data, in order to identify the strengths and weaknesses of each. The obtained knowledge could be useful for the future development of a new generation of Information Technology (IT) products.

Potentiometry is the main technique used to measure equilibrium constants for metal/ligand/proton systems in solution [18] or the stability constants of metal-ionophore complexes in the polymeric sensing membrane of ionophore-based ion selective electrodes [19,20]. This work aims to treat only the first application and the most used software by the NECTAR community, HYPERQUAD [14], SUPERQUAD [21], BSTAC [13], ReactLab™ pH PRO [22], and KEV [23], have been considered. In a first attempt to evaluate their performances, we processed with each aforementioned refinement programs an artificial dataset composed of six different titration curves simulated by the PyES program [24] for a hypothetical hexaprotic acid, with the aim of optimizing the six protonation constants. We have deliberately chosen to work with simulated data in order to rule out any possible source of uncertainty originating from experimental procedures and to attribute directly any bias in the refined protonation constants to the algorithm implemented in the different software. Moreover, the intercomparison test was performed on a chemical system for which the species involved are known, therefore the definition of a chemical model is not required. The features of each software are outlined and the results obtained, i.e. protonation constants and related uncertainty values, have also been analyzed and discussed. Furthermore, the same dataset but spoiled by different kinds of systematic errors, which could occur by carelessness during the process of data entry or by flaws in the experimental procedures, have also been processed in order to evaluate the impact of these errors on the refined parameters. The errors considered were: the use of (i) an incorrect formal electrode potential, (ii) a wrong titrant concentration, (iii) a wrong concentration of the solution components, and (iv) changes on the ionic strength during the titrations. The discussion of the results was also used to highlight the extent of the effects of errors on the refined parameters and to propose guidelines for the definition of reliable formation constants.

# 2. Software and dataset

## 2.1. Software used for the data analysis

The BSTAC [13] is a freeware written in BASIC and working in DOS environment. Recently the source code and the executable files were made available at https://github.com/Kastakin/BSTAC, though the graphic output only works partially due to incompatibility with modern video mode. BSTAC uses Newton-Raphson method for root finding, and the Marquardt-Levenberg algorithm for nonlinear least squares minimization, allowing the refinement of all the analytical parameters from one, or several, potentiometric titrations at constant or variable ionic
purpose. As such, a survey was run within the NECTAR COST Action – Network for Equilibria and Chemical Thermodynamics Advanced Research, to identify the most used software for the analysis of potentiometric data and to highlight their strengths and weaknesses. The features and the calculation processes of each software were analyzed and rationalized, and a simulated titration dataset of a hypothetic hexaprotic acid was processed by each software to compare and discuss the optimized protonation constants. Moreover, further data analysis was also carried out on the original dataset including some systematic errors from different sources, as some calibration parameters, the total analytical concentration of reagents and ionic strength variations during titrations, to evaluate their impact on the refined parameters. Results showed that differences on the protonation constants estimated by the tested software are not significant, while some of the considered systematic errors affect results. Overall, it emerged that software commonly used suffer from many limitations, highlighting the urgency of new dedicated and modern tools. In this context, some guidelines for data generation and treatment are also given.

# 1. Introduction

Chemical speciation studies [1] are of great importance in many fields: (i) geochemists and environmental chemists are interested in studying the real composition of all types of natural waters; (ii) biochemists are interested in the speciation in fluids such as blood, plasma, urine or saliva; (iii) for industrial chemists, it is important to define the species formed in formulations (detergents, personal care products, drugs, lubricants, etc.); (iv) chemists, in general, are often keen in the knowledge of the composition of multiphasic/multicomponent systems for kinetic, synthetic, or other kinds of studies. To estimate the abundance of complexes or adducts formed at equilibrium in solutions under well-defined conditions (e.g. solvent, pH, concentrations, temperature, pressure, ionic strength, salinity, etc.), it is fundamental to know the identity of the various species and the associated thermodynamic parameters of the chemical equilibria describing their formation. If this information is not available or the data are unreliable, the determination of the stoichiometric formation constant for each species prevailing in a given medium at constant temperature and pressure is a prerequisite. In thermodynamic terms, a chemical equilibrium analysis can be carried out by solving a set of nonlinear equations consisting of mass balance equations in which the unknowns to determine are commonly the concentrations of the formed species and/or their stability constants.

In this light, a significant boost to the study of chemical equilibria has been given by the development of Information Technology (IT) products. In 1961, the first general computer program for the evaluation of equilibrium constants from equilibrium data was announced [2]. From that date onwards, many research groups have developed their own software independently, with some being sold as a commercial product. In the 1980s, several reviews were written about computer programs thought for the evaluation of formation constants [3–5]. The main features for judging the advantages of applicable methods relate to their robustness in convergence for many different systems, their flexibility and effortlessness in designing a model to describe a certain chemical process, and the fast achievement of quadratic convergence [6].

To solve the system of mass balance equations, most of the software, such as LETAGROP [7], MICROQL [8], SUPERQUAD [9], ESTA [10], IMPACT [11], CHESS [12], BSTAC [13], HYPERQUAD [14], or PHREEQC [15], use the iterative Newton-Raphson method or some of its modifications. The Newton-Raphson method is subject to non-convergence when inaccurate initial conditions are selected and it is often implemented by techniques that apply a preliminary selection of initial solutions [12,16,17]. Other computer programs have been developed in the successive decades and, for those that are still in use, the code of the current versions has been updated or rewritten with different programming languages. In any case, the mathematical approach has not changed much in the last 20 years. Unfortunately, most of these software applications are no longer in an active development state and several of them have already lost compatibility with current operating systems.

Recently, under the frame of the aims of the COST Action CA18202 NECTAR – Network for Equilibria and Chemical Thermodynamics Advanced Research, a survey among several CA18202 NECTAR

participants was conducted regarding the most used software for potentiometric data analysis and the problems that users are facing. The outcomes show a high fragmentation in the software used, associated with a general dissatisfaction concerning the user experience. This highlights the need of a critical evaluation of the most used software for the analysis of potentiometric data, in order to identify the strengths and weaknesses of each. The obtained knowledge could be useful for the future development of a new generation of Information Technology (IT) products.

Potentiometry is the main technique used to measure equilibrium constants for metal/ligand/proton systems in solution [18] or the stability constants of metal-ionophore complexes in the polymeric sensing membrane of ionophore-based ion selective electrodes [19,20]. This work aims to treat only the first application and the most used software by the NECTAR community, HYPERQUAD [14], SUPERQUAD [21], BSTAC [13], ReactLab™ pH PRO [22], and KEV [23], have been considered. In a first attempt to evaluate their performances, we processed with each aforementioned refinement programs an artificial dataset composed of six different titration curves simulated by the PyES program [24] for a hypothetical hexaprotic acid, with the aim of optimizing the six protonation constants. We have deliberately chosen to work with simulated data in order to rule out any possible source of uncertainty originating from experimental procedures and to attribute directly any bias in the refined protonation constants to the algorithm implemented in the different software. Moreover, the intercomparison test was performed on a chemical system for which the species involved are known, therefore the definition of a chemical model is not required. The features of each software are outlined and the results obtained, i.e. protonation constants and related uncertainty values, have also been analyzed and discussed. Furthermore, the same dataset but spoiled by different kinds of systematic errors, which could occur by carelessness during the process of data entry or by flaws in the experimental procedures, have also been processed in order to evaluate the impact of these errors on the refined parameters. The errors considered were: the use of (i) an incorrect formal electrode potential, (ii) a wrong titrant concentration, (iii) a wrong concentration of the solution components, and (iv) changes on the ionic strength during the titrations. The discussion of the results was also used to highlight the extent of the effects of errors on the refined parameters and to propose guidelines for the definition of reliable formation constants.

# 2. Software and dataset

## 2.1. Software used for the data analysis

The BSTAC [13] is a freeware written in BASIC and working in DOS environment. Recently the source code and the executable files were made available at https://github.com/Kastakin/BSTAC, though the graphic output only works partially due to incompatibility with modern video mode. BSTAC uses Newton-Raphson method for root finding, and the Marquardt-Levenberg algorithm for nonlinear least squares minimization, allowing the refinement of all the analytical parameters from one, or several, potentiometric titrations at constant or variable ionic
strength, also within a single titration. BSTAC allows the refinement of different types of equilibrium constants (protonation, hydrolysis, complexation), together with all titration parameters (formal electrode potential, slope of the electrode, junction potential coefficients in the highly acidic and alkaline ranges), as well as the analytical concentration of all components. The software has also the capability to take into account the variation of the ionic strength occurring during a titration by refining the parameters expressing the ionic strength dependence of equilibrium constants and the electrode formal potential according to an Extended Debye–Hückel equation [13] (more details are reported in paragraph 2.3. Ionic strength dependence). Furthermore, it also allows the constrained refinement of analytical concentrations. The constrain between the component concentrations forces the system to vary their values, while maintaining their ratio constant during the calculations. This feature proves particularly useful for checking the concentration of stock solutions and their purities. Lastly, BSTAC can also deal with titrations performed by any kind of ion selective electrodes (ISE).

SUPERQUAD [21] is a FORTRAN program designed for fitting equilibrium constants based on potentiometric titration data. It is an improvement of its predecessor MINIQUAD [25], developed by the same authors. It is a text in/text out program in which a text file containing the model and the titration data is given as input and a text output is produced. SUPERQUAD allows the fitting of equilibrium constants as well as “dangerous parameters”, such as titrant concentration, initial amounts of reagents, and standard potential, although the use of this feature is discouraged. The fitting parameters are refined by means of the Marquardt–Levenberg algorithm. It assumes a Nernstian behavior of the electrode system, but ionic strength correction is not implemented. The weighting scheme is based on the inverse of the variance of the potential and volume calculated with the error propagation formula. Once a convergence criterion is met, if all the fitting parameters are positive and within the specified error threshold, the program stops. If not, the worst constant is removed and the fitting starts over with the reduced model.

HYPERQUAD [14] is the improved version of SUPERQUAD developed by the same authors. It is a program able to determine equilibrium constants from potentiometric data, and it is part of the homonym commercial suite [26] of various software for the data analysis from different experimental techniques (e.g. spectrophotometry, NMR, calorimetry). Over the years, different versions of HYPERQUAD have been released with some important differences among them. In particular, the HYPERQUAD2008 version refines the equilibrium constants and not their logarithmic counterparts as was originally done by SUPERQUAD, and even by BSTAC, but this functionality was reintroduced in the later HYPERQUAD2013 version. The change was based on the observation that using the logarithmic values, a natural non-negativity constraint is imposed on the equilibrium constants, improving at the same time the stability of the refinement procedure [27]. It is specified that this is not the only way to enforce a positive sign of the parameter, but it is difficult to envision the possibility of a “zero-cost” solution to the problem, e.g., some drawbacks will always be present, such as computational cost and/or physical relevance of obtained solutions. On the topic of protecting the calculation process from refinement discrepancies, HYPERQUAD uses the shift-cutting method, whereas SUPERQUAD uses the Marquardt method, once again the authors justify this change being a more generally applicable methodology. Last but not least, unlike SUPERQUAD, HYPERQUAD is presented through a graphical interface and it is generally compatible with 64 bit hardware.

KEV [23] is an online optimizer of the equilibrium constants. The data analysis applies the iterative Newton algorithm and the Hooke–Jeeves method. The web application requires the simple authentication, signing up is free but a letter has to be sent to the authors. KEV is able to treat potentiometric, spectrophotometric, NMR, and calorimetric data. By selecting the EMF tab, unknown equilibrium constants are derived from experimental potentiometric data, entered either manually or by uploading a text file having the comma separated values (csv)

format. The output file can be downloaded locally as a csv or as a Microsoft Excel® (in xlsx format) file. The input data are the stoichiometric indexes of reactions, the decimal logarithms of equilibrium constants (initial guesses for refined constants or best estimates for the already known constants treated as fixed values), the total analytical concentrations of reagents, the emf (electromotive force) values together with their experimental standard deviations obtained during the titration process, and the electrode calibration parameters, i.e. formal electrode potential and Nernstian slope. The volume of titrant added is not directly entered, as this information is extracted from the table of total concentrations for each component provided by the operator, each line of this table corresponding to a single titration point. This means that a pre-treatment of experimental data is required before uploading them in the calculation process. Moreover, it is not possible to analyze more than one titration within a single run.

ReactLab™ pH PRO [22] is a software launched in January 2020. In that respect, it is the most recent IT product evaluated in this work and the demo version of the software was used for data processing. ReactLab™ pH PRO can refine the concentration of reagents in solution and either global or stepwise equilibrium constants for single and multi-titration experiments, which can be recorded under different conditions. In the latter case, the program also offers the possibility to define and optimize so-called auxiliary parameters (e.g. the concentration of the stock solution of a reagent used in several titrations at different dilution levels, the volume of the stock solution used for each experiment being treated as a constant parameter). Both continuous (also known as automatic, i.e. steady increase of the total volume) or batch titration data (constant total volume) can be handled. The program has been developed in the Matlab® environment and compiled to run on Windows as a standalone application. Data entry is handled through a custom macro-powered Microsoft Excel® worksheet. Based on an adaptation of the Marquardt–Levenberg algorithm, the data-fitting procedure minimizes the unweighted sum of squared residuals between calculated and experimental emf or pH readings. For the second option, the pH-electrode calibration parameters (formal electrode potential and slope) can be also refined internally, according to the Nernst equation. An additional term, correcting junction potential errors in the acidic range only, can be entered by the user but cannot be refined. ReactLab™ pH PRO provides both a graphical and a textual output that includes standard deviations for each fitted parameter, as well as the total sum-of-squares. One peculiarity of the package is that the data analysis can also be based on activities rather than concentrations. In that instance, activity coefficients are approximated by the Debye–Hückel, extended Debye–Hückel, Davies, or Specific ion Inter-action Theory equations. As the mass action law is solved in terms of species activities rather than concentrations, the program returns the so-called thermodynamic equilibrium constants (constants at infinite dilution or zero ionic strength) instead of stoichiometric equilibrium constants, which are typically derived from titrations carried out in the presence of a supporting electrolyte taken in large excess over all other reagents. Considering that ionic strength variations during such experiments cannot be strictly avoided by the addition of a background salt that moreover might interfere with the investigated chemical system and/or the measuring probe (cf. Interferences of Li⁺, Na⁺, or K⁺ with the glass electrode in the alkaline solutions), the authors sought to introduce this calculation feature for processing potentiometric data collected in the absence of a supporting electrolyte.

## 2.2. Error minimization

BSTAC, HYPERQUAD, KEV, and SUPERQUAD refine the equilibrium constants by a weighted least squares method. The algorithm minimizes the weighted residual square sum, which is a measure of the difference between a measured observable and its approximation as predicted by the model. The observable can be the potential reading or the added volumes of titrant (directly related with the analytical concentration of
strength, also within a single titration. BSTAC allows the refinement of different types of equilibrium constants (protonation, hydrolysis, complexation), together with all titration parameters (formal electrode potential, slope of the electrode, junction potential coefficients in the highly acidic and alkaline ranges), as well as the analytical concentration of all components. The software has also the capability to take into account the variation of the ionic strength occurring during a titration by refining the parameters expressing the ionic strength dependence of equilibrium constants and the electrode formal potential according to an Extended Debye–Hückel equation [13] (more details are reported in paragraph 2.3. Ionic strength dependence). Furthermore, it also allows the constrained refinement of analytical concentrations. The constrain between the component concentrations forces the system to vary their values, while maintaining their ratio constant during the calculations. This feature proves particularly useful for checking the concentration of stock solutions and their purities. Lastly, BSTAC can also deal with titrations performed by any kind of ion selective electrodes (ISE).

SUPERQUAD [21] is a FORTRAN program designed for fitting equilibrium constants based on potentiometric titration data. It is an improvement of its predecessor MINIQUAD [25], developed by the same authors. It is a text in/text out program in which a text file containing the model and the titration data is given as input and a text output is produced. SUPERQUAD allows the fitting of equilibrium constants as well as “dangerous parameters”, such as titrant concentration, initial amounts of reagents, and standard potential, although the use of this feature is discouraged. The fitting parameters are refined by means of the Marquardt–Levenberg algorithm. It assumes a Nernstian behavior of the electrode system, but ionic strength correction is not implemented. The weighting scheme is based on the inverse of the variance of the potential and volume calculated with the error propagation formula. Once a convergence criterion is met, if all the fitting parameters are positive and within the specified error threshold, the program stops. If not, the worst constant is removed and the fitting starts over with the reduced model.

HYPERQUAD [14] is the improved version of SUPERQUAD developed by the same authors. It is a program able to determine equilibrium constants from potentiometric data, and it is part of the homonym commercial suite [26] of various software for the data analysis from different experimental techniques (e.g. spectrophotometry, NMR, calorimetry). Over the years, different versions of HYPERQUAD have been released with some important differences among them. In particular, the HYPERQUAD2008 version refines the equilibrium constants and not their logarithmic counterparts as was originally done by SUPERQUAD, and even by BSTAC, but this functionality was reintroduced in the later HYPERQUAD2013 version. The change was based on the observation that using the logarithmic values, a natural non-negativity constraint is imposed on the equilibrium constants, improving at the same time the stability of the refinement procedure [27]. It is specified that this is not the only way to enforce a positive sign of the parameter, but it is difficult to envision the possibility of a “zero-cost” solution to the problem, e.g., some drawbacks will always be present, such as computational cost and/or physical relevance of obtained solutions. On the topic of protecting the calculation process from refinement discrepancies, HYPERQUAD uses the shift-cutting method, whereas SUPERQUAD uses the Marquardt method, once again the authors justify this change being a more generally applicable methodology. Last but not least, unlike SUPERQUAD, HYPERQUAD is presented through a graphical interface and it is generally compatible with 64 bit hardware.

KEV [23] is an online optimizer of the equilibrium constants. The data analysis applies the iterative Newton algorithm and the Hooke–Jeeves method. The web application requires the simple authentication, signing up is free but a letter has to be sent to the authors. KEV is able to treat potentiometric, spectrophotometric, NMR, and calorimetric data. By selecting the EMF tab, unknown equilibrium constants are derived from experimental potentiometric data, entered either manually or by uploading a text file having the comma separated values (csv)

format. The output file can be downloaded locally as a csv or as a Microsoft Excel® (in xlxs format) file. The input data are the stoichiometric indexes of reactions, the decimal logarithms of equilibrium constants (initial guesses for refined constants or best estimates for the already known constants treated as fixed values), the total analytical concentrations of reagents, the emf (electromotive force) values together with their experimental standard deviations obtained during the titration process, and the electrode calibration parameters, i.e. formal electrode potential and Nernstian slope. The volume of titrant added is not directly entered, as this information is extracted from the table of total concentrations for each component provided by the operator, each line of this table corresponding to a single titration point. This means that a pre-treatment of experimental data is required before uploading them in the calculation process. Moreover, it is not possible to analyze more than one titration within a single run.

ReactLab™ pH PRO [22] is a software launched in January 2020. In that respect, it is the most recent IT product evaluated in this work and the demo version of the software was used for data processing. ReactLab™ pH PRO can refine the concentration of reagents in solution and either global or stepwise equilibrium constants for single and multi-titration experiments, which can be recorded under different conditions. In the latter case, the program also offers the possibility to define and optimize so-called auxiliary parameters (e.g. the concentration of the stock solution of a reagent used in several titrations at different dilution levels, the volume of the stock solution used for each experiment being treated as a constant parameter). Both continuous (also known as automatic, i.e. steady increase of the total volume) or batch titration data (constant total volume) can be handled. The program has been developed in the Matlab® environment and compiled to run on Windows as a standalone application. Data entry is handled through a custom macro-powered Microsoft Excel® worksheet. Based on an adaptation of the Marquardt–Levenberg algorithm, the data-fitting procedure minimizes the unweighted sum of squared residuals between calculated and experimental emf or pH readings. For the second option, the pH-electrode calibration parameters (formal electrode potential and slope) can be also refined internally, according to the Nernst equation. An additional term, correcting junction potential errors in the acidic range only, can be entered by the user but cannot be refined. ReactLab™ pH PRO provides both a graphical and a textual output that includes standard deviations for each fitted parameter, as well as the total sum-of-squares. One peculiarity of the package is that the data analysis can also be based on activities rather than concentrations. In that instance, activity coefficients are approximated by the Debye–Hückel, extended Debye–Hückel, Davies, or Specific ion Inter-action Theory equations. As the mass action law is solved in terms of species activities rather than concentrations, the program returns the so-called thermodynamic equilibrium constants (constants at infinite dilution or zero ionic strength) instead of stoichiometric equilibrium constants, which are typically derived from titrations carried out in the presence of a supporting electrolyte taken in large excess over all other reagents. Considering that ionic strength variations during such experiments cannot be strictly avoided by the addition of a background salt that moreover might interfere with the investigated chemical system and/or the measuring probe (cf. Interferences of Li⁺, Na⁺, or K⁺ with the glass electrode in the alkaline solutions), the authors sought to introduce this calculation feature for processing potentiometric data collected in the absence of a supporting electrolyte.

## 2.2. Error minimization

BSTAC, HYPERQUAD, KEV, and SUPERQUAD refine the equilibrium constants by a weighted least squares method. The algorithm minimizes the weighted residual square sum, which is a measure of the difference between a measured observable and its approximation as predicted by the model. The observable can be the potential reading or the added volumes of titrant (directly related with the analytical concentration of
the titrant). The effect of the two approaches on the error sum was widely discussed in Ref. [5]. Two examples of the minimized *U* function are given by Eqs. (1) and (2) [9,14,28].

$$U = \sqrt{\frac{\sum w_n (V_n^{\text{exp}} - V_n^{\text{calc}})^2}{N - p}} \quad (1)$$

or

$$U = \sqrt{\frac{\sum w_n (E_n^{\text{exp}} - E_n^{\text{calc}})^2}{N - p}} \quad (2)$$

In both equations, $V_n^{\text{exp}}$ and $V_n^{\text{calc}}$, and $E_n^{\text{exp}}$ and $E_n^{\text{calc}}$ are the experimental and calculated volumes and potentials readings, respectively, corresponding to each titration point; $N-p$ is the number of degrees of freedom calculated as the difference between the number of experimental points ($N$) and the number of refined parameters ($p$), less the eventual constrains; $w_n$ is the weight associated to each titration point as defined in Eq. (3).

$$w_n = \frac{1}{s^2} \quad (3)$$

with

$$s^2 = s_V^2 + \left(\frac{\partial V}{\partial E}\right)^2 \times s_E^2 \quad (4a)$$

or

$$s^2 = s_E^2 + \left(\frac{\partial E}{\partial V}\right)^2 \times s_V^2 \quad (4b)$$

$s^2$ is the calculated variance of each titration point; $s_E^2$ and $s_V^2$ are the estimated variances of the potential readings and added volume, respectively, taken individually; and $\frac{\partial E}{\partial V}$ is the slope of the titration curve, calculated numerically for each pair of titration points. The weight assigned to each titration point is therefore inversely proportional to the associated variance. Therefore, the data near an equivalence point commonly affected by a high uncertainty and thus of low significance, have a large $\frac{\partial E}{\partial V}$ value and have a lower weight in the calculations than the other data in a buffered pH region. The estimated *U* value corresponds to the uncertainty of the fit and defines the uncertainty of the refined parameters through the error propagation law. Consequently, the calculation process used for refining the parameters, along with the variances entered as known values (usually $s_E^2$ and $s_V^2$), affect the uncertainty of the refined parameters. An unweighted refinement can also be as easily performed by setting the first constant term of Eq. (4) to one and the second term to zero (or to a value equal or lower to $10^{-10}$ in the case of HYPERQUAD).

The tested software manage the weight assignment to each titration point differently. BSTAC allows working with both weighted and unweighted refinement. As an added peculiarity, BSTAC offers also the possibility to execute two refinement cycles. After the first one, in which assigned weights are either 1 or calculated by in Eq. (4), it may also perform a second cycle in which weights are defined as $w_n = 1/\delta^\epsilon$, where $\delta = |E_n^{\text{exp}} - E_n^{\text{calc}}|$ and $\epsilon$ is an empirical tunable damping factor (generally $0.5 \le \epsilon \le 1$). The possibility to perform two cycles allows minoring the weight of data points affected by higher errors ($\sigma$) during the second cycle. This approach proved to be very useful when dealing with particularly complex systems [29,30]. Slight modifications of the BSTAC routines, developed by the same research group, allow the same kind of refinement by minimization of errors on volume readings instead of emf (STACO, see details in Ref. [31]), or to deal with measurements at different temperatures (TSTACO [32]). The different minimization approach used in STACO may enhance, if results are comparable with those obtained by BSTAC, the confidence on their reliability. Reversely,

significant differences in results by these two software may help to evidence possible systematic errors in the analyzed dataset. Concerning TSTACO, the analysis of potentiometric titrations at different temperatures (at the same ionic strength value) allows the simultaneous refinement of formation enthalpy changes and/or heat capacities, in addition to equilibrium constants.

The weighting scheme implemented in SUPERQUAD is based on the inverse of the variance calculated with the error propagation formula. This scheme has later been extended to HYPERQUAD as well. The error in the potential readings and the error in the volume are input parameters of the modeling. The calculation of the slope of the titration curve $\frac{\partial E}{\partial V}$ is done with a five-point cubic first-derivative convolution filter. Then, the variance is calculated by Eq. (4b).

KEV does not allow choosing how the weights of the experimental points are assigned and takes into consideration only the errors on the potential values. ReactLab™ pH PRO refines the equilibrium constants exclusively by unweighted least squares.

## 2.3. Ionic strength dependence

For a generic equilibrium expressed by Eq. (5), in which $a, b, c, d$ are stoichiometric coefficients, the stoichiometric (or apparent) overall equilibrium constant $\beta$ is given by Eq. (6).

$$a A + b B \rightleftharpoons c C + d D \quad (5)$$

$$\beta = \frac{[C]^c [D]^d}{[A]^a [B]^b} \quad (6)$$

By expressing the activity $a_s$ of a generic species S as a function of the activity coefficient $\gamma_S$ (Eq. (7)), the relation between the stoichiometric constant $\beta$ and the thermodynamic constant $\beta^T$ can be expressed by Eq. (8).

$$a_S = \gamma_S [S] \quad (7)$$

$$\beta^T = \frac{a_C^c a_D^d}{a_A^a a_B^b} = \frac{\gamma_C^c \gamma_D^d [C]^c [D]^d}{\gamma_A^a \gamma_B^b [A]^a [B]^b} = \frac{\gamma_C^c \gamma_D^d}{\gamma_A^a \gamma_B^b} \beta \quad (8)$$

Since the activity coefficients depend on the experimental conditions, such as solvent, ionic strength, temperature and pressure, the same is true for the stoichiometric equilibrium constants. Semi-empirical models can reproduce to some extend the dependence of activity coefficients on ionic strength conditions. The best known theories are those of Debye-Hückel [33], Davies [34], Bromley [35], Pitzer [36, 37], and the Specific ion Interaction Theory (SIT) [38,39]. Each approach leads to a mathematical expression that relates the activity coefficients to the ionic strength of the solution, but the validity range and complexity of each model is highly variable. For a simple 1:1 strong electrolyte solution, the upper concentration limit varies from about 0.001 mol kg⁻¹ for the simplest Debye-Hückel equation, to about ca. 0.1 mol kg⁻¹ for the Davies equation, to 3.5 mol kg⁻¹ for the SIT and, finally, well above and up to saturation for the Pitzer model. To apply these equations, all concentrations should be expressed in the molal scale (moles of solute per kilogram of pure solvent). In practice, these equations can be safely used also in the most common molar scale (moles of solute per liter of solutions) up to 0.1 mol dm⁻³, as the introduced error is typically in the one percent range. Above this threshold value, molal concentrations should be used. Through these equations, it is possible to define the activity coefficients for the charged species and, therefore, define the thermodynamic constants, provided the selected equation applies throughout the entire concentration range explored during the titration.

Some software can work with two different modes.

i) in the concentration mode, the ionic strength has to be considered as constant during the titration process (i.e. a salt buffer is introduced in
the titrant). The effect of the two approaches on the error sum was widely discussed in Ref. [5]. Two examples of the minimized *U* function are given by Eqs. (1) and (2) [9,14,28].

$$U = \sqrt{\frac{\sum w_n (V_n^{\text{exp}} - V_n^{\text{calc}})^2}{N - p}} \quad (1)$$

or

$$U = \sqrt{\frac{\sum w_n (E_n^{\text{exp}} - E_n^{\text{calc}})^2}{N - p}} \quad (2)$$

In both equations, $V_n^{\text{exp}}$ and $V_n^{\text{calc}}$, and $E_n^{\text{exp}}$ and $E_n^{\text{calc}}$ are the experimental and calculated volumes and potentials readings, respectively, corresponding to each titration point; $N - p$ is the number of degrees of freedom calculated as the difference between the number of experimental points ($N$) and the number of refined parameters ($p$), less the eventual constrains; $w_n$ is the weight associated to each titration point as defined in Eq. (3).

$$w_n = \frac{1}{s^2} \quad (3)$$

with

$$s^2 = s_V^2 + \left(\frac{\partial V}{\partial E}\right)^2 \times s_E^2 \quad (4a)$$

or

$$s^2 = s_E^2 + \left(\frac{\partial E}{\partial V}\right)^2 \times s_V^2 \quad (4b)$$

$s^2$ is the calculated variance of each titration point; $s_E^2$ and $s_V^2$ are the estimated variances of the potential readings and added volume, respectively, taken individually; and $\frac{\partial E}{\partial V}$ is the slope of the titration curve, calculated numerically for each pair of titration points. The weight assigned to each titration point is therefore inversely proportional to the associated variance. Therefore, the data near an equivalence point commonly affected by a high uncertainty and thus of low significance, have a large $\frac{\partial E}{\partial V}$ value and have a lower weight in the calculations than the other data in a buffered pH region. The estimated *U* value corresponds to the uncertainty of the fit and defines the uncertainty of the refined parameters through the error propagation law. Consequently, the calculation process used for refining the parameters, along with the variances entered as known values (usually $s_E^2$ and $s_V^2$), affect the uncertainty of the refined parameters. An unweighted refinement can also be as easily performed by setting the first constant term of Eq. (4) to one and the second term to zero (or to a value equal or lower to $10^{-10}$ in the case of HYPERQUAD).

The tested software manage the weight assignment to each titration point differently. BSTAC allows working with both weighted and unweighted refinement. As an added peculiarity, BSTAC offers also the possibility to execute two refinement cycles. After the first one, in which assigned weights are either 1 or calculated by in Eq. (4), it may also perform a second cycle in which weights are defined as $w_n = 1/\delta^\epsilon$, where $\delta = |E_n^{\text{exp}} - E_n^{\text{calc}}|$ and $\epsilon$ is an empirical tunable damping factor (generally $0.5 \le \epsilon \le 1$). The possibility to perform two cycles allows minoring the weight of data points affected by higher errors ($\sigma$) during the second cycle. This approach proved to be very useful when dealing with particularly complex systems [29,30]. Slight modifications of the BSTAC routines, developed by the same research group, allow the same kind of refinement by minimization of errors on volume readings instead of emf (STACO, see details in Ref. [31]), or to deal with measurements at different temperatures (TSTACO [32]). The different minimization approach used in STACO may enhance, if results are comparable with those obtained by BSTAC, the confidence on their reliability. Reversely,

significant differences in results by these two software may help to evidence possible systematic errors in the analyzed dataset. Concerning TSTACO, the analysis of potentiometric titrations at different temperatures (at the same ionic strength value) allows the simultaneous refinement of formation enthalpy changes and/or heat capacities, in addition to equilibrium constants.

The weighting scheme implemented in SUPERQUAD is based on the inverse of the variance calculated with the error propagation formula. This scheme has later been extended to HYPERQUAD as well. The error in the potential readings and the error in the volume are input parameters of the modeling. The calculation of the slope of the titration curve $\frac{\partial E}{\partial V}$ is done with a five-point cubic first-derivative convolution filter. Then, the variance is calculated by Eq. (4b).

KEV does not allow choosing how the weights of the experimental points are assigned and takes into consideration only the errors on the potential values. ReactLab™ pH PRO refines the equilibrium constants exclusively by unweighted least squares.

## 2.3. Ionic strength dependence

For a generic equilibrium expressed by Eq. (5), in which $a, b, c, d$ are stoichiometric coefficients, the stoichiometric (or apparent) overall equilibrium constant $\beta$ is given by Eq. (6).

$$a A + b B \rightleftharpoons c C + d D \quad (5)$$

$$\beta = \frac{[C]^c [D]^d}{[A]^a [B]^b} \quad (6)$$

By expressing the activity $a_s$ of a generic species S as a function of the activity coefficient $\gamma_S$ (Eq. (7)), the relation between the stoichiometric constant $\beta$ and the thermodynamic constant $\beta^T$ can be expressed by Eq. (8).

$$a_S = \gamma_S [S] \quad (7)$$

$$\beta^T = \frac{a_C^c a_D^d}{a_A^a a_B^b} = \frac{\gamma_C^c \gamma_D^d [C]^c [D]^d}{\gamma_A^a \gamma_B^b [A]^a [B]^b} = \frac{\gamma_C^c \gamma_D^d}{\gamma_A^a \gamma_B^b} \beta \quad (8)$$

Since the activity coefficients depend on the experimental conditions, such as solvent, ionic strength, temperature and pressure, the same is true for the stoichiometric equilibrium constants. Semi-empirical models can reproduce to some extend the dependence of activity coefficients on ionic strength conditions. The best known theories are those of Debye-Hückel [33], Davies [34], Bromley [35], Pitzer [36, 37], and the Specific ion Interaction Theory (SIT) [38,39]. Each approach leads to a mathematical expression that relates the activity coefficients to the ionic strength of the solution, but the validity range and complexity of each model is highly variable. For a simple 1:1 strong electrolyte solution, the upper concentration limit varies from about 0.001 mol kg⁻¹ for the simplest Debye-Hückel equation, to about ca. 0.1 mol kg⁻¹ for the Davies equation, to 3.5 mol kg⁻¹ for the SIT and, finally, well above and up to saturation for the Pitzer model. To apply these equations, all concentrations should be expressed in the molal scale (moles of solute per kilogram of pure solvent). In practice, these equations can be safely used also in the most common molar scale (moles of solute per liter of solutions) up to 0.1 mol dm⁻³, as the introduced error is typically in the one percent range. Above this threshold value, molal concentrations should be used. Through these equations, it is possible to define the activity coefficients for the charged species and, therefore, define the thermodynamic constants, provided the selected equation applies throughout the entire concentration range explored during the titration.

Some software can work with two different modes.

i) in the concentration mode, the ionic strength has to be considered as constant during the titration process (i.e. a salt buffer is introduced in
the titrated solution in order to guarantee an approximately constant ionic strength), the activity coefficients are considered to remain unknown and are ignored, and thus only stoichiometric equilibrium constants ($\beta$), valid at the working ionic strength and for the considered supporting electrolyte, are refined;

ii) in the activity mode, the ionic strength cannot be considered constant during the titration process (i.e. the solution contains variable concentration of a charged species that strongly affect the ionic strength value) and the activity coefficients are estimated by one of the model equations listed above, therefore the thermodynamic constants ($\beta^T$) can be defined.

Among the evaluated software, only BSTAC and ReactLab™ pH PRO can handle both cases, but the application of the equations for the definition of the activity coefficients is quite different.

BSTAC uses the Extended Debye–Hückel equation (from here on denoted as EDH-1) reported as Eq. (9),

$$ \log \gamma = -z^2 A \frac{\sqrt{I}}{1 + B\sqrt{I}} + CI + DI^2 + EI^2 \quad (9) $$

where A and B are the so-called Debye–Hückel parameters (e.g. A = 0.51 mol⁻¹/² dm³/² and B = 1.5 mol⁻¹/² dm³/² at T = 298.15 K in aqueous solutions, with concentrations and ionic strengths expressed in the molar – mol dm⁻³ – concentration scale), and C, D and E are empirical parameters (the E term can generally be neglected for $I \le 1.0$ mol dm⁻³). The latter empirical parameters can be further split into two other parameters each (i.e. c₀, c₁, d₀, d₁, e₀, e₁), depending on both the stoichiometric coefficients and the charges of the species involved in the considered equilibrium (for a more detailed description of the calculation, see Ref. [24]). In this mode, the software needs additional inputs, namely the reference ionic strength (commonly defined as the one at which the experiments were conducted), the ionic strength values at which the equilibrium constants are introduced into the model, and the concentration of the ionic species that do not participate to the equilibrium, but affect the ionic strength of the solution. The calculation proceeds with the refinement of the equilibrium constants and of the empirical parameters (C, D, E in Eq. (9)) describing the ionic strength dependence. During the refinement process, the effect of ionic strength is thus accounted for $E^0$ (Eq. (10), see 2.4. Synthetic dataset used to evaluate the various programs paragraph), equilibrium constants and, consequently, species concentration. The output is constituted by the equilibrium constants at both the reference ionic strength and at a number of ionic strengths between infinite dilution and the maximum ionic strength calculated from the dataset, as well as the empirical ionic strength dependence parameters C, D, and E that were eventually refined, together with their uncertainties. It also calculates species concentrations at the actual ionic strength for each specific titration point.

ReactLab™ pH PRO offers the possibility to choose one among four models to define the activity coefficients, namely Debye–Hückel [33], Davies [34], Extended Debye–Hückel equation (EDH-2; different from the one implemented by BSTAC), and Specific ion Interaction Theory (SIT) [38,39]. Among these four possibilities, the Debye–Hückel model is the only one that does not require any empirical parameter for estimating the activity coefficient. For the EDH-2 model, a single empirical parameter per ion involved in an equilibrium is requested (radius of the hydrated ion). For Davies model the value of 0.3 for C parameter is the default value, but it can be changed by the user if desired. In case of the SIT model, the user has to provide the so-called specific ion interaction parameters $\epsilon(i,j)$ for each charged species i involved in the equilibrium, where j corresponds to any counter ion of opposite charge present in the solution at a relevant concentration level. If known, these parameters can be retrieved from the literature [40] and treated as constants. Unfortunately, the number of compiled $\epsilon(i,j)$ values is rather limited and those pertaining to new chemical systems under investigation are by

essence unknown. Furthermore, SIT model and its coefficients are referred to the molal concentration scale. If the titrations are conducted without addition of a supporting electrolyte to keep the ionic strength approximately constant and if the total ionic strength does not exceed 0.1 mol dm⁻³, it is advisable to use instead the non-specific Davies model that uses the same empirical coefficient for all considered ions for estimating their activity coefficient. The usefulness of the SIT model is therefore rather restricted to well-established chemical systems.

Given the equilibrium constants and the component concentration values initially introduced in the model, the software performs a fist evaluation of the species concentrations. Based on the calculated ion concentrations, the ionic strength of the system is estimated and its value is used to compute the activity coefficients of the species for each titration point. The calculation proceeds through the optimization of the thermodynamic equilibrium constants by the error minimization process. This procedure provides directly the thermodynamic equilibrium constants $\beta^T$ instead of the stoichiometric equilibrium constants $\beta$ at the working ionic strength. The calculation process is further described in Ref. [41].

## 2.4. Synthetic dataset used to evaluate the various programs

The dataset is composed of six simulated pH-metric titration curves describing the stepwise protonation equilibria of a hypothetic hexaprotic base (A⁶⁻) with protonation constant values of $log K_2^H = 10$; $log K_2^H = 8$; $log K_3^H = 6$; $log K_4^H = 4$; $log K_5^H = 3$, and $log K_6^H = 2$ (in $log K_5^H$, n = number of protons in the protonated form of A⁶⁻). The water self-dissociation constant $K_w$ was considered equal to $10^{-13.77}$ (a typical value for a 0.1 mol dm⁻³ alkali chloride or nitrate salt solution at 298.15 K). The curves were generated by PyES software [24], are composed of 95–101 points and are considered as the tested database for performing the calculations. The components of the chemical system were A⁶⁻ and H⁺. Three total concentrations of the acid H₆A were considered ($2 \times 10^{-3}$ mol dm⁻³, $4 \times 10^{-3}$ mol dm⁻³, $5 \times 10^{-3}$ mol dm⁻³) corresponding to concentrations of the proton in the titration vessel of $12 \times 10^{-3}$ mol dm⁻³, $24 \times 10^{-3}$ mol dm⁻³, and $30 \times 10^{-3}$ mol dm⁻³. For each concentration level, two replicated titration curves were considered. The replicates differ by a slight variation of the potential values. The titrant was defined as a generic strong base, with a concentration of 0.1 mol dm⁻³. The total initial volume was 25 cm³. A constant temperature of 298.15 K was considered and the ionic strength was set at 0.1 mol dm⁻³. The complete dataset can be downloaded from a data repository [42].

The Nernst law defines the relation between the emf recorded by the millivolt-meter and the proton concentration. In practice, the extended version given by Eq. (10) is often used, which includes two correction terms accounting for the liquid junction potentials in the acidic ($j_A[H⁺]$) and alkaline ($j_B[OH⁻] = j_BK_w/[H⁺]$) ranges.

$$ E = E^0 + S \log[H^+] + j_A[H^+] + j_B[OH^-] \quad (10) $$

In Eq. (10), $E^0$ and S stand for the formal potential of the measuring cell (combined pH electrode or simple glass electrode used in combination of a reference electrode) and the Nernstian slope, respectively. The glass electrode calibration, in terms of proton concentration (p[H] or pHc = -log [H⁺] to avoid any confusion with the IUPAC definition of pH, pH = -log $a_H^+$ = p[H] - log $\gamma_H^+$), is most conveniently achieved by titrating a strong acid solution of known concentration with a standardized strong base solution. It follows that the formal potential $E^0$ of Eq. (10) differs from the $E^0$ parameter of the Nernst equation, as it includes the activity coefficient of the proton ($E^0 = E^0 + S \log \gamma_H^+$) supposed to remain constant over the entire pH range. This procedure allows to define the values of the parameters of Eq. (10) by nonlinear least squares refinement, as reported by Braibanti et al. [43]. For the simulation considered herein, the calibration parameters were selected as follows: formal potential $E^0 = 405.0$ mV, Nernstian slope $S = -59.16$ mV, liquid junction potential coefficient $j_A = -64$ mV dm³ mol⁻¹. The $j_B$
the titrated solution in order to guarantee an approximately constant ionic strength), the activity coefficients are considered to remain unknown and are ignored, and thus only stoichiometric equilibrium constants (β), valid at the working ionic strength and for the considered supporting electrolyte, are refined;

ii) in the activity mode, the ionic strength cannot be considered constant during the titration process (i.e. the solution contains variable concentration of a charged species that strongly affect the ionic strength value) and the activity coefficients are estimated by one of the model equations listed above, therefore the thermodynamic constants (β⁻¹) can be defined.

Among the evaluated software, only BSTAC and ReactLab™ pH PRO can handle both cases, but the application of the equations for the definition of the activity coefficients is quite different.

BSTAC uses the Extended Debye–Hückel equation (from here on denoted as EDH–1) reported as Eq. (9),

$$ \log \gamma = -z^2 A \frac{\sqrt{I}}{1 + B\sqrt{I}} + CI + DI^3 + EI^2 \quad (9) $$

where A and B are the so-called Debye–Hückel parameters (e.g. A = 0.51 mol⁻¹/² dm³/² and B = 1.5 mol⁻¹/² dm³/² at T = 298.15 K in aqueous solutions, with concentrations and ionic strengths expressed in the molar – mol dm⁻³ – concentration scale), and C, D and E are empirical parameters (the E term can generally be neglected for I ≤ 1.0 mol dm⁻³). The latter empirical parameters can be further split into two other parameters each (i.e. c₀, c₁, d₀, d₁, e₀, e₁), depending on both the stoichiometric coefficients and the charges of the species involved in the considered equilibrium (for a more detailed description of the calculation, see Ref. [24]). In this mode, the software needs additional inputs, namely the reference ionic strength (commonly defined as the one at which the experiments were conducted), the ionic strength values at which the equilibrium constants are introduced into the model, and the concentration of the ionic species that do not participate to the equilibrium, but affect the ionic strength of the solution. The calculation proceeds with the refinement of the equilibrium constants and of the empirical parameters (C, D, E in Eq. (9)) describing the ionic strength dependence. During the refinement process, the effect of ionic strength is thus accounted for E⁰ (Eq. (10), see 2.4. Synthetic dataset used to evaluate the various programs paragraph), equilibrium constants and, consequently, species concentration. The output is constituted by the equilibrium constants at both the reference ionic strength and at a number of ionic strengths between infinite dilution and the maximum ionic strength calculated from the dataset, as well as the empirical ionic strength dependence parameters C, D, and E that were eventually refined, together with their uncertainties. It also calculates species concentrations at the actual ionic strength for each specific titration point.

ReactLab™ pH PRO offers the possibility to choose one among four models to define the activity coefficients, namely Debye–Hückel [33], Davies [34], Extended Debye–Hückel equation (EDH–2; different from the one implemented by BSTAC), and Specific ion Interaction Theory (SIT) [38,39]. Among these four possibilities, the Debye–Hückel model is the only one that does not require any empirical parameter for estimating the activity coefficient. For the EDH–2 model, a single empirical parameter per ion involved in an equilibrium is requested (radius of the hydrated ion). For Davies model the value of 0.3 for C parameter is the default value, but it can be changed by the user if desired. In case of the SIT model, the user has to provide the so-called specific ion interaction parameters ϵ(i,j) for each charged species i involved in the equilibrium, where j corresponds to any counter ion of opposite charge present in the solution at a relevant concentration level. If known, these parameters can be retrieved from the literature [40] and treated as constants. Unfortunately, the number of compiled ϵ(i,j) values is rather limited and those pertaining to new chemical systems under investigation are by

essence unknown. Furthermore, SIT model and its coefficients are referred to the molal concentration scale. If the titrations are conducted without addition of a supporting electrolyte to keep the ionic strength approximately constant and if the total ionic strength does not exceed 0.1 mol dm⁻³, it is advisable to use instead the non-specific Davies model that uses the same empirical coefficient for all considered ions for estimating their activity coefficient. The usefulness of the SIT model is therefore rather restricted to well-established chemical systems.

Given the equilibrium constants and the component concentration values initially introduced in the model, the software performs a fist evaluation of the species concentrations. Based on the calculated ion concentrations, the ionic strength of the system is estimated and its value is used to compute the activity coefficients of the species for each titration point. The calculation proceeds through the optimization of the thermodynamic equilibrium constants by the error minimization process. This procedure provides directly the thermodynamic equilibrium constants β⁻¹ instead of the stoichiometric equilibrium constants β at the working ionic strength. The calculation process is further described in Ref. [41].

## 2.4. Synthetic dataset used to evaluate the various programs

The dataset is composed of six simulated pH-metric titration curves describing the stepwise protonation equilibria of a hypothetic hexaprotic base (A⁶⁻) with protonation constant values of log Kₐ⁻¹ = 10; log K₂⁻¹ = 8; log K₃⁻¹ = 6; log K₄⁻¹ = 4; log K₅⁻¹ = 3, and log K₆⁻¹ = 2 (in log Kₙ⁻¹, n = number of protons in the protonated form of A⁶⁻). The water self-dissociation constant Kₗ was considered equal to 10⁻¹³.⁷⁷ (a typical value for a 0.1 mol dm⁻³ alkali chloride or nitrate salt solution at 298.15 K). The curves were generated by PyES software [24], are composed of 95–101 points and are considered as the tested database for performing the calculations. The components of the chemical system were A⁶⁻ and H⁺. Three total concentrations of the acid H₆A were considered (2 × 10⁻³ mol dm⁻³, 4 × 10⁻³ mol dm⁻³, 5 × 10⁻³ mol dm⁻³) corresponding to concentrations of the proton in the titration vessel of 12 × 10⁻³ mol dm⁻³, 24 × 10⁻³ mol dm⁻³, and 30 × 10⁻³ mol dm⁻³. For each concentration level, two replicated titration curves were considered. The replicates differ by a slight variation of the potential values. The titrant was defined as a generic strong base, with a concentration of 0.1 mol dm⁻³. The total initial volume was 25 cm³. A constant temperature of 298.15 K was considered and the ionic strength was set at 0.1 mol dm⁻³. The complete dataset can be downloaded from a data repository [42].

The Nernst law defines the relation between the emf recorded by the millivolt-meter and the proton concentration. In practice, the extended version given by Eq. (10) is often used, which includes two correction terms accounting for the liquid junction potentials in the acidic (jₐ[H⁺]) and alkaline (jₐ[OH⁻] = jₐKₗ/[H⁺]) ranges.

$$ E = E^0 + S \log[H^+] + j_A[H^+] + j_B[OH^-] \quad (10) $$

In Eq. (10), E⁰ and S stand for the formal potential of the measuring cell (combined pH electrode or simple glass electrode used in combination of a reference electrode) and the Nernstian slope, respectively. The glass electrode calibration, in terms of proton concentration (p[H] or pHₐ = –log [H⁺] to avoid any confusion with the IUPAC definition of pH, pH = –log aₐ⁺ = p[H] – log γₐ⁺), is most conveniently achieved by titrating a strong acid solution of known concentration with a standardized strong base solution. It follows that the formal potential E⁰ of Eq. (10) differs from the E⁰ parameter of the Nernst equation, as it includes the activity coefficient of the proton (E⁰ = E⁰ + S log γₐ⁺) supposed to remain constant over the entire pH range. This procedure allows to define the values of the parameters of Eq. (10) by nonlinear least squares refinement, as reported by Braibanti et al. [43]. For the simulation considered herein, the calibration parameters were selected as follows: formal potential E⁰ = 405.0 mV, Nernstian slope S = –59.16 mV, liquid junction potential coefficient jₐ = –64 mV dm³ mol⁻¹. The jₐ
value was considered as negligible. It should be noted that the calibration parameters, being estimated in the molar concentration scale, depends on the working ionic strength (i.e. $E^0$ depends on log $\gamma_H^+$). Therefore, the calibration process must be conducted at the same ionic strength and with the same background electrolyte as the associated titration.

If the data analysis explicitly takes into account the variation of the ionic strength during the titration process, the calibration parameters must also be set at the new ionic strength values, as done by the BSTAC calculation procedure (see 2.3. *Ionic strength dependence* paragraph). Differently, the ReactLab™ pH PRO authors propose to work in the activity mode. Since, in this case, the titrations are carried out without supporting electrolyte and the mass action law is used in the activity scale, the electrodes must be calibrated in activities, by using certified buffer solutions, or by acid/base titrations using the same equation as for the titration for estimating activity coefficients.

## 2.5. Data processing

The dataset was processed by the different software while keeping the concentration of the components constant and managing the errors as allowed by each IT tool. With BSTAC, SUPERQUAD, HYPERQUAD, and KEV, the equilibrium constants were refined by weighted least squares (see 2.2. *Error minimization* paragraph), whereas ReactLab™ pH PRO only allows unweighted least squares fitting. The standard deviations of the potential readings (emf values) and added volume plugged in the weighing scheme were as follows: $s_v = 0.005 \text{ cm}^{-3}$ and $s_E = 0.15 \text{ mV}$ for BSTAC and HYPERQUAD; $s_v = 0.003 \text{ cm}^{-3}$ and $s_E = 0.3 \text{ mV}$ for SUPERQUAD, and $s_E = 0.15 \text{ mV}$ for KEV (KEV considers only the errors on the potential values). These values correspond to those commonly used by the research groups involved in this work.

## 2.6. Perturbation of the dataset

As stated beforehand, the same dataset was also analyzed after introducing some perturbations in the input file that reflect common systematic errors or typesetting mistakes with the aim of evaluate their effect on the refined parameters:

(i) The first perturbation considered concerns the calibration parameters: the data were analyzed using two different formal potentials, which deviate by 1 mV.

(ii) The other perturbations considered are related to the purity of the reagents. Data were first processed by considering a base concentration slightly lower than the actual value of 0.1 mol dm⁻³ used to generate the artificial titration curves. This error mainly derives from inaccurate standardization of the titrant solution. The carbonation of the base solution due to the dissolution of atmospheric CO₂ is another important source of error, as this ageing-related process cannot be avoided even by keeping the flask under a nitrogen or argon atmosphere. Polyethylene flasks and tubing being not gastight, CO₂ slowly diffuses into the solution, while alkaline solutions should not be stored in glass bottles to prevent their contamination by silicates that behave as weak bases, like carbonate. The carbonate content of the base solution can be determined by a Gran titration, while the solution should be discarded for concentrations higher than 0.5–1%. The general effect of carbonate contamination is that the pH is slightly higher than it should be and this shall introduce an obvious error in the refined stability constant values. However, the influence of the CO₂⁻ concentration on the refined protonation constants has deliberately not been evaluated herein for the following reasons. In most cases, the effect of carbonate impurity becomes perceptible only near the end-point of a titration curve. If so,

protonation constants should not be significantly affected. Moreover, the amount of dissolved CO₂ and HCO₃⁻ cannot be strictly controlled in the acidic region, as CO₂ may escape from the solution and there is no way to consider the gas/liquid equilibrium during the modeling process. Hence, a systematic variation of the total CO₂⁻ concentration in the titrant solution might not accurately predict the consequences on the lowest protonation constants of A⁶⁻. It should be stressed here that the issues of carbonate contamination cannot be resolved by simply including the HCO₃⁻/CO₂·H₂O and CO₂⁻/HCO₃⁻ equilibria into the chemical model. Moreover, the effect is different from that resulting from an error in the total amount of H⁺ or OH⁻. In the latter case, a small adjustment of the concentration or number of millimoles will bring the observed and calculated end-points into coincidence with each other, but this does not happen in case of carbonate contamination. The same conclusions hold also for silicate contamination. In both cases, great care has to be devoted to avoid the presence of these interfering anions.

(iii) Next, the data were analyzed by considering concentrations of the hexaprotic acid in the titrated solution that were 2% lower than those used for simulating the titration curves. In order to prevent divergence during the fit, the total proton concentration was refined too, as usual in these conditions.

(iv) Lastly, the effect of ionic strength variations occurring during the titration process was considered. The data was processed with all programs by assuming a constant ionic strength value during the titration. In addition, in the case of BSTAC and ReactLab™ pH PRO, the calculation of the actual ionic strength and the estimation of the activity coefficients was performed, in order to correct for the ionic strength variations upon neutralization of H₆A. Table 1 summarizes the simulation conditions used in each case.

# 3. Results

## 3.1. Software features

Some key features were selected to better compare the merits of the different IT tools. Table 2 lists available functionalities for each software included in this review. Unfortunately, only the most modern applications, HYPERQUAD, KEV, and ReactLab™ pH PRO are compatible with modern PC's operating systems and have graphical interfaces that simplify their use. Most of the software are able to perform a global fit, this means that they can handle several experimental titration curves derived from different experimental conditions. This ability is not so crucial in the studied case, but it is essential for more complex chemical systems, especially when the speciation model is non-trivial nor firmly established. In the latter case, it is often necessary to perform titrations using different experimental conditions and to process all curves simultaneously with a unique model to be able to identify all the possible species formed in solution and to define the corresponding formation constants. BSTAC and ReactLab™ pH PRO are the only two software that can optimize both cumulative and stepwise constants together with the corresponding uncertainties. HYPERQUAD has a dialog box for the calculation of stepwise constants, and their errors, from the cumulative ones. The calculation of the standard deviation on the stepwise constant takes into account the covariance (COV) between the overall constants: $s^2(\log K) = s^2(\log \beta_1) + s^2(\log \beta_2) - 2\text{COV}(\log \beta_1, \log \beta_2)$ [26]. BSTAC and ReactLab™ pH PRO can take into account correction terms for liquid junction potentials in the pH electrode calibration function (Eq. (10)), although ReactLab™ pH PRO can only treat $j_A$ values, taking by default $j_B$ as equal to zero. Moreover, only BSTAC and ReactLab™ pH PRO can handle ionic strength variations along an experiment, but, as already stressed before, the calculation processes and final aims are quite different.
value was considered as negligible. It should be noted that the calibration parameters, being estimated in the molar concentration scale, depends on the working ionic strength (i.e. $E^0$ depends on $log \gamma_H^+$). Therefore, the calibration process must be conducted at the same ionic strength and with the same background electrolyte as the associated titration.

If the data analysis explicitly takes into account the variation of the ionic strength during the titration process, the calibration parameters must also be set at the new ionic strength values, as done by the BSTAC calculation procedure (see 2.3. *Ionic strength dependence* paragraph). Differently, the ReactLab™ pH PRO authors propose to work in the activity mode. Since, in this case, the titrations are carried out without supporting electrolyte and the mass action law is used in the activity scale, the electrodes must be calibrated in activities, by using certified buffer solutions, or by acid/base titrations using the same equation as for the titration for estimating activity coefficients.

## 2.5. Data processing

The dataset was processed by the different software while keeping the concentration of the components constant and managing the errors as allowed by each IT tool. With BSTAC, SUPERQUAD, HYPERQUAD, and KEV, the equilibrium constants were refined by weighted least squares (see 2.2. *Error minimization* paragraph), whereas ReactLab™ pH PRO only allows unweighted least squares fitting. The standard deviations of the potential readings (emf values) and added volume plugged in the weighing scheme were as follows: $s_v = 0.005 \text{ cm}^{-3}$ and $s_E = 0.15 \text{ mV}$ for BSTAC and HYPERQUAD; $s_v = 0.003 \text{ cm}^{-3}$ and $s_E = 0.3 \text{ mV}$ for SUPERQUAD, and $s_E = 0.15 \text{ mV}$ for KEV (KEV considers only the errors on the potential values). These values correspond to those commonly used by the research groups involved in this work.

## 2.6. Perturbation of the dataset

As stated beforehand, the same dataset was also analyzed after introducing some perturbations in the input file that reflect common systematic errors or typesetting mistakes with the aim of evaluate their effect on the refined parameters:

(i) The first perturbation considered concerns the calibration parameters: the data were analyzed using two different formal potentials, which deviate by 1 mV.

(ii) The other perturbations considered are related to the purity of the reagents. Data were first processed by considering a base concentration slightly lower than the actual value of 0.1 mol dm$^{-3}$ used to generate the artificial titration curves. This error mainly derives from inaccurate standardization of the titrant solution. The carbonation of the base solution due to the dissolution of atmospheric CO$_2$ is another important source of error, as this ageing-related process cannot be avoided even by keeping the flask under a nitrogen or argon atmosphere. Polyethylene flasks and tubing being not gastight, CO$_2$ slowly diffuses into the solution, while alkaline solutions should not be stored in glass bottles to prevent their contamination by silicates that behave as weak bases, like carbonate. The carbonate content of the base solution can be determined by a Gran titration, while the solution should be discarded for concentrations higher than 0.5–1%. The general effect of carbonate contamination is that the pH is slightly higher than it should be and this shall introduce an obvious error in the refined stability constant values. However, the influence of the CO$_3^{2-}$ concentration on the refined protonation constants has deliberately not been evaluated herein for the following reasons. In most cases, the effect of carbonate impurity becomes perceptible only near the end-point of a titration curve. If so,

protonation constants should not be significantly affected. Moreover, the amount of dissolved CO$_2$ and HCO$_3^-$ cannot be strictly controlled in the acidic region, as CO$_2$ may escape from the solution and there is no way to consider the gas/liquid equilibrium during the modeling process. Hence, a systematic variation of the total CO$_3^{2-}$ concentration in the titrant solution might not accurately predict the consequences on the lowest protonation constants of A$^6^-$. It should be stressed here that the issues of carbonate contamination cannot be resolved by simply including the HCO$_3^-$ / CO$_2$·H$_2$O and CO$_3^{2-}$ / HCO$_3^-$ equilibria into the chemical model. Moreover, the effect is different from that resulting from an error in the total amount of H$^+$ or OH$^-$. In the latter case, a small adjustment of the concentration or number of millimoles will bring the observed and calculated end-points into coincidence with each other, but this does not happen in case of carbonate contamination. The same conclusions hold also for silicate contamination. In both cases, great care has to be devoted to avoid the presence of these interfering anions.

(iii) Next, the data were analyzed by considering concentrations of the hexaprotic acid in the titrated solution that were 2% lower than those used for simulating the titration curves. In order to prevent divergence during the fit, the total proton concentration was refined too, as usual in these conditions.

(iv) Lastly, the effect of ionic strength variations occurring during the titration process was considered. The data was processed with all programs by assuming a constant ionic strength value during the titration. In addition, in the case of BSTAC and ReactLab™ pH PRO, the calculation of the actual ionic strength and the estimation of the activity coefficients was performed, in order to correct for the ionic strength variations upon neutralization of H$_6$A. Table 1 summarizes the simulation conditions used in each case.

# 3. Results

## 3.1. Software features

Some key features were selected to better compare the merits of the different IT tools. Table 2 lists available functionalities for each software included in this review. Unfortunately, only the most modern applications, HYPERQUAD, KEV, and ReactLab™ pH PRO are compatible with modern PC's operating systems and have graphical interfaces that simplify their use. Most of the software are able to perform a global fit, this means that they can handle several experimental titration curves derived from different experimental conditions. This ability is not so crucial in the studied case, but it is essential for more complex chemical systems, especially when the speciation model is non-trivial nor firmly established. In the latter case, it is often necessary to perform titrations using different experimental conditions and to process all curves simultaneously with a unique model to be able to identify all the possible species formed in solution and to define the corresponding formation constants. BSTAC and ReactLab™ pH PRO are the only two software that can optimize both cumulative and stepwise constants together with the corresponding uncertainties. HYPERQUAD has a dialog box for the calculation of stepwise constants, and their errors, from the cumulative ones. The calculation of the standard deviation on the stepwise constant takes into account the covariance (COV) between the overall constants: $s^2(\log K) = s^2(\log \beta_1) + s^2(\log \beta_2) - 2\text{COV}(\log \beta_1, \log \beta_2)$ [26]. BSTAC and ReactLab™ pH PRO can take into account correction terms for liquid junction potentials in the pH electrode calibration function (Eq. (10)), although ReactLab™ pH PRO can only treat $j_A$ values, taking by default $j_B$ as equal to zero. Moreover, only BSTAC and ReactLab™ pH PRO can handle ionic strength variations along an experiment, but, as already stressed before, the calculation processes and final aims are quite different.
<table><caption>**Table 1**
Reference and perturbed conditions used to analyze the dataset with the different software.</caption><thead><tr><th></th><th>Perturbation</th><th>Affected parameter</th><th>Unit</th><th>Reference condition</th><th>Perturbed condition</th></tr></thead><tbody><tr><td>Calibration parameters</td><td>Wrong formal potential</td><td>E<sup>0</sup></td><td>mV</td><td>405.0</td><td>406.0</td></tr><tr><td rowspan="2">Purity of the reagents</td><td>Wrong titrant titer</td><td>Concentration of the titrant</td><td>mol dm<sup>-3</sup></td><td>0.1000</td><td>0.0980</td></tr><tr><td>Wrong concentration of the weak acid</td><td>Concentration of the weak acid</td><td>mol dm<sup>-3</sup></td><td>2 &times; 10<sup>-3</sup></td><td>1.96 &times; 10<sup>-3</sup></td></tr><tr><td rowspan="3">Data processing</td><td rowspan="3">Ionic strength changes</td><td rowspan="3">I</td><td>mol dm<sup>-3</sup></td><td>4 &times; 10<sup>-3</sup></td><td>3.92 &times; 10<sup>-3</sup></td></tr><tr><td>mol dm<sup>-3</sup></td><td>5 &times; 10<sup>-3</sup></td><td>4.90 &times; 10<sup>-3</sup></td></tr><tr><td>mol dm<sup>-3</sup></td><td>0.1 fixed</td><td>0.1 variable</td></tr></tbody></table>

<table><caption>**Table 2**
Comparison of the key features of the different tested software.</caption><thead><tr><th>Features</th><th>BSTAC</th><th>HYPERQUAD</th><th>KEV</th><th>ReactLab™</th><th>SUPERQUAD</th></tr></thead><tbody><tr><td>Handle more than one titration (global fit)</td><td>✓</td><td>✓</td><td>X</td><td>✓</td><td>✓</td></tr><tr><td>Optimize the concentration of reagents</td><td>✓</td><td>✓</td><td>X</td><td>✓</td><td>✓</td></tr><tr><td>Optimize both global and stepwise constants</td><td>✓</td><td>X</td><td>X</td><td>✓</td><td>X</td></tr><tr><td>Allow to constrain the total concentration values</td><td>✓</td><td>✓<sup>a</sup></td><td>X</td><td>X</td><td>✓</td></tr><tr><td>Treat all the calibration parameters (Eq. (10))</td><td>✓</td><td>X</td><td>X</td><td>✓<sup>b</sup></td><td>X</td></tr><tr><td>Optimize the calibration parameters</td><td>✓</td><td>✓</td><td>X</td><td>✓<sup>b</sup></td><td>✓</td></tr><tr><td>Work with variable ionic strength</td><td>✓</td><td>X</td><td>X</td><td>✓</td><td>X</td></tr><tr><td>Deal with titrations performed by any kind of ISEs</td><td>✓</td><td>✓</td><td>✓</td><td>X</td><td>✓</td></tr><tr><td>Freeware</td><td>✓</td><td>X</td><td>✓</td><td>X</td><td>X</td></tr><tr><td>Graphical interface</td><td>X<sup>c</sup></td><td>✓</td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>Natively compatible with modern PC and operating systems</td><td>X</td><td>✓</td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>On-line support</td><td>✓<sup>d</sup></td><td>X</td><td>✓<sup>d</sup></td><td>✓<sup>d</sup></td><td>X</td></tr><tr><td>User manual</td><td>X</td><td>X<sup>e</sup></td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>Help file</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>Year of the last issued version</td><td>2000</td><td>2013</td><td>2021</td><td>2023</td><td>1985</td></tr></tbody></table>

<sup>a</sup> Up to the 2006 version, constraining the total concentration changes of two or more reagents was allowed, but this option is no longer available in the latest Hyperquad 2013 release.

<sup>b</sup> It is possible to provide but not to refine a j<sub>A</sub> value, moreover it has to be the same for all the analyzed titration curves. In turn, j<sub>B</sub> is not considered. A different set of electrode potential (E<sup>0</sup>) and slope (S) values can be defined or refined for each titration curve (see Eq. (10)).

<sup>c</sup> BSTAC allows a very old, DOS-type graphical presentation of experimental/calculated titration curves.

<sup>d</sup> BSTAC: https://github.com/Kastakin/BSTAC; KEV: https://k-ev.org/; ReactLab™: http://jplusconsulting.com/Contact-Us/.

<sup>e</sup> On-line instructions.

The computation time for processing the tested dataset varies between hundredths of a second and a few seconds for most software, with the exception of KEV, that requires longer computation times, with the difference that the computation resources used are those provided by the hosting party and not those of the local machine.

## 3.2. Protonation constants of the hexaprotic acid in the absence of systematic errors

The protonation constants of the hexaprotic acid were refined by the different software using the simulated dataset containing no systematic errors, in order to compare their performances. The adjusted values and the corresponding uncertainties are reported in Table 3. The uncertainty is expressed on both the overall log β<sup>H</sup> and stepwise log K<sup>H</sup> formation constants. In case the software only provides the uncertainty on log β<sup>H</sup>, or log K<sup>H</sup> values, the undefined uncertainties have been estimated by the error propagation law (Eq. (11)) while neglecting the covariances, although the refined log β<sup>H</sup> are often highly correlated. Since KEV does not allow handling more than one titration curve, the uncertainties considered in this case are the standard deviations of log β<sup>H</sup> values estimated independently for each titration.

$$s_{\log K_n^H} = \sqrt{\sum_{i}^{n} s_{\log \beta_i^H}^2} \qquad (11a)$$

$$s_{\log \beta_i^H} = \sqrt{\sum_{i}^{n} s_{\log K_i^H}^2} \qquad (11b)$$

The log K<sup>H</sup> values are graphically displayed in Fig. 1. All the values are in excellent agreement with each other. The width of the confidence

interval estimated by the standard deviation (ν = 4 degrees of freedom, probability level of 95%, α = 0.05) of log K<sup>H</sup> values obtained by the different software (Fig. 1 – red dotted lines) is lower or comparable to the uncertainty of the constant, expressed as the standard deviation of the fit, excepted for log K<sup>H</sup><sub>1</sub>. Therefore, the dispersion of the log K<sup>H</sup><sub>n</sub> values due to the use of different software can be considered as negligible in most cases. However, log K<sup>H</sup><sub>1</sub> and log K<sup>H</sup><sub>6</sub> values show greater dispersions than the other four protonation constants, as expected for such high and low constants with real values set at 10 and 2, respectively. Indeed, only few titration points show a percentage of H<sub>6</sub>A (with respect to the total concentration of the acid) higher than 20% software (Fig. 1S – Supplementary material file). Since at least 10–15 titration points that cover 20–80% of the maximal abundance of each particular species in solution should be analyzed to provide reliable equilibrium constants, the high uncertainty observed for log K<sup>H</sup><sub>6</sub> values is not unexpected.

The uncertainties estimated by the different software are quite similar, also in case of KEV, although the uncertainty calculation was carried out by a different process (see 2.2. Error minimization paragraph). Only the uncertainty estimated by HYPERQUAD and ReactLab™ pH PRO are, respectively, slightly higher and lower than those estimated by the other software (Fig. 2S – Supplementary material file). However, the order of magnitude is the same. It has to be noted that ReactLab™ pH PRO is the only software that assigns by default an equal weight (w = 1) to each data point. Hence, the effects of using a weighted vs. an un-weighted least squares minimization on the fitted parameters and their associated uncertainty was evaluated by BSTAC, which conveniently allows working with both weighing schemes. The uncertainty returned by the unweighted least squares method are generally lower, but the trend is the same as that observed for the results obtained by ReactLab™ pH PRO. The values are graphically displayed in Fig. 1S and listed in Table 3 (column “reference conditions”, weighted least squares), and
<table><caption>**Table 1**
Reference and perturbed conditions used to analyze the dataset with the different software.</caption><thead><tr><th></th><th>Perturbation</th><th>Affected parameter</th><th>Unit</th><th>Reference condition</th><th>Perturbed condition</th></tr></thead><tbody><tr><td rowspan="2">Calibration parameters</td><td>Wrong formal potential</td><td>E<sup>0</sup></td><td>mV</td><td>405.0</td><td>406.0</td></tr><tr><td>Wrong titrant titer</td><td>Concentration of the titrant</td><td>mol dm<sup>-3</sup></td><td>0.1000</td><td>0.0980</td></tr><tr><td rowspan="2">Purity of the reagents</td><td rowspan="2">Wrong concentration of the weak acid</td><td rowspan="2">Concentration of the weak acid</td><td>mol dm<sup>-3</sup></td><td>2 × 10<sup>-3</sup></td><td>1.96 × 10<sup>-3</sup></td></tr><tr><td>mol dm<sup>-3</sup></td><td>4 × 10<sup>-3</sup></td><td>3.92 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">Data processing</td><td rowspan="2">Ionic strength changes</td><td rowspan="2">I</td><td>mol dm<sup>-3</sup></td><td>0.1 fixed</td><td>0.1 variable</td></tr><tr><td>mol dm<sup>-3</sup></td><td>0.1 fixed</td><td>0.1 variable</td></tr></tbody></table>

<table><caption>**Table 2**
Comparison of the key features of the different tested software.</caption><thead><tr><th>Features</th><th>BSTAC</th><th>HYPERQUAD</th><th>KEV</th><th>ReactLab™</th><th>SUPERQUAD</th></tr></thead><tbody><tr><td>Handle more than one titration (global fit)</td><td>✓</td><td>✓</td><td>X</td><td>✓</td><td>✓</td></tr><tr><td>Optimize the concentration of reagents</td><td>✓</td><td>✓</td><td>X</td><td>✓</td><td>✓</td></tr><tr><td>Optimize both global and stepwise constants</td><td>✓</td><td>X</td><td>X</td><td>✓</td><td>X</td></tr><tr><td>Allow to constrain the total concentration values</td><td>✓</td><td>✓<sup>a</sup></td><td>X</td><td>X</td><td>✓</td></tr><tr><td>Treat all the calibration parameters (Eq. (10))</td><td>✓</td><td>X</td><td>X</td><td>✓<sup>b</sup></td><td>X</td></tr><tr><td>Optimize the calibration parameters</td><td>✓</td><td>✓</td><td>X</td><td>✓<sup>b</sup></td><td>✓</td></tr><tr><td>Work with variable ionic strength</td><td>✓</td><td>X</td><td>X</td><td>✓</td><td>X</td></tr><tr><td>Deal with titrations performed by any kind of ISEs</td><td>✓</td><td>✓</td><td>✓</td><td>X</td><td>✓</td></tr><tr><td>Freeware</td><td>✓</td><td>X</td><td>✓</td><td>X</td><td>X</td></tr><tr><td>Graphical interface</td><td>X<sup>c</sup></td><td>✓</td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>Natively compatible with modern PC and operating systems</td><td>X</td><td>✓</td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>On-line support</td><td>✓<sup>d</sup></td><td>X</td><td>✓<sup>d</sup></td><td>✓<sup>d</sup></td><td>X</td></tr><tr><td>User manual</td><td>X</td><td>X<sup>e</sup></td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>Help file</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>X</td></tr><tr><td>Year of the last issued version</td><td>2000</td><td>2013</td><td>2021</td><td>2023</td><td>1985</td></tr></tbody></table>

<sup>a</sup> Up to the 2006 version, constraining the total concentration changes of two or more reagents was allowed, but this option is no longer available in the latest Hyperquad 2013 release.

<sup>b</sup> It is possible to provide but not to refine a j<sub>A</sub> value, moreover it has to be the same for all the analyzed titration curves. In turn, j<sub>B</sub> is not considered. A different set of electrode potential (E<sup>0</sup>) and slope (S) values can be defined or refined for each titration curve (see Eq. (10)).

<sup>c</sup> BSTAC allows a very old, DOS-type graphical presentation of experimental/calculated titration curves.

<sup>d</sup> BSTAC: https://github.com/Kastakin/BSTAC; KEV: https://k-ev.org/; ReactLab™: http://jplusconsulting.com/Contact-Us/.

<sup>e</sup> On-line instructions.

The computation time for processing the tested dataset varies between hundredths of a second and a few seconds for most software, with the exception of KEV, that requires longer computation times, with the difference that the computation resources used are those provided by the hosting party and not those of the local machine.

## 3.2. Protonation constants of the hexaprotic acid in the absence of systematic errors

The protonation constants of the hexaprotic acid were refined by the different software using the simulated dataset containing no systematic errors, in order to compare their performances. The adjusted values and the corresponding uncertainties are reported in Table 3. The uncertainty is expressed on both the overall log β<sup>H</sup> and stepwise log K<sup>H</sup> formation constants. In case the software only provides the uncertainty on log β<sup>H</sup>, or log K<sup>H</sup> values, the undefined uncertainties have been estimated by the error propagation law (Eq. (11)) while neglecting the covariances, although the refined log β<sup>H</sup> are often highly correlated. Since KEV does not allow handling more than one titration curve, the uncertainties considered in this case are the standard deviations of log β<sup>H</sup> values estimated independently for each titration.

$$s_{\log K_n^H} = \sqrt{\sum_{i}^{n} s_{\log \beta_i^H}^2} \qquad (11a)$$

$$s_{\log \beta_i^H} = \sqrt{\sum_{i}^{n} s_{\log K_i^H}^2} \qquad (11b)$$

The log K<sup>H</sup> values are graphically displayed in Fig. 1. All the values are in excellent agreement with each other. The width of the confidence

interval estimated by the standard deviation (ν = 4 degrees of freedom, probability level of 95%, α = 0.05) of log K<sup>H</sup> values obtained by the different software (Fig. 1 – red dotted lines) is lower or comparable to the uncertainty of the constant, expressed as the standard deviation of the fit, excepted for log K<sup>H</sup><sub>1</sub>. Therefore, the dispersion of the log K<sup>H</sup><sub>n</sub> values due to the use of different software can be considered as negligible in most cases. However, log K<sup>H</sup><sub>1</sub> and log K<sup>H</sup><sub>6</sub> values show greater dispersions than the other four protonation constants, as expected for such high and low constants with real values set at 10 and 2, respectively. Indeed, only few titration points show a percentage of H<sub>6</sub>A (with respect to the total concentration of the acid) higher than 20% software (Fig. 1S – Supplementary material file). Since at least 10–15 titration points that cover 20–80% of the maximal abundance of each particular species in solution should be analyzed to provide reliable equilibrium constants, the high uncertainty observed for log K<sup>H</sup><sub>6</sub> values is not unexpected.

The uncertainties estimated by the different software are quite similar, also in case of KEV, although the uncertainty calculation was carried out by a different process (see 2.2. Error minimization paragraph). Only the uncertainty estimated by HYPERQUAD and ReactLab™ pH PRO are, respectively, slightly higher and lower than those estimated by the other software (Fig. 2S – Supplementary material file). However, the order of magnitude is the same. It has to be noted that ReactLab™ pH PRO is the only software that assigns by default an equal weight (w = 1) to each data point. Hence, the effects of using a weighted vs. an un-weighted least squares minimization on the fitted parameters and their associated uncertainty was evaluated by BSTAC, which conveniently allows working with both weighing schemes. The uncertainty returned by the unweighted least squares method are generally lower, but the trend is the same as that observed for the results obtained by ReactLab™ pH PRO. The values are graphically displayed in Fig. 1S and listed in Table 3 (column “reference conditions”, weighted least squares), and
<table><caption>**Table 3**
Protonation constants obtained by the different software in reference and perturbed conditions.</caption><thead><tr><th colspan="10">BSTAC</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>10.010 ± 0.003</td><td>10.010 ± 0.003</td><td>10.049 ± 0.006</td><td>10.049 ± 0.006</td><td>10.303 ± 0.009</td><td>10.303 ± 0.009</td><td>9.979 ± 0.005</td><td>9.979 ± 0.005</td></tr><tr><td>2</td><td>18.014 ± 0.005</td><td>8.004 ± 0.004</td><td>18.066 ± 0.008</td><td>8.018 ± 0.007</td><td>18.47 ± 0.01</td><td>8.17 ± 0.01</td><td>18.002 ± 0.009</td><td>8.022 ± 0.005</td></tr><tr><td>3</td><td>24.017 ± 0.006</td><td>6.003 ± 0.004</td><td>24.091 ± 0.001</td><td>6.020 ± 0.007</td><td>24.62 ± 0.02</td><td>6.14 ± 0.01</td><td>24.06 ± 0.01</td><td>6.061 ± 0.005</td></tr><tr><td>4</td><td>28.018 ± 0.007</td><td>4.002 ± 0.004</td><td>28.11 ± 0.01</td><td>4.020 ± 0.007</td><td>28.70 ± 0.02</td><td>4.08 ± 0.01</td><td>28.14 ± 0.02</td><td>4.075 ± 0.004</td></tr><tr><td>5</td><td>31.021 ± 0.007</td><td>3.003 ± 0.005</td><td>31.13 ± 0.01</td><td>3.022 ± 0.009</td><td>31.70 ± 0.02</td><td>3.01 ± 0.01</td><td>31.20 ± 0.02</td><td>3.058 ± 0.004</td></tr><tr><td>6</td><td>33.01 ± 0.01</td><td>1.99 ± 0.01</td><td>33.19 ± 0.02</td><td>2.06 ± 0.02</td><td>33.70 ± 0.03</td><td>1.99 ± 0.03</td><td>33.30 ± 0.02</td><td>2.104 ± 0.008</td></tr><tr><th colspan="10">HYPERQUAD</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.010 ± 0.003</td><td>10.010 ± 0.003</td><td>10.049 ± 0.006</td><td>10.049 ± 0.006</td><td>10.303 ± 0.009</td><td>10.303 ± 0.009</td><td>9.980 ± 0.005</td><td>9.980 ± 0.005</td></tr><tr><td>2</td><td>18.014 ± 0.005</td><td>8.004 ± 0.006</td><td>18.066 ± 0.008</td><td>8.02 ± 0.01</td><td>18.47 ± 0.01</td><td>8.17 ± 0.01</td><td>18.002 ± 0.009</td><td>8.02 ± 0.01</td></tr><tr><td>3</td><td>24.017 ± 0.006</td><td>6.003 ± 0.008</td><td>24.09 ± 0.01</td><td>6.02 ± 0.01</td><td>24.62 ± 0.02</td><td>6.14 ± 0.02</td><td>24.06 ± 0.01</td><td>6.06 ± 0.01</td></tr><tr><td>4</td><td>28.019 ± 0.007</td><td>4.002 ± 0.009</td><td>28.11 ± 0.01</td><td>4.02 ± 0.01</td><td>28.70 ± 0.02</td><td>4.08 ± 0.03</td><td>28.14 ± 0.02</td><td>4.08 ± 0.02</td></tr><tr><td>5</td><td>31.020 ± 0.007</td><td>3.00 ± 0.01</td><td>31.13 ± 0.01</td><td>3.02 ± 0.01</td><td>31.71 ± 0.02</td><td>3.01 ± 0.03</td><td>31.20 ± 0.02</td><td>3.06 ± 0.03</td></tr><tr><td>6</td><td>33.03 ± 0.01</td><td>2.01 ± 0.01</td><td>33.21 ± 0.02</td><td>2.08 ± 0.02</td><td>33.72 ± 0.03</td><td>2.01 ± 0.04</td><td>33.32 ± 0.02</td><td>2.12 ± 0.03</td></tr><tr><th colspan="10">KEV</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.004 ± 0.002</td><td>10.004 ± 0.002</td><td>10.030 ± 0.004</td><td>10.030 ± 0.004</td><td>10.25 ± 0.01</td><td>10.25 ± 0.01</td><td>–</td><td>–</td></tr><tr><td>2</td><td>18.007 ± 0.002</td><td>8.003 ± 0.003</td><td>18.049 ± 0.005</td><td>8.019 ± 0.006</td><td>18.43 ± 0.02</td><td>8.18 ± 0.02</td><td>–</td><td>–</td></tr><tr><td>3</td><td>24.009 ± 0.003</td><td>6.002 ± 0.004</td><td>24.069 ± 0.006</td><td>6.020 ± 0.008</td><td>24.58 ± 0.02</td><td>6.15 ± 0.02</td><td>–</td><td>–</td></tr><tr><td>4</td><td>28.011 ± 0.003</td><td>4.002 ± 0.004</td><td>28.088 ± 0.007</td><td>4.019 ± 0.009</td><td>28.67 ± 0.02</td><td>4.09 ± 0.03</td><td>–</td><td>–</td></tr><tr><td>5</td><td>31.012 ± 0.005</td><td>3.001 ± 0.006</td><td>31.11 ± 0.01</td><td>3.02 ± 0.01</td><td>31.70 ± 0.03</td><td>3.02 ± 0.04</td><td>–</td><td>–</td></tr><tr><td>6</td><td>33.02 ± 0.01</td><td>2.00 ± 0.01</td><td>33.18 ± 0.03</td><td>2.07 ± 0.03</td><td>33.65 ± 0.06</td><td>1.95 ± 0.07</td><td>–</td><td>–</td></tr><tr><th colspan="10">ReactLab™ pH PRO</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.005 ± 0.002</td><td>10.005 ± 0.002</td><td>10.031 ± 0.002</td><td>10.031 ± 0.002</td><td>10.250 ± 0.006</td><td>10.250 ± 0.006</td><td>9.963 ± 0.004</td><td>9.963 ± 0.004</td></tr><tr><td>2</td><td>18.007 ± 0.002</td><td>8.002 ± 0.001</td><td>18.049 ± 0.003</td><td>8.018 ± 0.002</td><td>18.431 ± 0.008</td><td>8.181 ± 0.005</td><td>17.967 ± 0.005</td><td>8.004 ± 0.003</td></tr><tr><td>3</td><td>24.009 ± 0.002</td><td>6.002 ± 0.001</td><td>24.068 ± 0.003</td><td>6.019 ± 0.002</td><td>24.577 ± 0.009</td><td>6.146 ± 0.005</td><td>24.012 ± 0.006</td><td>6.045 ± 0.003</td></tr><tr><td>4</td><td>28.010 ± 0.003</td><td>4.001 ± 0.002</td><td>28.087 ± 0.005</td><td>4.019 ± 0.003</td><td>28.68 ± 0.01</td><td>4.099 ± 0.006</td><td>28.082 ± 0.007</td><td>4.070 ± 0.003</td></tr><tr><td>5</td><td>31.012 ± 0.005</td><td>3.002 ± 0.004</td><td>31.109 ± 0.008</td><td>3.022 ± 0.006</td><td>31.67 ± 0.01</td><td>2.99 ± 0.01</td><td>31.122 ± 0.008</td><td>3.040 ± 0.004</td></tr><tr><td>6</td><td>32.99 ± 0.01</td><td>1.98 ± 0.01</td><td>33.17 ± 0.02</td><td>2.06 ± 0.02</td><td>33.69 ± 0.03</td><td>2.02 ± 0.03</td><td>33.23 ± 0.01</td><td>2.108 ± 0.009</td></tr><tr><th colspan="10">SUPERQUAD</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.007 ± 0.002</td><td>10.007 ± 0.002</td><td>10.035 ± 0.003</td><td>10.035 ± 0.003</td><td>10.260 ± 0.006</td><td>10.260 ± 0.006</td><td>9.974 ± 0.004</td><td>9.974 ± 0.004</td></tr><tr><td>2</td><td>18.010 ± 0.003</td><td>8.003 ± 0.003</td><td>18.054 ± 0.004</td><td>8.019 ± 0.005</td><td>18.437 ± 0.008</td><td>8.18 ± 0.01</td><td>17.989 ± 0.007</td><td>8.016 ± 0.008</td></tr><tr><td>3</td><td>24.013 ± 0.003</td><td>6.003 ± 0.004</td><td>24.074 ± 0.005</td><td>6.020 ± 0.007</td><td>24.58 ± 0.01</td><td>6.14 ± 0.01</td><td>24.04 ± 0.01</td><td>6.06 ± 0.01</td></tr><tr><td>4</td><td>28.014 ± 0.004</td><td>4.002 ± 0.005</td><td>28.094 ± 0.006</td><td>4.020 ± 0.008</td><td>28.67 ± 0.01</td><td>4.09 ± 0.02</td><td>28.12 ± 0.01</td><td>4.08 ± 0.02</td></tr><tr><td>5</td><td>31.015 ± 0.004</td><td>3.001 ± 0.006</td><td>31.114 ± 0.007</td><td>3.02 ± 0.01</td><td>31.68 ± 0.01</td><td>3.00 ± 0.02</td><td>31.17 ± 0.01</td><td>3.05 ± 0.02</td></tr><tr><td>6</td><td>33.021 ± 0.009</td><td>2.01 ± 0.01</td><td>33.19 ± 0.01</td><td>2.08 ± 0.02</td><td>33.70 ± 0.03</td><td>2.03 ± 0.03</td><td>33.29 ± 0.02</td><td>2.13 ± 0.02</td></tr></tbody></table>

Table 4 (column “constant ionic strength”, unweighted least squares). A discussion on what is the best approach to be used for the regression analysis of potentiometric data fitting was reported by E. Casassas et al. [5]. The authors recommended weighted least squares refinement, whether the dependent variable corresponds to the potential readings or the added titrant volumes. However, when component concentrations are optimized, we would advise to work with unweighted least squares to give an equal weight to all titration points, including those close to the equivalent points that strictly relate to the component concentrations.

It should be noticed that the uncertainty estimation needs to be carefully handled if the software optimizes the values of the equilibrium constants rather than their logarithms. Assuming that the same titration process is repeated and that the titration data are subject to only random errors, a defined distribution of optimized log K values will be obtained. The same is true if the K value are optimized instead, but the type of the

distributions for the estimated log K and K values would not be the same, and the uncertainty evaluation should take into account the data distribution law.

### 3.3. Consequences of systematic errors

The bias introduced by an incorrect formal potential (E⁰) is quite pronounced. Indeed, the discrepancies in the protonation constants are already significant when the formal potential is 1 mV higher (equivalent to a 0.017 pH unit shift) than the correct value used to generate the artificial dataset. The differences between the log K<sup>H</sup><sub>n</sub> estimated by the different programs for both cases are reported in Fig. 2. The protonation constants more sensitive to the calibration conditions are those at the borders of the considered pH range, log K<sup>H</sup><sub>6</sub> and log K<sup>H</sup><sub>1</sub>; in any case, all the log K<sup>H</sup><sub>n</sub> values are overestimated on average by about 0.03 units,

<footnote>a The reported uncertainties are the standard deviations obtained by the calculation process (see 2.5. Data processing paragraph).</footnote>
<table><caption>**Table 3**
Protonation constants obtained by the different software in reference and perturbed conditions.</caption><thead><tr><th colspan="10">BSTAC</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>10.010 ± 0.003</td><td>10.010 ± 0.003</td><td>10.049 ± 0.006</td><td>10.049 ± 0.006</td><td>10.303 ± 0.009</td><td>10.303 ± 0.009</td><td>9.979 ± 0.005</td><td>9.979 ± 0.005</td></tr><tr><td>2</td><td>18.014 ± 0.005</td><td>8.004 ± 0.004</td><td>18.066 ± 0.008</td><td>8.018 ± 0.007</td><td>18.47 ± 0.01</td><td>8.17 ± 0.01</td><td>18.002 ± 0.009</td><td>8.022 ± 0.005</td></tr><tr><td>3</td><td>24.017 ± 0.006</td><td>6.003 ± 0.004</td><td>24.091 ± 0.001</td><td>6.020 ± 0.007</td><td>24.62 ± 0.02</td><td>6.14 ± 0.01</td><td>24.06 ± 0.01</td><td>6.061 ± 0.005</td></tr><tr><td>4</td><td>28.018 ± 0.007</td><td>4.002 ± 0.004</td><td>28.11 ± 0.01</td><td>4.020 ± 0.007</td><td>28.70 ± 0.02</td><td>4.08 ± 0.01</td><td>28.14 ± 0.02</td><td>4.075 ± 0.004</td></tr><tr><td>5</td><td>31.021 ± 0.007</td><td>3.003 ± 0.005</td><td>31.13 ± 0.01</td><td>3.022 ± 0.009</td><td>31.70 ± 0.02</td><td>3.01 ± 0.01</td><td>31.20 ± 0.02</td><td>3.058 ± 0.004</td></tr><tr><td>6</td><td>33.01 ± 0.01</td><td>1.99 ± 0.01</td><td>33.19 ± 0.02</td><td>2.06 ± 0.02</td><td>33.70 ± 0.03</td><td>1.99 ± 0.03</td><td>33.30 ± 0.02</td><td>2.104 ± 0.008</td></tr><tr><th colspan="10">HYPERQUAD</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.010 ± 0.003</td><td>10.010 ± 0.003</td><td>10.049 ± 0.006</td><td>10.049 ± 0.006</td><td>10.303 ± 0.009</td><td>10.303 ± 0.009</td><td>9.980 ± 0.005</td><td>9.980 ± 0.005</td></tr><tr><td>2</td><td>18.014 ± 0.005</td><td>8.004 ± 0.006</td><td>18.066 ± 0.008</td><td>8.02 ± 0.01</td><td>18.47 ± 0.01</td><td>8.17 ± 0.01</td><td>18.002 ± 0.009</td><td>8.02 ± 0.01</td></tr><tr><td>3</td><td>24.017 ± 0.006</td><td>6.003 ± 0.008</td><td>24.09 ± 0.01</td><td>6.02 ± 0.01</td><td>24.62 ± 0.02</td><td>6.14 ± 0.02</td><td>24.06 ± 0.01</td><td>6.06 ± 0.01</td></tr><tr><td>4</td><td>28.019 ± 0.007</td><td>4.002 ± 0.009</td><td>28.11 ± 0.01</td><td>4.02 ± 0.01</td><td>28.70 ± 0.02</td><td>4.08 ± 0.03</td><td>28.14 ± 0.02</td><td>4.08 ± 0.02</td></tr><tr><td>5</td><td>31.020 ± 0.007</td><td>3.00 ± 0.01</td><td>31.13 ± 0.01</td><td>3.02 ± 0.01</td><td>31.71 ± 0.02</td><td>3.01 ± 0.03</td><td>31.20 ± 0.02</td><td>3.06 ± 0.03</td></tr><tr><td>6</td><td>33.03 ± 0.01</td><td>2.01 ± 0.01</td><td>33.21 ± 0.02</td><td>2.08 ± 0.02</td><td>33.72 ± 0.03</td><td>2.01 ± 0.04</td><td>33.32 ± 0.02</td><td>2.12 ± 0.03</td></tr><tr><th colspan="10">KEV</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.004 ± 0.002</td><td>10.004 ± 0.002</td><td>10.030 ± 0.004</td><td>10.030 ± 0.004</td><td>10.25 ± 0.01</td><td>10.25 ± 0.01</td><td>–</td><td>–</td></tr><tr><td>2</td><td>18.007 ± 0.002</td><td>8.003 ± 0.003</td><td>18.049 ± 0.005</td><td>8.019 ± 0.006</td><td>18.43 ± 0.02</td><td>8.18 ± 0.02</td><td>–</td><td>–</td></tr><tr><td>3</td><td>24.009 ± 0.003</td><td>6.002 ± 0.004</td><td>24.069 ± 0.006</td><td>6.020 ± 0.008</td><td>24.58 ± 0.02</td><td>6.15 ± 0.02</td><td>–</td><td>–</td></tr><tr><td>4</td><td>28.011 ± 0.003</td><td>4.002 ± 0.004</td><td>28.088 ± 0.007</td><td>4.019 ± 0.009</td><td>28.67 ± 0.02</td><td>4.09 ± 0.03</td><td>–</td><td>–</td></tr><tr><td>5</td><td>31.012 ± 0.005</td><td>3.001 ± 0.006</td><td>31.11 ± 0.01</td><td>3.02 ± 0.01</td><td>31.70 ± 0.03</td><td>3.02 ± 0.04</td><td>–</td><td>–</td></tr><tr><td>6</td><td>33.02 ± 0.01</td><td>2.00 ± 0.01</td><td>33.18 ± 0.03</td><td>2.07 ± 0.03</td><td>33.65 ± 0.06</td><td>1.95 ± 0.07</td><td>–</td><td>–</td></tr><tr><th colspan="10">ReactLab™ pH PRO</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.005 ± 0.002</td><td>10.005 ± 0.002</td><td>10.031 ± 0.002</td><td>10.031 ± 0.002</td><td>10.250 ± 0.006</td><td>10.250 ± 0.006</td><td>9.963 ± 0.004</td><td>9.963 ± 0.004</td></tr><tr><td>2</td><td>18.007 ± 0.002</td><td>8.002 ± 0.001</td><td>18.049 ± 0.003</td><td>8.018 ± 0.002</td><td>18.431 ± 0.008</td><td>8.181 ± 0.005</td><td>17.967 ± 0.005</td><td>8.004 ± 0.003</td></tr><tr><td>3</td><td>24.009 ± 0.002</td><td>6.002 ± 0.001</td><td>24.068 ± 0.003</td><td>6.019 ± 0.002</td><td>24.577 ± 0.009</td><td>6.146 ± 0.005</td><td>24.012 ± 0.006</td><td>6.045 ± 0.003</td></tr><tr><td>4</td><td>28.010 ± 0.003</td><td>4.001 ± 0.002</td><td>28.087 ± 0.005</td><td>4.019 ± 0.003</td><td>28.68 ± 0.01</td><td>4.099 ± 0.006</td><td>28.082 ± 0.007</td><td>4.070 ± 0.003</td></tr><tr><td>5</td><td>31.012 ± 0.005</td><td>3.002 ± 0.004</td><td>31.109 ± 0.008</td><td>3.022 ± 0.006</td><td>31.67 ± 0.01</td><td>2.99 ± 0.01</td><td>31.122 ± 0.008</td><td>3.040 ± 0.004</td></tr><tr><td>6</td><td>32.99 ± 0.01</td><td>1.98 ± 0.01</td><td>33.17 ± 0.02</td><td>2.06 ± 0.02</td><td>33.69 ± 0.03</td><td>2.02 ± 0.03</td><td>33.23 ± 0.01</td><td>2.108 ± 0.009</td></tr><tr><th colspan="10">SUPERQUAD</th></tr><tr><th rowspan="2">n</th><th colspan="2">Reference condition</th><th colspan="2">Wrong formal potential</th><th colspan="2">Wrong titrant concentration</th><th colspan="2">Wrong acid concentration</th></tr><tr><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr><tr><td>1</td><td>10.007 ± 0.002</td><td>10.007 ± 0.002</td><td>10.035 ± 0.003</td><td>10.035 ± 0.003</td><td>10.260 ± 0.006</td><td>10.260 ± 0.006</td><td>9.974 ± 0.004</td><td>9.974 ± 0.004</td></tr><tr><td>2</td><td>18.010 ± 0.003</td><td>8.003 ± 0.003</td><td>18.054 ± 0.004</td><td>8.019 ± 0.005</td><td>18.437 ± 0.008</td><td>8.18 ± 0.01</td><td>17.989 ± 0.007</td><td>8.016 ± 0.008</td></tr><tr><td>3</td><td>24.013 ± 0.003</td><td>6.003 ± 0.004</td><td>24.074 ± 0.005</td><td>6.020 ± 0.007</td><td>24.58 ± 0.01</td><td>6.14 ± 0.01</td><td>24.04 ± 0.01</td><td>6.06 ± 0.01</td></tr><tr><td>4</td><td>28.014 ± 0.004</td><td>4.002 ± 0.005</td><td>28.094 ± 0.006</td><td>4.020 ± 0.008</td><td>28.67 ± 0.01</td><td>4.09 ± 0.02</td><td>28.12 ± 0.01</td><td>4.08 ± 0.02</td></tr><tr><td>5</td><td>31.015 ± 0.004</td><td>3.001 ± 0.006</td><td>31.114 ± 0.007</td><td>3.02 ± 0.01</td><td>31.68 ± 0.01</td><td>3.00 ± 0.02</td><td>31.17 ± 0.01</td><td>3.05 ± 0.02</td></tr><tr><td>6</td><td>33.021 ± 0.009</td><td>2.01 ± 0.01</td><td>33.19 ± 0.01</td><td>2.08 ± 0.02</td><td>33.70 ± 0.03</td><td>2.03 ± 0.03</td><td>33.29 ± 0.02</td><td>2.13 ± 0.02</td></tr></tbody></table>

Table 4 (column “constant ionic strength”, unweighted least squares). A discussion on what is the best approach to be used for the regression analysis of potentiometric data fitting was reported by E. Casassas et al. [5]. The authors recommended weighted least squares refinement, whether the dependent variable corresponds to the potential readings or the added titrant volumes. However, when component concentrations are optimized, we would advise to work with unweighted least squares to give an equal weight to all titration points, including those close to the equivalent points that strictly relate to the component concentrations.

It should be noticed that the uncertainty estimation needs to be carefully handled if the software optimizes the values of the equilibrium constants rather than their logarithms. Assuming that the same titration process is repeated and that the titration data are subject to only random errors, a defined distribution of optimized log K values will be obtained. The same is true if the K value are optimized instead, but the type of the

distributions for the estimated log K and K values would not be the same, and the uncertainty evaluation should take into account the data distribution law.

### 3.3. Consequences of systematic errors

The bias introduced by an incorrect formal potential (E⁰) is quite pronounced. Indeed, the discrepancies in the protonation constants are already significant when the formal potential is 1 mV higher (equivalent to a 0.017 pH unit shift) than the correct value used to generate the artificial dataset. The differences between the log K<sup>H</sup> estimated by the different programs for both cases are reported in Fig. 2. The protonation constants more sensitive to the calibration conditions are those at the borders of the considered pH range, log K<sup>H</sup><sub>6</sub> and log K<sup>H</sup><sub>1</sub>; in any case, all the log K<sup>H</sup> values are overestimated on average by about 0.03 units,

<footnote>a The reported uncertainties are the standard deviations obtained by the calculation process (see 2.5. *Data processing* paragraph).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** The log K<sub>n</sub><sup>H</sup> values estimated by the different software. The error bars represent the standard deviations estimated as explained in the paragraph 2.2 *Error minimization* and Eq. (11). Red dotted lines correspond to the mean values and their confidence interval (v = 4, α = 0.05). The titration curves were generated with log K<sub>n</sub><sup>H</sup> values 10, 8, 6, 4, 3 and 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

corresponding to relative errors of about 7% for K<sub>n</sub><sup>H</sup>. The consequences of
this perturbation are quite similar for the different tested software,
although log K<sub>1</sub><sup>H</sup> seems to be a little more sensitive (cf. The best estimates
returned by HYPERQUAD and BSTAC).

Decreasing the titrant concentration by 2% has even more critical consequences on the estimation of the protonation constants. The differences in log K<sub>n</sub><sup>H</sup> values are shown in Fig. 3. The deviation from the reference value increases from log K<sub>6</sub><sup>H</sup> to log K<sub>1</sub><sup>H</sup>, coherently with the volume of titrant added (see Fig. 1S – Supplementary material file). The protonation constants that can be determined upon addition of two equivalents of titrant (with respect to the acid concentration) are quite similar to the reference ones, whereas those that require larger amounts to be estimated are progressively affected by steadily increasing errors.

As for the condition described before, the sensitivity of the log K<sub>n</sub><sup>H</sup> to this perturbation is quite similar for the different tested software: only HYPERQUAD and BSTAC seem to slightly standout, as log K<sub>1</sub><sup>H</sup> is overestimated by 0.04 log units compared to the other software.

In order to show how the log K<sub>n</sub><sup>H</sup> changes with the titrant concentration error, the log K<sub>n</sub><sup>H</sup> values where optimized progressively increasing the percentage error. As an example, the results obtained by BSTAC are shown in Fig. 3S – Supplementary material file.

The effect of the purity of the reagent was evaluated by setting, for each fit, the effective total concentration of the hexaprotic acid as 2% lower than the real value. As mentioned above, the calculations were performed by refining simultaneously both the protonation constants and the total proton concentration for each curve. Therefore, the com-
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** The log K<sub>n</sub><sup>H</sup> values estimated by the different software. The error bars represent the standard deviations estimated as explained in the paragraph 2.2 *Error minimization* and Eq. (11). Red dotted lines correspond to the mean values and their confidence interval (v = 4, α = 0.05). The titration curves were generated with log K<sub>n</sub><sup>H</sup> values 10, 8, 6, 4, 3 and 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

corresponding to relative errors of about 7% for K<sub>n</sub><sup>H</sup>. The consequences of
this perturbation are quite similar for the different tested software,
although log K<sub>1</sub><sup>H</sup> seems to be a little more sensitive (cf. The best estimates
returned by HYPERQUAD and BSTAC).

Decreasing the titrant concentration by 2% has even more critical consequences on the estimation of the protonation constants. The differences in log K<sub>n</sub><sup>H</sup> values are shown in Fig. 3. The deviation from the reference value increases from log K<sub>6</sub><sup>H</sup> to log K<sub>1</sub><sup>H</sup>, coherently with the volume of titrant added (see Fig. 1S – Supplementary material file). The protonation constants that can be determined upon addition of two equivalents of titrant (with respect to the acid concentration) are quite similar to the reference ones, whereas those that require larger amounts to be estimated are progressively affected by steadily increasing errors.

As for the condition described before, the sensitivity of the log K<sub>n</sub><sup>H</sup> to this perturbation is quite similar for the different tested software: only HYPERQUAD and BSTAC seem to slightly standout, as log K<sub>1</sub><sup>H</sup> is overestimated by 0.04 log units compared to the other software.

In order to show how the log K<sub>n</sub><sup>H</sup> changes with the titrant concentration error, the log K<sub>n</sub><sup>H</sup> values where optimized progressively increasing the percentage error. As an example, the results obtained by BSTAC are shown in Fig. 3S – Supplementary material file.

The effect of the purity of the reagent was evaluated by setting, for each fit, the effective total concentration of the hexaprotic acid as 2% lower than the real value. As mentioned above, the calculations were performed by refining simultaneously both the protonation constants and the total proton concentration for each curve. Therefore, the com-
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Differences between the stepwise protonation constants log K<sub>n</sub><sup>H</sup> refined by considering E<sup>0</sup> = 405.0 mV (real value) and 406.0 mV.</figcaption></figure>

**Table 4**

Protonation constants (I = 0.1 mol dm<sup>-3</sup>) obtained by BSTAC and ReactLab™ pH PRO at constant and variable ionic strength. The weight assigned to each titration point was set equal to 1 and the Debye–Hückel equation (DH) was selected for estimating the activity coefficients.

<table><thead><tr><th rowspan="3"></th><th rowspan="3">n</th><th colspan="6">BSTAC</th></tr><tr><th colspan="4">I = 0.1 mol dm<sup>-3</sup></th><th colspan="2">I = 0.0 mol dm<sup>-3</sup></th></tr><tr><th colspan="2">Constant ionic strength</th><th colspan="2">Variable ionic strength</th><th colspan="2">Variable ionic strength</th></tr><tr><th></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>10.005 ± 0.002</td><td>10.005 ± 0.002</td><td>10.219 ± 0.005</td><td>10.219 ± 0.005</td><td>11.689 ± 0.006</td><td>11.689 ± 0.006</td></tr><tr><td>2</td><td>18.008 ± 0.002</td><td>8.003 ± 0.001</td><td>18.323 ± 0.007</td><td>8.104 ± 0.004</td><td>21.019 ± 0.008</td><td>9.330 ± 0.005</td></tr><tr><td>3</td><td>24.010 ± 0.002</td><td>6.002 ± 0.001</td><td>24.375 ± 0.006</td><td>6.052 ± 0.004</td><td>28.051 ± 0.007</td><td>7.032 ± 0.005</td></tr><tr><td>4</td><td>28.012 ± 0.003</td><td>4.002 ± 0.002</td><td>28.392 ± 0.007</td><td>4.017 ± 0.006</td><td>32.803 ± 0.009</td><td>4.753 ± 0.007</td></tr><tr><td>5</td><td>31.013 ± 0.004</td><td>3.012 ± 0.004</td><td>31.41 ± 0.01</td><td>3.02 ± 0.01</td><td>36.31 ± 0.02</td><td>3.50 ± 0.01</td></tr><tr><td>6</td><td>33.018 ± 0.009</td><td>2.01 ± 0.01</td><td>33.20 ± 0.05</td><td>1.79 ± 0.05</td><td>38.34 ± 0.05</td><td>2.04 ± 0.05</td></tr></tbody></table>

<table><thead><tr><th rowspan="3"></th><th rowspan="3">n</th><th colspan="6">ReactLab™ pH PRO</th></tr><tr><th colspan="4">I = 0.1 mol dm<sup>-3</sup></th><th colspan="2">I = 0.0 mol dm<sup>-3</sup></th></tr><tr><th colspan="2">Constant ionic strength</th><th colspan="2">Variable ionic strength<sup>a</sup></th><th colspan="2">Variable ionic strength</th></tr><tr><th></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>10.005 ± 0.002</td><td>10.005 ± 0.002</td><td>10.142 ± 0.003</td><td>10.142 ± 0.003</td><td>11.612 ± 0.003</td><td>11.612 ± 0.003</td></tr><tr><td>2</td><td>18.007 ± 0.002</td><td>8.002 ± 0.001</td><td>18.204 ± 0.003</td><td>8.062 ± 0.002</td><td>20.899 ± 0.004</td><td>9.287 ± 0.002</td></tr><tr><td>3</td><td>24.009 ± 0.002</td><td>6.002 ± 0.001</td><td>24.213 ± 0.002</td><td>6.009 ± 0.002</td><td>27.888 ± 0.003</td><td>6.989 ± 0.002</td></tr><tr><td>4</td><td>28.010 ± 0.003</td><td>4.001 ± 0.002</td><td>28.189 ± 0.003</td><td>3.976 ± 0.003</td><td>32.599 ± 0.004</td><td>4.711 ± 0.003</td></tr><tr><td>5</td><td>31.012 ± 0.005</td><td>3.002 ± 0.004</td><td>31.152 ± 0.006</td><td>2.963 ± 0.004</td><td>36.052 ± 0.006</td><td>3.453 ± 0.005</td></tr><tr><td>6</td><td>32.99 ± 0.01</td><td>1.98 ± 0.01</td><td>32.89 ± 0.02</td><td>1.73 ± 0.02</td><td>38.03 ± 0.02</td><td>1.98 ± 0.02</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Differences between the stepwise protonation constants log K<sub>n</sub><sup>H</sup> refined by considering a correct titrant concentration of 0.100 mol dm<sup>-3</sup> and an erroneous value of 0.098 mol dm<sup>-3</sup>.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Differences between the stepwise protonation constants log K<sub>n</sub><sup>H</sup> refined considering the correct and a 2% lower hexaprotic acid concentration.</figcaption></figure>

parison was done between the protonation constants estimated with the same calculation procedure for the datasets corresponding to the nominal and erroneously lower ligand concentrations. Fig. 4 shows the differences between the stepwise protonation constants obtained using these two calculations. The effect of this perturbation influences the most acidic constants and this likely depends on the fact that the error on the total ligand concentration induces an error on the total acid concentration (i.e. c<sub>H</sub> = 6 c<sub>L</sub>). Indeed, the deviations from the corresponding reference values increase from log K<sub>1</sub><sup>H</sup> to log K<sub>6</sub><sup>H</sup>, with HYPERQUAD and BSTAC providing the highest overestimations (Δlog K<sub>n</sub><sup>H</sup> = 0.05–0.15) among the four tested programs, except for log K<sub>1</sub><sup>H</sup>.

Finally, the unspoiled titration data were analyzed by considering the variation of the ionic strength of the solution during the titration process. Among the five-reviewed software, only BSTAC and ReactLab™ pH PRO offer that possibility. The approaches used for estimating the activity coefficients are not exactly the same, as already outlined in the

<footnote><sup>a</sup> Protonation constants recalculated at I = 0.1 mol dm<sup>-3</sup> with the DH equation from the refined thermodynamic equilibrium constants returned by the software.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Differences between the stepwise protonation constants log K<sub>n</sub><sup>H</sup> refined by considering E<sup>0</sup> = 405.0 mV (real value) and 406.0 mV.</figcaption></figure>

**Table 4**

Protonation constants (I = 0.1 mol dm<sup>-3</sup>) obtained by BSTAC and ReactLab™ pH PRO at constant and variable ionic strength. The weight assigned to each titration point was set equal to 1 and the Debye–Hückel equation (DH) was selected for estimating the activity coefficients.

<table><thead><tr><th rowspan="3"></th><th rowspan="3">n</th><th colspan="6">BSTAC</th></tr><tr><th colspan="4">I = 0.1 mol dm<sup>-3</sup></th><th colspan="2">I = 0.0 mol dm<sup>-3</sup></th></tr><tr><th colspan="2">Constant ionic strength</th><th colspan="2">Variable ionic strength</th><th colspan="2">Variable ionic strength</th></tr><tr><th></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>10.005 ± 0.002</td><td>10.005 ± 0.002</td><td>10.219 ± 0.005</td><td>10.219 ± 0.005</td><td>11.689 ± 0.006</td><td>11.689 ± 0.006</td></tr><tr><td>2</td><td>18.008 ± 0.002</td><td>8.003 ± 0.001</td><td>18.323 ± 0.007</td><td>8.104 ± 0.004</td><td>21.019 ± 0.008</td><td>9.330 ± 0.005</td></tr><tr><td>3</td><td>24.010 ± 0.002</td><td>6.002 ± 0.001</td><td>24.375 ± 0.006</td><td>6.052 ± 0.004</td><td>28.051 ± 0.007</td><td>7.032 ± 0.005</td></tr><tr><td>4</td><td>28.012 ± 0.003</td><td>4.002 ± 0.002</td><td>28.392 ± 0.007</td><td>4.017 ± 0.006</td><td>32.803 ± 0.009</td><td>4.753 ± 0.007</td></tr><tr><td>5</td><td>31.013 ± 0.004</td><td>3.012 ± 0.004</td><td>31.41 ± 0.01</td><td>3.02 ± 0.01</td><td>36.31 ± 0.02</td><td>3.50 ± 0.01</td></tr><tr><td>6</td><td>33.018 ± 0.009</td><td>2.01 ± 0.01</td><td>33.20 ± 0.05</td><td>1.79 ± 0.05</td><td>38.34 ± 0.05</td><td>2.04 ± 0.05</td></tr></tbody></table>

<table><thead><tr><th rowspan="3"></th><th rowspan="3">n</th><th colspan="6">ReactLab™ pH PRO</th></tr><tr><th colspan="4">I = 0.1 mol dm<sup>-3</sup></th><th colspan="2">I = 0.0 mol dm<sup>-3</sup></th></tr><tr><th colspan="2">Constant ionic strength</th><th colspan="2">Variable ionic strength<sup>a</sup></th><th colspan="2">Variable ionic strength</th></tr><tr><th></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th><th>log β<sup>H</sup><sub>n</sub></th><th>log K<sup>H</sup><sub>n</sub></th></tr></thead><tbody><tr><td>1</td><td>10.005 ± 0.002</td><td>10.005 ± 0.002</td><td>10.142 ± 0.003</td><td>10.142 ± 0.003</td><td>11.612 ± 0.003</td><td>11.612 ± 0.003</td></tr><tr><td>2</td><td>18.007 ± 0.002</td><td>8.002 ± 0.001</td><td>18.204 ± 0.003</td><td>8.062 ± 0.002</td><td>20.899 ± 0.004</td><td>9.287 ± 0.002</td></tr><tr><td>3</td><td>24.009 ± 0.002</td><td>6.002 ± 0.001</td><td>24.213 ± 0.002</td><td>6.009 ± 0.002</td><td>27.888 ± 0.003</td><td>6.989 ± 0.002</td></tr><tr><td>4</td><td>28.010 ± 0.003</td><td>4.001 ± 0.002</td><td>28.189 ± 0.003</td><td>3.976 ± 0.003</td><td>32.599 ± 0.004</td><td>4.711 ± 0.003</td></tr><tr><td>5</td><td>31.012 ± 0.005</td><td>3.002 ± 0.004</td><td>31.152 ± 0.006</td><td>2.963 ± 0.004</td><td>36.052 ± 0.006</td><td>3.453 ± 0.005</td></tr><tr><td>6</td><td>32.99 ± 0.01</td><td>1.98 ± 0.01</td><td>32.89 ± 0.02</td><td>1.73 ± 0.02</td><td>38.03 ± 0.02</td><td>1.98 ± 0.02</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Differences between the stepwise protonation constants log K<sub>n</sub><sup>H</sup> refined by considering a correct titrant concentration of 0.100 mol dm<sup>-3</sup> and an erroneous value of 0.098 mol dm<sup>-3</sup>.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Differences between the stepwise protonation constants log K<sub>n</sub><sup>H</sup> refined considering the correct and a 2% lower hexaprotic acid concentration.</figcaption></figure>

parison was done between the protonation constants estimated with the same calculation procedure for the datasets corresponding to the nominal and erroneously lower ligand concentrations. Fig. 4 shows the differences between the stepwise protonation constants obtained using these two calculations. The effect of this perturbation influences the most acidic constants and this likely depends on the fact that the error on the total ligand concentration induces an error on the total acid concentration (i.e. c<sub>H</sub> = 6 c<sub>L</sub>). Indeed, the deviations from the corresponding reference values increase from log K<sub>1</sub><sup>H</sup> to log K<sub>6</sub><sup>H</sup>, with HYPERQUAD and BSTAC providing the highest overestimations (Δlog K<sub>n</sub><sup>H</sup> = 0.05–0.15) among the four tested programs, except for log K<sub>1</sub><sup>H</sup>.

Finally, the unspoiled titration data were analyzed by considering the variation of the ionic strength of the solution during the titration process. Among the five-reviewed software, only BSTAC and ReactLab™ pH PRO offer that possibility. The approaches used for estimating the activity coefficients are not exactly the same, as already outlined in the

<footnote><sup>a</sup> Protonation constants recalculated at I = 0.1 mol dm<sup>-3</sup> with the DH equation from the refined thermodynamic equilibrium constants returned by the software.</footnote>
2.3. *Ionic strength dependence* paragraph. In order to make the results of the two software directly comparable, the calculations were conducted by selecting identical refinement parameters: the weight assigned to each titration point was set equal to 1, the same liquid junction potential value (jA = -64 mV dm³ mol⁻¹) was used, and the Debye–Hückel equation (DH) was applied to define the activity coefficients. To satisfy that last condition, the parameters of the EDH–1 equation implemented in BSTAC were set as follows: A = 0.51, B = 1, C = 0.0, D = 0.0, and E = 0.0, reducing it to the same DH expression used by ReactLab™ pH PRO. Moreover, E⁰ has been fitted too, as this parameters depends on the proton activity coefficient and, thus, on the ionic strength. Depending on the titration curve, the optimized E⁰ values ranged between 408.5 and 411.2 mV, which are quite different from the real one used to simulate the titration datasets (405 mV), but are in agreement with the activity coefficients of H⁺ at the working ionic strength [44] (E⁰ = E⁰ + S log γH⁺). As aforementioned in paragraph 2.3. *Ionic strength dependence*, the ReactLab™ pH PRO returns the thermodynamic constants (log K<sub>n</sub><sup>T,H</sup> at I = 0.0 mol dm<sup>-3</sup>) if the activity modality is used. Therefore, these values were extrapolated to I = 0.1 mol dm<sup>-3</sup> by the same DH equation to facilitate the comparison with the stoichiometric protonation constants obtained for a strictly constant ionic strength of I = 0.1 mol dm<sup>-3</sup> throughout the titrations. The results are reported in Table 4. It turns out that the medium effects between I = 0 and 0.1 mol dm<sup>-3</sup>, as roughly approximated by the DH model, are quite large for highly charged species, since log K<sub>n</sub><sup>T,H</sup> – log K<sub>n</sub><sup>H</sup> = 0.122 × z* = 1.46, 1.22, 0.98, 0.73, 0.49, 0.25 for 1 < n < 6 (z* = Σz<sub>reac</sub><sup>2</sup> – Σz<sub>prod</sub><sup>2</sup>). Davies model would be better for the considered ionic strength. Indeed log K<sub>n</sub><sup>T,H</sup> – log K<sub>n</sub><sup>H</sup> = 0.107 × z* = 1.28, 1.07, 0.86, 0.64, 0.43, 0.21 for 1 < n < 6 (z* = Σz<sub>reac</sub><sup>2</sup> – Σz<sub>prod</sub><sup>2</sup>). However, the effect of the ionic strength variations on the protonation constants would remain significant.

The deprotonation of a hexaprotic acid by a 1:1 strong base such as NaOH or KOH greatly affects the ionic strength of the solution along the titration. Indeed, the actual ionic strength values calculated by both software (an example is reported in Fig. 5a) range between 0.102 and 0.156 mol dm⁻³ along the titrations, showing maximum variation for the most concentrated solutions. Hence, the differences between the refined log Kₙ>H values, obtained by taking into account or not the ionic strength changes at each titration point, are far from being negligible for highly charged polyanions such as A⁶⁻ and AH⁵⁻ (Fig. 5b).

4. Discussion

Among the evaluated software, only the more modern or the com-
mercial packages are interactive through a graphical interface. User-
friendly tools facilitate their use and make them more appealing to
newcomers in the field, but also reduce the possibility of gross type-
setting errors at the input stage. It was one of the major purpose of this
work to assess the consequences of such errors in the input file on the
optimized equilibrium constants. Based on the results obtained herein, it
is possible to suggest a list of good practices that might limit the risk of
introducing systematic errors. The results highlight the relevance of a
correct electrode calibration [43] and, as a matter of fact, the necessity
to calibrate the electrode quite often. Since a variation of ca. 1 mV of the
formal electrode potential value within few consecutive days is not so
unusual, the suggestion is to alternate the calibration process to the
planned titrations during the working day and to use for data analysis
the calibration parameters obtained just before, or just after, the titra-
tion. Ideally, calibration should be run before and after each titration in
order to detect drifts. If some occur, the entire dataset should of course
be discarded. Moreover, the use of a control chart to monitor the trend of
the calibration parameters during the use of the electrode may help to
check the sensor performances.

Errors in titrant or component concentrations affect significantly the equilibrium constants values. Therefore, experimentalists should devote great attention to the purity of all reagents, and standardize the concentration of the stock and titrant solutions. Notably, the carbonation of alkaline titrants is the occurrence that most affects the estimated protonation constants (this can be even more important when determining stability constants of metal complexes, since metal carbonate species that may form are usually quite stable). It is therefore recommended to prepare brand new solutions frequently using freshly prepared ultrapure water (in the most critical cases, to boil it too) and to standardize the hydroxide solution used as titrant against dried, high purity potassium hydrogen phthalate. The Gran's method [45] can be used to check the carbonation level of alkaline titrants. Highly important too is the purity of the titrated compounds that has to be ascertained carefully before any solution equilibrium studies by combining a series of classical analysis (e.g. elemental C–H–N–S–O analyses, IR, mass, NMR, HPLC, ...).

Although the calculation procedures used by the refinement pro-
grams tested herein are not exactly the same, this work highlights no
significant differences between the sets of computed protonation con-
stants. Nevertheless, it is worth to note that each one has its own pros
and cons, some weaknesses and some desirable features, but often

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** (a) Blue line: simulated titration curve of 4 × 10⁻³ mol dm⁻³ hexaprotic acid solution, background salt 0.1 mol dm⁻³; red points: ionic strength calculated at each titration point as a function of the charged species in solution. (b) Differences between the stepwise protonation constants log Kₙ>H estimated at constant ionic strength (the weight assigned to each titration point was set equal to 1) and considering the actual ionic strength of each titration point. Activity coefficients were estimated by the Debye–Hückel equation. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>
2.3. *Ionic strength dependence* paragraph. In order to make the results of the two software directly comparable, the calculations were conducted by selecting identical refinement parameters: the weight assigned to each titration point was set equal to 1, the same liquid junction potential value (jA = -64 mV dm³ mol⁻¹) was used, and the Debye–Hückel equation (DH) was applied to define the activity coefficients. To satisfy that last condition, the parameters of the EDH–1 equation implemented in BSTAC were set as follows: A = 0.51, B = 1, C = 0.0, D = 0.0, and E = 0.0, reducing it to the same DH expression used by ReactLab™ pH PRO. Moreover, E⁰ has been fitted too, as this parameters depends on the proton activity coefficient and, thus, on the ionic strength. Depending on the titration curve, the optimized E⁰ values ranged between 408.5 and 411.2 mV, which are quite different from the real one used to simulate the titration datasets (405 mV), but are in agreement with the activity coefficients of H⁺ at the working ionic strength [44] (E⁰ = E⁰ + S log γH⁺). As aforementioned in paragraph 2.3. *Ionic strength dependence*, the ReactLab™ pH PRO returns the thermodynamic constants (log K<sub>n</sub><sup>T,H</sup> at I = 0.0 mol dm<sup>-3</sup>) if the activity modality is used. Therefore, these values were extrapolated to I = 0.1 mol dm<sup>-3</sup> by the same DH equation to facilitate the comparison with the stoichiometric protonation constants obtained for a strictly constant ionic strength of I = 0.1 mol dm<sup>-3</sup> throughout the titrations. The results are reported in Table 4. It turns out that the medium effects between I = 0 and 0.1 mol dm<sup>-3</sup>, as roughly approximated by the DH model, are quite large for highly charged species, since log K<sub>n</sub><sup>T,H</sup> – log K<sub>n</sub><sup>H</sup> = 0.122 × z* = 1.46, 1.22, 0.98, 0.73, 0.49, 0.25 for 1 < n < 6 (z* = Σz<sub>reac</sub><sup>2</sup> – Σz<sub>prod</sub><sup>2</sup>). Davies model would be better for the considered ionic strength. Indeed log K<sub>n</sub><sup>T,H</sup> – log K<sub>n</sub><sup>H</sup> = 0.107 × z* = 1.28, 1.07, 0.86, 0.64, 0.43, 0.21 for 1 < n < 6 (z* = Σz<sub>reac</sub><sup>2</sup> – Σz<sub>prod</sub><sup>2</sup>). However, the effect of the ionic strength variations on the protonation constants would remain significant.

The deprotonation of a hexaprotic acid by a 1:1 strong base such as NaOH or KOH greatly affects the ionic strength of the solution along the titration. Indeed, the actual ionic strength values calculated by both software (an example is reported in Fig. 5a) range between 0.102 and 0.156 mol dm⁻³ along the titrations, showing maximum variation for the most concentrated solutions. Hence, the differences between the refined log Kₙ>H values, obtained by taking into account or not the ionic strength changes at each titration point, are far from being negligible for highly charged polyanions such as A⁶⁻ and AH⁵⁻ (Fig. 5b).

4. Discussion

Among the evaluated software, only the more modern or the com-
mercial packages are interactive through a graphical interface. User-
friendly tools facilitate their use and make them more appealing to
newcomers in the field, but also reduce the possibility of gross type-
setting errors at the input stage. It was one of the major purpose of this
work to assess the consequences of such errors in the input file on the
optimized equilibrium constants. Based on the results obtained herein, it
is possible to suggest a list of good practices that might limit the risk of
introducing systematic errors. The results highlight the relevance of a
correct electrode calibration [43] and, as a matter of fact, the necessity
to calibrate the electrode quite often. Since a variation of ca. 1 mV of the
formal electrode potential value within few consecutive days is not so
unusual, the suggestion is to alternate the calibration process to the
planned titrations during the working day and to use for data analysis
the calibration parameters obtained just before, or just after, the titra-
tion. Ideally, calibration should be run before and after each titration in
order to detect drifts. If some occur, the entire dataset should of course
be discarded. Moreover, the use of a control chart to monitor the trend of
the calibration parameters during the use of the electrode may help to
check the sensor performances.

Errors in titrant or component concentrations affect significantly the equilibrium constants values. Therefore, experimentalists should devote great attention to the purity of all reagents, and standardize the concentration of the stock and titrant solutions. Notably, the carbonation of alkaline titrants is the occurrence that most affects the estimated protonation constants (this can be even more important when determining stability constants of metal complexes, since metal carbonate species that may form are usually quite stable). It is therefore recommended to prepare brand new solutions frequently using freshly prepared ultrapure water (in the most critical cases, to boil it too) and to standardize the hydroxide solution used as titrant against dried, high purity potassium hydrogen phthalate. The Gran's method [45] can be used to check the carbonation level of alkaline titrants. Highly important too is the purity of the titrated compounds that has to be ascertained carefully before any solution equilibrium studies by combining a series of classical analysis (e.g. elemental C–H–N–S–O analyses, IR, mass, NMR, HPLC, ...).

Although the calculation procedures used by the refinement pro-
grams tested herein are not exactly the same, this work highlights no
significant differences between the sets of computed protonation con-
stants. Nevertheless, it is worth to note that each one has its own pros
and cons, some weaknesses and some desirable features, but often

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** (a) Blue line: simulated titration curve of 4 × 10⁻³ mol dm⁻³ hexaprotic acid solution, background salt 0.1 mol dm⁻³; red points: ionic strength calculated at each titration point as a function of the charged species in solution. (b) Differences between the stepwise protonation constants log Kₙ>H estimated at constant ionic strength (the weight assigned to each titration point was set equal to 1) and considering the actual ionic strength of each titration point. Activity coefficients were estimated by the Debye–Hückel equation. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>
missing functionalities. In particular, only few software allow for a
complete treatment of input data, i.e., optimization of calibration pa-
rameters, introduction of constraints between the total concentrations to
be refined, management of the ionic strength fluctuations that invari-
ably occur during a continuous titration experiment. Each of these key
features is discussed in more detail hereafter.

i) The ability to optimize the electrode calibration parameters using the same software for analyzing titration data recorded for that purpose and for determining formation constants. In case of internal calibration procedure [46], this option becomes particularly useful, provided that the liquid junction potential parameters jA and jB can be refined together with the formal potential E0 and the slope S, according to Eq. (10).

ii) The possibility to impose refinement constrains between the total amounts (mmol) or concentrations of some reagents, such as those of the proton and protogenic species, ensures to reach convergence conditions in agreement with a hard modeling defined by stoichiometric laws. Moreover, the comparison between the results obtained by working with free and constrained concentrations may give information on the reagent's purity. However, the use of such constraints needs to be adequately pondered: it has to be taken into account that, in case the protogenic component is partially neutralized or residual acid is present as an impurity (cf. Formic or trifluoroacetic acid are often used as co-eluent in HPLC but they are difficult to remove in vacuo), this kind of modeling cannot be applied.

iii) On the basis of the results obtained by BSTAC and ReactLab™ pH PRO, it was possible to evaluate that the variation of the ionic strength during the titration process is not completely negligible even in the presence of a strong electrolyte. Ionic strength variation should be maintained within 10% of the actual ionic strength. It turns out that the most frequently used concentration of 0.1 mol dm⁻³ is often too low, particularly when working with highly charged species in solution. Therefore, having the ability to account for these variations by applying equations able to reliably estimate activity coefficients, should in principle afford more accurate equilibrium constant values. However, the selection of these equations in the calculation process requires further examinations and cautiousness in their application. Working at variable ionic strength requires selecting the best equation(s), depending on the ionic strength range in which the data are collected. Furthermore, the values of the numerical coefficients of the chosen equation must be selected with great care too, because their values change as a function of the species involved (charges and/or ionic radii) and of the medium [38]. For a background electrolyte concentration of 0.1 mol dm⁻³, the classical extended DH or Scatchard equations are no longer valid, while this concentration corresponds almost to the upper validity limit of the Davies equation assuming singly charged species. While none of these models includes unknown empirical parameters specific for the investigated species, this is unfortunately no longer true for the more complex models having a wider applicability range, like the SIT or any other virial expansions of the classical DH equation, to which the EDH-1 model belongs. Therefore, treating properly the ionic strength dependence of equilibrium constants for an unknown chemical system implies the experimental determination of these empirical coefficients by performing a series of measurements at different ionic strengths. Moreover, it must be taken into account that: i) the background salt, added to maintain the ionic strength more or less constant during the titration, allows to work in a conducting solution providing more stable signals and, thus, better electrode performances; ii) the calibration parameters of the electrode depend on the ionic

strength and on the chemical nature of background electrolytes
(for a glass electrode, E⁰ includes the activity coefficient of H⁺;
iii) the introduction of a new mathematical relation in the
computation process, as in the case of the equations for activity
coefficients estimations, increases the uncertainty of the opti-
mized parameters.

Performing potentiometric titrations in the absence of supporting electrolyte and thus variable but low ionic strength, as suggested in the manual of ReactLab™ pH PRO, should provide in principle a convenient access to thermodynamic equilibrium constants, as the activity coefficients should be more accurately estimated by the DH or Davies equation than in the presence of a high amount of salt. However, the titrated solution might have a low and variable electrical conductivity under such conditions, giving rise to large junction potential variations which are deleterious for stable and thus precise p[H] or pH<sub>c</sub> (±0.003 units), or emf measurements (±0.1 mV). Actually, only few commercially available ion-selective electrodes (ISEs) can work in very dilute solutions, while the most common devices require the addition of a supporting electrolyte to achieve the best performances. A discussion on the factors affecting the reproducibility of potentiometric sensors was reported by Lisak et al. [47]. In the activity mode, addition of a supporting electrolyte is also expected to become more and more problematic as the total ionic strength increases, because the DH correction terms applied to equilibrium constants for extrapolating them to I = 0 become larger and thus less reliable. For these reasons, this procedure, although appealing, cannot be recommended. A more careful validation endeavor has to be sought before envisaging the popularization of this approach.

All considered, the best practice in our view still remains to work with a medium as representative as possible of the real system under study, in order to express stoichiometric equilibrium constants in well-defined ionic medium and ionic strength conditions. If, during the titration process, significant ionic strength changes cannot be avoided, the equilibrium constants should be optimized by taking into account this variation. In the latter case, different considerations have to be taken into account: i) the best equation for the activity coefficient estimation has to be chosen (for more details see Ref. [38]), ii) the electrode calibration parameters should take into account the ionic strength changes (i.e. the expression of E⁰ should integrate the variations of the activity coefficient for the ion detected by the electrode), iii) experimentalists should process their data using both the variable and constant ionic strength mode and check to which extent the values of equilibrium constants are affected, and iv) declare what equation was used for estimating the activity coefficient and how it was parametrized.

5. Conclusions

In this work, the functionalities offered by five commercial or open-
access software commonly used to optimize protonation constants from
potentiometric data have been compared. To assess their performances,
the protonation constants of a hypothetical hexaanionic base of well-
defined acid-base properties have been refined by processing simu-
lated datasets. Despite the simplicity of the considered system, it was
possible to highlight some outcomes. The tested software presents
several differences, such as the type of interface and the refinement
possibilities, but all provided very similar sets of optimized protonation
constants. The analysis of the perturbed dataset highlights the impact of
some common errors in data entry or experimental processes on the
fitted parameters. In particular, the effect of ionic strength variations
during the titration process has to be stressed because it can most
significantly affect the optimized equilibrium constants values. This
aspect is important to be considered because large changes in the ionic
strength during the titration cannot be avoided in some instances, even
missing functionalities. In particular, only few software allow for a
complete treatment of input data, i.e., optimization of calibration pa-
rameters, introduction of constraints between the total concentrations to
be refined, management of the ionic strength fluctuations that invari-
ably occur during a continuous titration experiment. Each of these key
features is discussed in more detail hereafter.

i) The ability to optimize the electrode calibration parameters using the same software for analyzing titration data recorded for that purpose and for determining formation constants. In case of internal calibration procedure [46], this option becomes particularly useful, provided that the liquid junction potential parameters jA and jB can be refined together with the formal potential E0 and the slope S, according to Eq. (10).

ii) The possibility to impose refinement constrains between the total amounts (mmol) or concentrations of some reagents, such as those of the proton and protogenic species, ensures to reach convergence conditions in agreement with a hard modeling defined by stoichiometric laws. Moreover, the comparison between the results obtained by working with free and constrained concentrations may give information on the reagent's purity. However, the use of such constraints needs to be adequately pondered: it has to be taken into account that, in case the protogenic component is partially neutralized or residual acid is present as an impurity (cf. Formic or trifluoroacetic acid are often used as co-eluent in HPLC but they are difficult to remove in vacuo), this kind of modeling cannot be applied.

iii) On the basis of the results obtained by BSTAC and ReactLab™ pH PRO, it was possible to evaluate that the variation of the ionic strength during the titration process is not completely negligible even in the presence of a strong electrolyte. Ionic strength variation should be maintained within 10% of the actual ionic strength. It turns out that the most frequently used concentration of 0.1 mol dm⁻³ is often too low, particularly when working with highly charged species in solution. Therefore, having the ability to account for these variations by applying equations able to reliably estimate activity coefficients, should in principle afford more accurate equilibrium constant values. However, the selection of these equations in the calculation process requires further examinations and cautiousness in their application. Working at variable ionic strength requires selecting the best equation(s), depending on the ionic strength range in which the data are collected. Furthermore, the values of the numerical coefficients of the chosen equation must be selected with great care too, because their values change as a function of the species involved (charges and/or ionic radii) and of the medium [38]. For a background electrolyte concentration of 0.1 mol dm⁻³, the classical extended DH or Scatchard equations are no longer valid, while this concentration corresponds almost to the upper validity limit of the Davies equation assuming singly charged species. While none of these models includes unknown empirical parameters specific for the investigated species, this is unfortunately no longer true for the more complex models having a wider applicability range, like the SIT or any other virial expansions of the classical DH equation, to which the EDH-1 model belongs. Therefore, treating properly the ionic strength dependence of equilibrium constants for an unknown chemical system implies the experimental determination of these empirical coefficients by performing a series of measurements at different ionic strengths. Moreover, it must be taken into account that: i) the background salt, added to maintain the ionic strength more or less constant during the titration, allows to work in a conducting solution providing more stable signals and, thus, better electrode performances; ii) the calibration parameters of the electrode depend on the ionic

strength and on the chemical nature of background electrolytes
(for a glass electrode, E⁰ includes the activity coefficient of H⁺;
iii) the introduction of a new mathematical relation in the
computation process, as in the case of the equations for activity
coefficients estimations, increases the uncertainty of the opti-
mized parameters.

Performing potentiometric titrations in the absence of supporting electrolyte and thus variable but low ionic strength, as suggested in the manual of ReactLab™ pH PRO, should provide in principle a convenient access to thermodynamic equilibrium constants, as the activity coefficients should be more accurately estimated by the DH or Davies equation than in the presence of a high amount of salt. However, the titrated solution might have a low and variable electrical conductivity under such conditions, giving rise to large junction potential variations which are deleterious for stable and thus precise p[H] or pH<sub>c</sub> (±0.003 units), or emf measurements (±0.1 mV). Actually, only few commercially available ion-selective electrodes (ISEs) can work in very dilute solutions, while the most common devices require the addition of a supporting electrolyte to achieve the best performances. A discussion on the factors affecting the reproducibility of potentiometric sensors was reported by Lisak et al. [47]. In the activity mode, addition of a supporting electrolyte is also expected to become more and more problematic as the total ionic strength increases, because the DH correction terms applied to equilibrium constants for extrapolating them to I = 0 become larger and thus less reliable. For these reasons, this procedure, although appealing, cannot be recommended. A more careful validation endeavor has to be sought before envisaging the popularization of this approach.

All considered, the best practice in our view still remains to work with a medium as representative as possible of the real system under study, in order to express stoichiometric equilibrium constants in well-defined ionic medium and ionic strength conditions. If, during the titration process, significant ionic strength changes cannot be avoided, the equilibrium constants should be optimized by taking into account this variation. In the latter case, different considerations have to be taken into account: i) the best equation for the activity coefficient estimation has to be chosen (for more details see Ref. [38]), ii) the electrode calibration parameters should take into account the ionic strength changes (i.e. the expression of E⁰ should integrate the variations of the activity coefficient for the ion detected by the electrode), iii) experimentalists should process their data using both the variable and constant ionic strength mode and check to which extent the values of equilibrium constants are affected, and iv) declare what equation was used for estimating the activity coefficient and how it was parametrized.

5. Conclusions

In this work, the functionalities offered by five commercial or open-
access software commonly used to optimize protonation constants from
potentiometric data have been compared. To assess their performances,
the protonation constants of a hypothetical hexaanionic base of well-
defined acid-base properties have been refined by processing simu-
lated datasets. Despite the simplicity of the considered system, it was
possible to highlight some outcomes. The tested software presents
several differences, such as the type of interface and the refinement
possibilities, but all provided very similar sets of optimized protonation
constants. The analysis of the perturbed dataset highlights the impact of
some common errors in data entry or experimental processes on the
fitted parameters. In particular, the effect of ionic strength variations
during the titration process has to be stressed because it can most
significantly affect the optimized equilibrium constants values. This
aspect is important to be considered because large changes in the ionic
strength during the titration cannot be avoided in some instances, even
in the presence of a supporting electrolyte. This issue is very often overlooked, especially for chemical systems involving highly charged metal cations and/or polyelectrolyte ligands when investigated at rather low ionic strengths ($I = 0.1–0.2 \text{ mol dm}^{-3}$ range). Therefore, having numerical tools at our hands that are able to manage the variation of this parameter throughout the analyzed titration curves would be highly valuable and desirable as more accurate stoichiometric equilibrium constants should be obtained. Currently, this option is only implemented in two programs, the freeware but not widely spread BSTAC and the commercially available ReactLab™ pH PRO package. Although further developments and software improvements are sought, implementation of the ionic strength correction into the general minimization algorithm opens new perspectives in the field of solution complexation thermodynamics and might profoundly change some firmly established experimental practices. For decades now, titrations are performed in the presence of a supposedly inert background salt to keep the ionic strength and thus the activity coefficients constants. Working with a very low salt concentration ($I < 0.05 \text{ mol dm}^{-3}$), or even without, might appear as a straightforward mean for determining thermodynamic equilibrium constants if the mathematical treatment accurately estimates activity coefficients for each species. However, at the present stage of the knowledge, this unusual procedure cannot be recommended, as it requires a careful validation by expert laboratories in the frame of round-robin tests. A first weakness of this approach is that activity coefficients for neutral species are intrinsically assumed to equal 1, because only DH-type equations have so far been implemented in both BSTAC and ReactLab™ pH PRO. Another caveat of this approach, specific for potentiometric measurements, is related to the low salinity of the titrated solution if no extra salt is added. In the absence of supporting electrolyte, the liquid junction potential to the reference electrode is not kept constant, which gives rise to unstable and thus less accurate emf readings. Moreover, the entire set of emf data can no longer be treated by assuming a constant formal potential $E^0$, as the latter also depends on the proton activity coefficient.

It is clear that the scientific community working on the determination of equilibrium constants needs modern tools for data analysis. Considering the results obtained in this collegial project, the mathematical approach can be derived from the currently available IT products, but it would be desirable to improve their usability by proposing new tools offering more possibilities. Five major improvement directions have been identified. The future software should i) be most user-friendly with an interactive graphical interface, ii) be able to quickly handle a large amount of data, iii) refine simultaneously both calibration and titration data, iv) be supported by quality control tools, and v) properly manage the variation of the ionic strength during each titration.

### CRediT authorship contribution statement

**Silvia Berto:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Conceptualization.
**Salvador Blasco:** Software, Investigation.
**Lorenzo Castellino:** Writing – review & editing, Formal analysis, Data curation.
**Aleksandar Cvetkovski:** Writing – original draft.
**Concetta De Stefano:** Supervision, Software.
**Sofia Gama:** Software, Investigation.
**Enrique García-España:** Supervision, Software.
**Petr Hermann:** Supervision, Gabriele Lando: Writing – original draft, Software, Investigation.
**Matteo Marafante:** Visualization, Data curation.
**Michel Meyer:** Writing – review & editing, Supervision, Software.
**Winfried Plass:** Supervision.
**Lauryn Quinodo:** Investigation, Data curation.
**Demetrio Milea:** Writing – review & editing, Supervision, Resources, Project administration.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

The complete dataset can be downloaded at https://doi.org/10.17632/cfhkhmk6rt1

### Acknowledgments

This publication is based upon work from COST Action CA18202 NECTAR – Network for Equilibria and Chemical Thermodynamics Advanced Research, supported by COST (European Cooperation in Science and Technology, www_cost.eu, https://cost-nectar.eu/). SB, LC and MM also acknowledge support from the Project CH4.0 under the MUR program “Dipartimenti di Eccellenza 2023–2027” (CUP: D13C22003520001). SBI and EGE also acknowledge support from Conselleria de Innovación, Universidades, Ciencia y Sociedad Digital, Generalitat Valenciana (PROMETEO Grant CIPROM/2021/030).

### Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.aca.2024.342476.

### References

[1] D.M. Templeton, F. Ariese, R. Cornelis, L.-G. Danielsson, H. Muntau, H.P. Van Leeuwen, R. Łobiński, Clinical chemistry section, commission on toxicology guidelines for terms related to chemical speciation and fractionation of elements. Definitions, structural aspects, and methodological approaches, Pure Appl. Chem. 72 (2000) 1453–1470, https://doi.org/10.1351/pac200072081453.

[2] D. David, I. Nils, L.G. Sillén, Pit-mapping – a general approach for computer refining of equilibrium constants, Acta Chem. Scand. 15 (1961) 694–696, https://doi.org/10.3891/acta_chem_scand.15-0694.

[3] D.K. Nordstrom, J.W. Ball, Chemical models, computer programs and metal complexation in natural waters, in: Complexation Trace Met, Nat. Waters, 1984, pp. 149–164, https://doi.org/10.1007/978-94-009-6167-8_12.

[4] David J. Leggett (Ed.), Computational Methods for the Determination of Formation Constants, first ed., Springer, NY, New York, 2012 https://doi.org/10.1007/978-1-4684-4934-1.

[5] E. Casassas, R. Tauler, M. Filella, A critical comparison of computer programs for the potentiometric determination of stability constants, Anal. Chim. Acta 191 (1986) 399–411, https://doi.org/10.1016/S0003-2670(00)86325-4.

[6] K.A. Morin, Simplified explanations and examples of computerized methods for calculating chemical equilibrium in water, Comput. Geosci. 11 (1985) 409–416, https://doi.org/10.1016/0098-3004(85)90024-X.

[7] L.G. Sillen, High-speed computers as a supplement to graphical methods III. Twist matrix methods for minimizing the error-square sum in problems with many unknown constants, Acta Chem. Scand. 18 (1964) 1085–1098, https://doi.org/10.3891/acta_chem_scand.18-1085.

[8] J. Westall, MICROQL: I. A Chemical Equilibrium Program in Basic: II. Computation of Adsorption Equilibria in Basic, Swiss Federal Institute of Technology EAWAG, Dübendorf, Switzerland, 1979. https://www.worldcat.org/title/microql-i-a-chemical-equilibrium-program-in-basic-ii-computation-of-adsorption-equilibria-in-basic/oclc/636737294.

[9] P. Gans, A. Sabatini, A. Vacca, Superquad - a new computer program for determination of stability constants of complexes by potentiometric titration, Inorg. Chim. Acta. 79 (1983) 219–220, https://doi.org/10.1016/S0020-1693(00)95255-1.

[10] P.M. May, K. Murray, R.D. Willimas, The use of glass electrodes for the determination of formation constants-II Simulation of titration data, Talanta 32 (1985) 483–489, https://doi.org/10.1016/0039-9140(85)80261-7.

[11] R. Krebs, M. Sardin, D. Schweich, Mineral dissolution, precipitation, and ion exchange in surfactant flooding, AIChE J. 33 (1987) 1371–1378, https://doi.org/10.1002/AIC.690330814.

[12] J. van der Lee, Thermodynamic and mathematical concepts of chess stands for chemical equilibrium speciation with surfaces, Technical Report LHM/RD/98/39, CIG, École des Mines de Paris, Fontainebleau, France, September. https://api semanticscholar.org/CorpusID:18011323.
in the presence of a supporting electrolyte. This issue is very often overlooked, especially for chemical systems involving highly charged metal cations and/or polyelectrolyte ligands when investigated at rather low ionic strengths ($I = 0.1–0.2 \text{ mol dm}^{-3}$ range). Therefore, having numerical tools at our hands that are able to manage the variation of this parameter throughout the analyzed titration curves would be highly valuable and desirable as more accurate stoichiometric equilibrium constants should be obtained. Currently, this option is only implemented in two programs, the freeware but not widely spread BSTAC and the commercially available ReactLab™ pH PRO package. Although further developments and software improvements are sought, implementation of the ionic strength correction into the general minimization algorithm opens new perspectives in the field of solution complexation thermodynamics and might profoundly change some firmly established experimental practices. For decades now, titrations are performed in the presence of a supposedly inert background salt to keep the ionic strength and thus the activity coefficients constants. Working with a very low salt concentration ($I < 0.05 \text{ mol dm}^{-3}$), or even without, might appear as a straightforward mean for determining thermodynamic equilibrium constants if the mathematical treatment accurately estimates activity coefficients for each species. However, at the present stage of the knowledge, this unusual procedure cannot be recommended, as it requires a careful validation by expert laboratories in the frame of round-robin tests. A first weakness of this approach is that activity coefficients for neutral species are intrinsically assumed to equal 1, because only DH-type equations have so far been implemented in both BSTAC and ReactLab™ pH PRO. Another caveat of this approach, specific for potentiometric measurements, is related to the low salinity of the titrated solution if no extra salt is added. In the absence of supporting electrolyte, the liquid junction potential to the reference electrode is not kept constant, which gives rise to unstable and thus less accurate emf readings. Moreover, the entire set of emf data can no longer be treated by assuming a constant formal potential $E^0$, as the latter also depends on the proton activity coefficient.

It is clear that the scientific community working on the determination of equilibrium constants needs modern tools for data analysis. Considering the results obtained in this collegial project, the mathematical approach can be derived from the currently available IT products, but it would be desirable to improve their usability by proposing new tools offering more possibilities. Five major improvement directions have been identified. The future software should i) be most user-friendly with an interactive graphical interface, ii) be able to quickly handle a large amount of data, iii) refine simultaneously both calibration and titration data, iv) be supported by quality control tools, and v) properly manage the variation of the ionic strength during each titration.

### CRediT authorship contribution statement

**Silvia Berto:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Conceptualization.
**Salvador Blasco:** Software, Investigation.
**Lorenzo Castellino:** Writing – review & editing, Formal analysis, Data curation.
**Aleksandar Cvetkovski:** Writing – original draft.
**Concetta De Stefano:** Supervision, Software.
**Sofia Gama:** Software, Investigation.
**Enrique García-España:** Supervision, Software.
**Petr Hermann:** Supervision, Gabriele Lando: Writing – original draft, Software, Investigation.
**Matteo Marafante:** Visualization, Data curation.
**Michel Meyer:** Writing – review & editing, Supervision, Software.
**Winfried Plass:** Supervision.
**Lauryn Quinodo:** Investigation, Data curation.
**Demetrio Milea:** Writing – review & editing, Supervision, Resources, Project administration.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

The complete dataset can be downloaded at https://doi.org/10.17632/cfhkhmk6rt1

### Acknowledgments

This publication is based upon work from COST Action CA18202 NECTAR — Network for Equilibria and Chemical Thermodynamics Advanced Research, supported by COST (European Cooperation in Science and Technology, www_cost.eu, https://cost-nectar.eu/). SB, LC and MM also acknowledge support from the Project CH4.0 under the MUR program "Dipartimenti di Eccellenza 2023-2027" (CUP: D13C22003520001). SBI and EGE also acknowledge support from Conselleria de Innovación, Universidades, Ciencia y Sociedad Digital, Generalitat Valenciana (PROMETEO Grant CIPROM/2021/030).

### Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.aca.2024.342476.

### References

[1] D.M. Templeton, F. Ariese, R. Cornelis, L.-G. Danielsson, H. Muntau, H.P. Van Leeuwen, R. Łobiński, Clinical chemistry section, commission on toxicology guidelines for terms related to chemical speciation and fractionation of elements. Definitions, structural aspects, and methodological approaches, Pure Appl. Chem. 72 (2000) 1453–1470, https://doi.org/10.1351/pac200072081453.

[2] D. David, I. Nils, L.G. Sillén, Pit-mapping – a general approach for computer refining of equilibrium constants, Acta Chem. Scand. 15 (1961) 694–696, https://doi.org/10.3891/acta_chem_scand.15-0694.

[3] D.K. Nordstrom, J.W. Ball, Chemical models, computer programs and metal complexation in natural waters, in: Complexation Trace Met, Nat. Waters, 1984, pp. 149–164, https://doi.org/10.1007/978-94-009-6167-8_12.

[4] David J. Leggett (Ed.), Computational Methods for the Determination of Formation Constants, first ed., Springer, NY, New York, 2012 https://doi.org/10.1007/978-1-4684-4934-1.

[5] E. Casassas, R. Tauler, M. Filella, A critical comparison of computer programs for the potentiometric determination of stability constants, Anal. Chim. Acta 191 (1986) 399–411, https://doi.org/10.1016/S0003-2670(00)86325-4.

[6] K.A. Morin, Simplified explanations and examples of computerized methods for calculating chemical equilibrium in water, Comput. Geosci. 11 (1985) 409–416, https://doi.org/10.1016/0098-3004(85)90024-X.

[7] L.G. Sillen, High-speed computers as a supplement to graphical methods III. Twist matrix methods for minimizing the error-square sum in problems with many unknown constants, Acta Chem. Scand. 18 (1964) 1085–1098, https://doi.org/10.3891/acta_chem_scand.18-1085.

[8] J. Westall, MICROQL: I. A Chemical Equilibrium Program in Basic: II. Computation of Adsorption Equilibria in Basic, Swiss Federal Institute of Technology EAWAG, Dübendorf, Switzerland, 1979. https://www.worldcat.org/title/microql-i-a-chemical-equilibrium-program-in-basic-ii-computation-of-adsorption-equilibria-in-basic/oclc/636737294.

[9] P. Gans, A. Sabatini, A. Vacca, Superquad - a new computer program for determination of stability constants of complexes by potentiometric titration, Inorg. Chim. Acta. 79 (1983) 219–220, https://doi.org/10.1016/S0020-1693(00)95255-1.

[10] P.M. May, K. Murray, R.D. Willimas, The use of glass electrodes for the determination of formation constants-II Simulation of titration data, Talanta 32 (1985) 483–489, https://doi.org/10.1016/0039-9140(85)80261-7.

[11] R. Krebs, M. Sardin, D. Schweich, Mineral dissolution, precipitation, and ion exchange in surfactant flooding, AIChE J. 33 (1987) 1371–1378, https://doi.org/10.1002/AIC.690330814.

[12] J. van der Lee, Thermodynamic and mathematical concepts of chess stands for chemical equilibrium speciation with surfaces, Technical Report LHM/RD/98/39, CIG, École des Mines de Paris, Fontainebleau, France, September. https://api semanticscholar.org/CorpusID:18011323.