# Geometric-algebraic approach to aqueous solutions of diprotic acids and its
buffer mixtures

Juan C. Morales  [  Carlos A. Arango  [<EMAIL_ADDRESS>
###### Abstract

A closed-form analytical expression for $\ce{[H3O+]}$ has been obtained for
aqueous solutions of diprotic acids and its soluble salts. This formula allows
to calculate the pH of aqueous solutions of diprotic acids, their buffer
solutions, and the titrations of these two by a strong base, from the values
of p$K_{1}$, p$K_{2}$, and the effective concentrations of the acid and the
base, $\bar{C}_{\mathrm{a}}$ and $\bar{C}_{\mathrm{b}}$ respectively. It is
shown that a strong base titration of an acid, or its buffer solutions, is
always a linear path in the $\bar{C}_{\mathrm{a}}$–$\bar{C}_{\mathrm{b}}$
plane, which allows a simple analysis of the pH stability of buffer solutions.
The mathematical analysis of the equilibrium equations of the dissolution of a
diprotic acid in water and the physical constraints allowed to obtain two
approximate equations for the diprotic acids. One of the approximations is
useful for acids with $\mathrm{p}K_{2}-\mathrm{p}K_{1}\leq\log_{10}4$, the
other for acids with $\mathrm{p}K_{2}-\mathrm{p}K_{1}\leq-\log_{10}4$.

###### keywords:

diprotic weak acids, diprotic acid buffer solutions, acid-base titration,
buffer stability

Icesi]Department of Pharmaceutical and Chemical Sciences, Universidad Icesi,
Cali, Colombia Icesi]Department of Pharmaceutical and Chemical Sciences,
Universidad Icesi, Cali, Colombia

##  1 Introduction

Diprotic acids are of central importance in biochemistry, physiology, and
industrial and environmental chemistry. In biochemistry, several amino acids
behave as diprotic acids with two dissociated protons: one proton on the
$\alpha$ amino group and one on the $\alpha$ carboxyl group [1]. In
physiology, the regulation of blood pH cannot be understood without
considering the buffer made by carbonic acid, $\ce{H2CO3}$, and the
bicarbonate ion, $\ce{HCO3-}$ [2]. In environmental chemistry, the current
model for understanding ocean acidification is based on the aqueous chemical
equilibrium between $\ce{CO2}$, $\ce{H2CO3}$, $\ce{HCO3-}$, and $\ce{CO3-}$
[3].

A Brønsted diprotic acid is a chemical substance $\ce{H2B}$ that partially
dissolves in water producing hydronium ion, $\ce{H3O+}$, and the conjugate
base, $\ce{HB-}$. This conjugate base further dissociates partially producing
the second conjugate base, $\ce{B^{2-}}$. In the state of equilibrium, the
concentrations of the chemical species are constant [4, 5]. The equilibrium
concentrations of the chemical species are given by the equations of chemical
equilibrium, and the chemical and electric balance [6]. The aqueous
dissociation of a diprotic acid and, its soluble salts, involves five chemical
species and five mathematical relations between these species, therefore, in
principle is possible to obtain the equilibrium concentrations of all the
chemical species by solving this system of equations. In practice, the system
of equations involves nonlinear terms making difficult to obtain exact
mathematical expression for the equilibrium concentrations. For the
dissociation of a diprotic acid and its soluble salts, the algebraic
manipulation of the system of equations gives a quartic equation for the
concentration of $\ce{H3O+}$, $\ce{[H3O+]}$. The equilibrium concentration of
the hydronium ion is obtained by finding the roots of its quartic equation.
Although there is a quartic formula that gives the explicit roots of a quartic
equation, it is not practical to use due to its complexity. The use of the
quartic formulas gives four roots, each of them implies the execution of at
least 57 mathematical operations. Although these type of calculation is a
simple task for modern computers, the formulas obtained from the quartic
equation are not simplified which causes accumulation of computational error.
On the other hand, graphical and numerical solutions are easily obtained using
numerical calculators and computer software [7]. Although the graphical-
numerical approach is fast and efficient to calculate concentrations as
function of the pH, it has some disadvantages against an analytical closed
form solution. The analytical solution can be differentiated to study buffers
and buffer stability, or can be easily function-composed to study titrations
of acids and buffers by strong bases [8]. Another advantage of an analytical
closed form is the possibility to analyze mathematically the effect on the pH
of parameters such as the concentrations and the acid dissociation constants.
In this work it has been found that that the constrain
p$K_{2}-$p$K_{1}\geq\log_{10}{4}$, on the p$K$s of the acid, has an important
effect on the nature of the roots of the quartic polynomial for $\ce{[H3O+]}$.
This constrain has been previously obtained by considering isomerization of
the ionization products and a detailed equilibrium scheme in which the micro-
equilibrium constants correspond to partial equilibria [9, 10, 11]. Direct
observation of the experimental values of p$K_{1}$ and p$K_{2}$ of a large set
of diprotic acids allowed to find that several compounds, in particular the
nitrogenous organics, follow the constrain p$K_{2}-$p$K_{1}\leq-\log_{10}{4}$.

The main result of this paper is a closed form analytical solutions for [H3O+]
for the full chemical equilibrium of the aqueous dissociation of a diprotic
acid, and its monobasic and dibasic salts. The use of effective acid and base
concentrations allow to have only one mathematical expression for aqueous
dissolutions of diprotic acids, and buffer dissolutions of diprotic acids and
its soluble salts. In this work it is shown how this unified approach to
diprotic acids and its buffers allows to study the pH stability of buffer
solutions in relation with the equivalent acid dissolution.

This article is organized as follows: In the Theory and Methods section, the
first subsection is dedicated to establishing the notation and fundamental
equation of chemical equilibrium and physical constraints. In this subsection
a unified notation is introduced, allowing the same equations to be used for:
aqueous solutions of diprotic acids, buffers of diprotic acids, and titrations
with strong bases. In the second subsection of Theory and Methods, a
mathematical expression for $\ce{[H3O+]}$ is obtained and analyzed, showing
the complete expression for $\ce{[H3O+]}$. The final expressions of
$\ce{[H3O+]}$ are written in algebraic terms using basic arithmetic operations
and radicals (square and cube roots). These expressions can be used with the
help of computer algebra software to obtain the pH without approximations. The
reader interested in the mathematical details and the procedures used to
obtain these equations is referred to the Appendix. The third and final
subsection introduces exact expressions for the titration functions of the
aqueous solution and the buffer solution of diprotic acids. The Results and
Discussion section shows the results obtained using the expressions obtained
previously. The first section shows that although most diprotic acids obey the
condition $\mathrm{p}K_{2}-\mathrm{p}K_{1}\geq\log_{10}4$, there are some
acids that follow the condition
$\mathrm{p}K_{2}-\mathrm{p}K_{1}\leq-\log_{10}{4}$. In the next subsection,
the physical constraints of diprotic acids are used to obtain two
approximations for $\ce{[H3O+]}$, these approximations are used to obtain
analytical expressions for the upper and lower limits of pH. In the next
subsection, we discuss the common approach of monitoring the second
dissociation constant and show that in the case of micro-molar concentrations
this approach fails. The following subsection shows the use of exact closed
forms of pH and titration functions to analyze the neutralization of aqueous
solutions of diprotic acids and their buffer mixtures. The differences between
the exact expressions of this work and the approximate results of recent works
for two cases are shown in detail: maleic acid and 1,8-Octanediamine. Finally,
the last subsection of Results and Discussion shows an analysis of the pH
stability of diprotic acid buffer solutions. In this subsection, the pH
stability is analyzed as a parametric curve in the plane formed by the pH of
the acid and the pH of the corresponding buffer solution.

##  2 Theory and Methods

###  2.1 Aqueous solutions of weak diprotic acids and their salts

The aqueous dissociation equilibrium of a diprotic weak acid $\ce{H2B}$ is
given by the chemical equations

Relevant chemical species are $\ce{H3O+}$, $\ce{OH-}$, $\ce{H2B}$, $\ce{HB-}$,
and $\ce{B^{2-}}$ with equilibrium molar concentrations $\ce{[H3O+]}$,
$\ce{[OH^{-}]}$, $\ce{[H2B]}$, $\ce{[HB^{-}]}$, and $\ce{[B^{2-}]}$,
respectively. The equilibria displayed in equations (2.1)–(2.1) are effective
equilibria since the two protons of $\ce{H2B}$ can dissociate separately and
not necessarily consecutively [9, 1, 11].

A solution of the acid $\ce{H2B}$ is prepared in water at analytical molar
concentration $C_{\mathrm{a}}$. Once the system reaches chemical equilibrium,
the concentrations of the chemical species are given by five physical
conditions: two weak acid dissociation constant $K_{\mathrm{1}}$ and
$K_{\mathrm{2}}$, the water auto-ionization constant $K_{\mathrm{w}}$, the
electric neutrality, and the mass balance,

$\displaystyle K_{\mathrm{1}}$
$\displaystyle=\frac{\ce{[H3O+]}}{C^{\circ}}\frac{\ce{[HB^{-}]}}{C^{\circ}}\left(\frac{\ce{[H2B]}}{C^{\circ}}\right)^{-1},$
(1) $\displaystyle K_{\mathrm{2}}$
$\displaystyle=\frac{\ce{[H3O+]}}{C^{\circ}}\frac{\ce{[B^{2-}]}}{C^{\circ}}\left(\frac{\ce{[HB^{-}]}}{C^{\circ}}\right)^{-1},$
(2) $\displaystyle K_{\mathrm{w}}$
$\displaystyle=\frac{\ce{[H3O^{+}]}}{C^{\circ}}\frac{\ce{[OH^{-}]}}{C^{\circ}},$
(3) $\displaystyle\ce{[H3O+]}$
$\displaystyle=\ce{[OH^{-}]}+\ce{[HB^{-}]}+2\ce{[B^{2-}]},$ (4) $\displaystyle
C_{\mathrm{a}}$ $\displaystyle=\ce{[HB]}+\ce{[HB^{-}]}+\ce{[B^{2-}]},$ (5)

respectively. The standard molar concentration is ${C^{\circ}=1\,\mathrm{M}}$.
The acid constants $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ are dimensionless,
and their value range typically between $10^{-10}$ and $10^{-1}$. In this
work, the biochemical standard state
$C^{\standardstate}=C^{\circ}\sqrt{K_{\mathrm{w}}}$ is used to define the
dimensionless variables: ${x=\ce{[H3O^{+}]}/C^{\standardstate}}$,
${y=\ce{[OH^{-}]}/C^{\standardstate}}$,
${z_{0}=\ce{[H2B]}/C^{\standardstate}}$,
${z_{1}=\ce{[HB^{-}]}/C^{\standardstate}}$,
${z_{2}=\ce{[B^{2-}]}/C^{\standardstate}}$, and the parameter
${c_{\mathrm{a}}=C_{\mathrm{a}}/C^{\standardstate}}$. These definitions make
the equilibrium constants ${k_{1}=K_{\mathrm{1}}/\sqrt{K_{\mathrm{w}}}}$,
${k_{2}=K_{\mathrm{2}}/\sqrt{K_{\mathrm{w}}}}$, and ${k_{\mathrm{w}}=1}$. In
terms of the new variables and constants, equations (1)–(5) are replaced by

$\displaystyle k_{\mathrm{1}}$ $\displaystyle=\frac{xz_{1}}{z_{0}},$ (6)
$\displaystyle k_{\mathrm{2}}$ $\displaystyle=\frac{xz_{2}}{z_{1}},$ (7)
$\displaystyle k_{\mathrm{w}}$ $\displaystyle=xy=1,$ (8) $\displaystyle x$
$\displaystyle=y+z_{1}+2z_{2},$ (9) $\displaystyle c_{\mathrm{a}}$
$\displaystyle=z_{0}+z_{1}+z_{2}.$ (10)

The equations for electric neutrality (9) and mass balance (10) are explicitly
affected by the presence of a strong base and salts of the conjugate bases
$\ce{HB-}$ and $\ce{B^{2-}}$, _e.g._ $\ce{NaOH}$, $\ce{NaHB}$ and $\ce{Na2B}$
respectively. If the dimensionless concentrations of the strong base and salts
are ${c_{\mathrm{b}}=\ce{[NaOH]}/C^{\standardstate}}$,
${s_{\mathrm{1}}=\ce{[NaHB]}/C^{\standardstate}}$ and
${s_{\mathrm{2}}=\ce{[Na2B]}/C^{\standardstate}}$, the charge and mass balance
equations are modified to

$\displaystyle x+\bar{c}_{\mathrm{b}}$ $\displaystyle=y+z_{1}+2z_{2},$ (11)
$\displaystyle\bar{c}_{\mathrm{a}}$ $\displaystyle=z_{0}+z_{1}+z_{2},$ (12)

with effective concentrations
$\bar{c}_{\mathrm{a}}=c_{\mathrm{a}}+s_{\mathrm{1}}+s_{\mathrm{2}}$ and
$\bar{c}_{\mathrm{b}}=c_{\mathrm{b}}+s_{\mathrm{1}}+2s_{\mathrm{2}}$. These
effective dimensionless variables are related to the effective molar
concentrations by
$\bar{C}_{\mathrm{a}}=C^{\standardstate}\bar{c}_{\mathrm{a}}$, and
$\bar{C}_{\mathrm{b}}=C^{\standardstate}\bar{c}_{\mathrm{b}}$.

The use of $y=1/x$, obtained from equation (8), in the equations (6), (7),
(11) and (12), gives a non-linear system of four equation $\mathcal{S}_{4}$
with four unknowns $x$, $z_{0}$, $z_{1}$ and $z_{2}$.

###  2.2 Mathematical expression for $\ce{[H3O+]}$

Before obtaining the full solution of the non-linear system $\mathcal{S}_{4}$
is useful to analyze the linear subsystem $\mathcal{S}_{3}$ made by equations
(6), (7) and (12). This subsystem can be easily solved to obtain the
concentrations $z_{0}$, $z_{1}$, $z_{2}$, in terms of $x$. The linear system
$\mathcal{S}_{3}$ can be expressed as $\mathsf{K}\bm{z}=\bm{c}$ with
$\mathsf{K}=\mathsf{K}(x)$ given by

$\mathsf{K}=\begin{pmatrix}k_{1}&-x&0\\\ 0&k_{2}&-x\\\ 1&1&1\end{pmatrix},$
(13)

$\bm{z}=\left(z_{0},z_{1},z_{2}\right)^{\intercal}$, and
$\bm{c}=\left(0,0,\bar{c}_{\mathrm{a}}\right)^{\intercal}$. Solving for
$\bm{z}$ gives

$\bm{z}=\frac{\bar{c}_{\mathrm{a}}}{\det{\mathsf{K}}}\begin{pmatrix}x^{2}\\\
k_{1}x\\\ k_{1}k_{2}\end{pmatrix},$ (14)

with $\det{\mathsf{K}}=x^{2}+k_{1}x+k_{1}k_{2}$ as the determinant of
$\mathsf{K}$. It is convenient to write this determinant as
${\det{\mathsf{K}}=(x-\kappa_{1})(x-\kappa_{2})}$ with

$\kappa_{1,2}=\tfrac{k_{1}}{2}\left(-1\pm\sqrt{1-4\kappa}\right),$ (15)

and $\kappa=k_{2}/k_{1}$ as the ratio of the diprotic dissociation constants.
These $\kappa_{1,2}$ are related to the Simms constants [12, 7], $g_{1,2}$, by
$g_{1,2}=-\kappa_{1,2}$.

Given the condition $\kappa\leq 1/4$, _i.e._ $k_{1}\geq 4k_{2}$, the roots
$\kappa_{1,2}$ are both no positive real numbers, $\kappa_{1,2}\leq 0$,
otherwise these roots are a pair of complex conjugate numbers with negative
real part, _i.e._ $\kappa_{1}=\kappa_{2}^{*}$ and
$\mathrm{re}{(\kappa_{1})}=\mathrm{re}{(\kappa_{2})}<0$. The inequality
$k_{1}\geq 4k_{2}$ have been obtained previously by Adams in his analysis of
polyprotic acid dissociations [9].

The solution of $\mathcal{S}_{3}$ gives the concentrations
$\bm{z}=\left(z_{0},z_{1},z_{2}\right)^{\intercal}$ as functions of $x$,
$k_{1}$, $k_{2}$ and $\bar{c}_{\mathrm{a}}$. Although $\bm{z}$ does not depend
explicitly on $\bar{c}_{\mathrm{b}}$, it depends implicitly through $x$. This
dependency is specified by using $y=1/x$ in equation (11), which gives

$x-\tfrac{1}{x}+\bar{c}_{\mathrm{b}}=-\bm{q}\cdot\bm{z},$ (16)

with $\bm{q}=\left(0,-1,-2\right)^{\intercal}$ as the vector of electric
charges of $z_{0}$, $z_{1}$, and $z_{2}$. This equation keeps some similarity
with equation (41) used in the work of Kalka [7]. However, unlike that
article, in this paper closed analytic solutions for $x$ are obtained instead
of graphical or numerical solutions.

Multiplying equation (16) by $x\det{\mathsf{K}}$, and expanding the scalar
product, produces

$\left(x-\kappa_{1}\right)\left(x-\kappa_{2}\right)\left(x^{2}+\bar{c}_{\mathrm{b}}x-1\right)=\bar{c}_{\mathrm{a}}k_{1}x\left(x+2k_{2}\right),$
(17)

which can be written

$\left(x-\kappa_{1}\right)\left(x-\kappa_{2}\right)\left(x-\sigma_{1}\right)\left(x-\sigma_{2}\right)=\bar{c}_{\mathrm{a}}k_{1}x\left(x+2k_{2}\right),$
(18)

with
$\sigma_{1,2}=\tfrac{1}{2}\left(-\bar{c}_{\mathrm{b}}\pm\sqrt{\bar{c}_{\mathrm{b}}^{2}+4}\right)$.
The roots $\sigma_{1,2}$ are both real numbers with $0<\sigma_{1}\leq 1$ and
$\sigma_{2}\leq-1$. The case of $\bar{c}_{\mathrm{b}}=0$ gives
$\sigma_{1,2}=\pm 1$.

Expansion of equation (18) gives $P=0$ with

$P=x^{4}+c_{3}x^{3}+c_{2}x^{2}+c_{1}x+c_{0},$ (19)

and

$\displaystyle c_{3}$ $\displaystyle=\bar{c}_{\mathrm{b}}+k_{1},$ (20)
$\displaystyle c_{2}$
$\displaystyle=-\left(1+k_{1}\left(\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}}-k_{2}\right)\right),$
(21) $\displaystyle c_{1}$
$\displaystyle=-k_{1}\left(1+k_{2}\left(2\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}}\right)\right),$
(22) $\displaystyle c_{0}$ $\displaystyle=-k_{1}k_{2},$ (23)

with $c_{4}=1$.

Before finding the roots of the equation $P=0$ it is helpful to analyze the
nature of its roots, which is studied by considering the 5-tuple of its
coefficients,

$\operatorname{coef}[P]=\left(c_{4},c_{3},c_{2},c_{1},c_{0}\right),$ (24)

and its signs

$\operatorname{sgn}(\operatorname{coef}[P])=\left(\operatorname{sgn}{c_{4}},\operatorname{sgn}{c_{3}},\operatorname{sgn}{c_{2}},\operatorname{sgn}{c_{1}},\operatorname{sgn}{c_{0}}\right).$
(25)

It is straightforward to see that $\operatorname{sgn}{c_{4}}=+$,
$\operatorname{sgn}{c_{3}}=+$, and $\operatorname{sgn}{c_{0}}=-$. The sign of
$c_{2}$ and $c_{1}$ requires a careful analysis. There are four possible
cases: $\left(\operatorname{sgn}{c_{2}},\operatorname{sgn}{c_{1}}\right)$:
$\left(+,+\right)$, $\left(+,-\right)$, $\left(-,+\right)$, and
$\left(-,-\right)$. The 5-tuple $\operatorname{sgn}(\operatorname{coef}[P])$
can have four possible outcomes: $\left(+,+,+,+,-\right)$,
$\left(+,+,+,-,-\right)$, $\left(+,+,-,+,-\right)$, and
$\left(+,+,-,-,-\right)$. These 5-tuples display one or three changes of sign
along the sequence of elements. Descartes’s rule of signs states that the
number of positive roots of a polynomial, $P$, is either equal to the number
of sign changes of $\operatorname{sgn}(\operatorname{coef}[P])$, or is less
than it by an even number. The application of Descartes’ rule to $P$ gives
either one or three positive roots. It can be proved that the polynomial $P$
has only one positive real root by a careful analysis of equation (18). The
left hand side of equation (18) is a fourth degree polynomial $P_{L}(x)$ with
only one positive root $\sigma_{1}$, one negative root $\sigma_{2}$ and two
roots $\kappa_{1,2}$ that could be either negative or a complex conjugate
pair. The right hand side of equation (18) is an upward parabola $P_{R}(x)$
with roots at zero and $-2k_{2}$. The coefficients of the quartic term of
$P_{L}(x)$ and the quadratic term of $P_{R}$ are both positive numbers,
therefore $P_{L}(x)$ and $P_{R}(x)$ must tend to infinity as $x$ goes to
positive or negative infinity. Since the quartic function grows always faster
than the quadratic function, and regarding that $P_{L}(0)=-k_{1}k_{1}<0$, the
polynomials $P_{L}(x)$ and $P_{R}(x)$ must be equal at only one positive $x$.

In the Appendix is shown that using Ferrari’s method [13, 14], the quartic
equation $P=0$ can be written as an associated resolvent cubic equation $R=0$
with

$R=y^{3}-c_{2}y^{2}+\left(c_{1}c_{3}-4c_{0}\right)y+\left(4c_{0}c_{2}-c_{0}c_{3}^{2}-c_{1}^{2}\right).$
(26)

The cubic equation $R=0$ can be solved by Cardano’s method [8], for which the
change of variable $y=\bar{y}+\frac{c_{2}}{3}$ is necessary to obtain a
depressed cubic equation $R_{\mathrm{dc}}=0$ with

$R_{\mathrm{dc}}=\bar{y}^{3}+\bar{p}\bar{y}+\bar{q},$ (27)

where

$\displaystyle\bar{p}$ $\displaystyle=c_{1}c_{3}-\frac{c_{2}^{2}}{3}-4c_{0},$
(28) $\displaystyle\bar{q}$
$\displaystyle=\frac{8c_{0}c_{2}}{3}+\frac{c_{1}c_{2}c_{3}}{3}-\frac{2c_{2}^{3}}{27}-c_{1}^{2}-c_{0}c_{3}^{2},$
(29)

and the discriminant of $P$ is $\Delta=-4\bar{p}^{3}-27\bar{q}^{2}$ [13, 14].

The positive root of $P=0$ is given by three cases depending on the sign of
the $\Delta$, and the sign of the functions $\xi_{1,2}$,

$\xi_{1,2}=-\frac{\bar{q}}{2}\pm\frac{1}{2}\sqrt{-\frac{\Delta}{27}}.$ (30)

The quantities $\Delta$, $\xi_{1}$, and $\xi_{2}$ are functions of the
equilibrium constants, $k_{1}$ and $k_{2}$, and the effective concentrations,
$\bar{c}_{\mathrm{a}}$ and $\bar{c}_{\mathrm{b}}$. Explicitly, the positive
root of $P=0$ is given by

$x=\begin{cases}x_{1},&\Delta>0,\\\ x_{1}&\Delta<0,\;\xi_{1}>0,\;\xi_{2}>0,\\\
x_{3},&\Delta<0,\;\xi_{1}<0,\;\xi_{2}<0,\\\
x_{3},&\Delta<0,\;\xi_{1}>0,\;\xi_{2}<0.\end{cases}$ (31)

The roots $x_{1}$ and $x_{3}$ are:

$\displaystyle x_{1}$
$\displaystyle=\tfrac{1}{2}\left(-\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}-t_{1}\right)+\sqrt{\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}-t_{1}\right)^{2}-2y_{1}+\tfrac{(\bar{c}_{\mathrm{b}}+k_{1})y_{1}+2k_{1}\left(1+(2\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}})k_{2}\right)}{t_{1}}}\right),$
(32) $\displaystyle x_{3}$
$\displaystyle=\tfrac{1}{2}\left(-\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}+t_{1}\right)+\sqrt{\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}+t_{1}\right)^{2}-2y_{1}-\tfrac{(\bar{c}_{\mathrm{b}}+k_{1})y_{1}+2k_{1}\left(1+(2\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}})k_{2}\right)}{t_{1}}}\right),$
(33)

with $y_{1}$ and $t_{1}$:

$\displaystyle t_{1}$
$\displaystyle=\sqrt{1+\tfrac{1}{4}(\bar{c}_{\mathrm{b}}+k_{1})^{2}+k_{1}\left(\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}}-k_{2}\right)+y_{1}},$
(34) $\displaystyle y_{1}$
$\displaystyle=\bar{y}_{1}-\tfrac{1+k_{1}\left(\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}}-k_{2}\right)}{3},$
(35)

and $\bar{y}_{1}$:

$\bar{y}_{1}=\begin{cases}\tfrac{2}{3}\sqrt{1+k_{1}Q_{1}+k_{1}^{2}Q_{2}}\cos{\left(\tfrac{\theta}{3}\right)},&\Delta>0,\\\
\sqrt[3]{|\xi_{1}|}+\sqrt[3]{|\xi_{2}|},&\Delta<0,\;\xi_{1}>0,\;\xi_{2}>0,\\\
-(\sqrt[3]{|\xi_{1}|}+\sqrt[3]{|\xi_{2}|}),&\Delta<0,\;\xi_{1}<0,\;\xi_{2}<0,\\\
\sqrt[3]{|\xi_{1}|}-\sqrt[3]{|\xi_{2}|},&\Delta<0,\;\xi_{1}>0,\;\xi_{2}<0.\end{cases}$
(36)

The functions $\theta$, $Q_{1}$, and $Q_{2}$ are given in the Appendix. The
concentration $x$ for the most common case of diprotic acid, with
$k_{1}>4k_{2}$, obeys $x=x_{1}$ for most of the concentrations.

###  2.3 Strong base titration of dissolutions of diprotic acids and its
buffer mixtures

The titration by a strong base of an acid solution, or an acid buffer
solution, can be analyzed by using the effective concentrations
$\bar{c}_{\mathrm{a}}$ and $\bar{c}_{\mathrm{b}}$. Recall that the effective
concentrations are related to the molar analytical concentrations by
$\bar{c}_{\mathrm{a}}=C_{\mathrm{a}}/C^{\standardstate}$ and
$\bar{c}_{\mathrm{b}}=C_{\mathrm{b}}/C^{\standardstate}$. A buffer solution is
made of a volume $V_{\mathrm{a0}}$ of an acid solution with analytical
concentration $C_{\mathrm{a0}}$, and volumes $V_{10}$ and $V_{20}$ of salt
solutions with analytical concentrations $C_{10}$ and $C_{20}$, respectively.
The total volume of the buffer solution is
$V_{\mathrm{B0}}=V_{\mathrm{a0}}+V_{\mathrm{10}}+V_{\mathrm{20}}$. The case of
and acid solution is obtained by using $V_{10}=0$, and $V_{20}=0$, in the
volume of the buffer: $V_{\mathrm{B0}}=V_{\mathrm{a0}}$. The buffer effective
concentrations are given by

$\displaystyle\bar{c}_{\mathrm{a0}}$
$\displaystyle=\left(c_{\mathrm{a0}}V_{\mathrm{a0}}+s_{\mathrm{10}}V_{\mathrm{10}}+s_{\mathrm{20}}V_{\mathrm{20}}\right)/V_{\mathrm{B0}},\text{
and}$ (37) $\displaystyle\bar{c}_{\mathrm{b0}}$
$\displaystyle=\left(s_{\mathrm{10}}V_{\mathrm{10}}+2s_{\mathrm{20}}V_{\mathrm{20}}\right)/V_{\mathrm{B0}}.$
(38)

These expression for the effective concentrations are obtained from the
analytical concentrations simply by using the scaling factor
$1/C^{\standardstate}$.

A volume $V_{\mathrm{b}}$ of a strong base with analytical concentration
$C_{\mathrm{b0}}$ is added to the buffer of volume $V_{\mathrm{B0}}$ and
effective concentrations $\bar{c}_{\mathrm{a0}}$ and $\bar{c}_{\mathrm{b0}}$.
The addition of the strong base changes the volume of the buffer to
$V_{\mathrm{B}}=V_{\mathrm{B0}}+V_{\mathrm{b}}$, and the effective
concentrations to the titration effective concentrations

$\displaystyle\bar{c}_{\mathrm{a}}$
$\displaystyle=\left(c_{\mathrm{a0}}V_{\mathrm{a0}}+s_{\mathrm{10}}V_{\mathrm{10}}+s_{\mathrm{20}}V_{\mathrm{20}}\right)/V_{\mathrm{B}},\text{
and}$ (39) $\displaystyle\bar{c}_{\mathrm{b}}$
$\displaystyle=\left(c_{\mathrm{b0}}V_{\mathrm{b}}+s_{\mathrm{10}}V_{\mathrm{10}}+2s_{\mathrm{20}}V_{\mathrm{20}}\right)/V_{\mathrm{B}}.$
(40)

The use of the buffer effective concentrations, equations (37) and (38), in
the titration effective concentrations, equations (39) and (40), gives

$\displaystyle\frac{\bar{c}_{\mathrm{a0}}}{\bar{c}_{\mathrm{a}}}-1$
$\displaystyle=\frac{V_{\mathrm{b}}}{V_{\mathrm{B0}}},\text{ and}$ (41)
$\displaystyle\left(\frac{c_{\mathrm{b0}}}{\bar{c}_{\mathrm{b}}}-1\right)\frac{V_{\mathrm{b}}}{V_{\mathrm{B0}}}$
$\displaystyle=1-\frac{\bar{c}_{\mathrm{b0}}}{\bar{c}_{\mathrm{b}}}.$ (42)

The titration effective concentrations can be combined to obtain an equation
for $\bar{c}_{\mathrm{b}}$ in terms of $\bar{c}_{\mathrm{a}}$,
$c_{\mathrm{b0}}$ and the buffer effective concentrations,

$\bar{c}_{\mathrm{b}}=\frac{\bar{c}_{\mathrm{a}}}{\bar{c}_{\mathrm{a0}}}\left(\bar{c}_{\mathrm{b0}}-c_{\mathrm{b0}}\right)+c_{\mathrm{b0}}.$
(43)

This is the equation of a straight line with slope
$\left(\bar{c}_{\mathrm{b0}}-c_{\mathrm{b0}}\right)/\bar{c}_{\mathrm{a0}}$ and
ordinate intercept $c_{\mathrm{b0}}$. The slope of this straight line could be
be negative, zero, or positive, depending on the concentrations of the buffer,
$\bar{c}_{\mathrm{b0}}=s_{\mathrm{10}}+2s_{\mathrm{20}}$, and the titrating
base, $c_{\mathrm{b0}}$.

In addition to the case of the titration of buffer solutions, this equation
can be used for the titration of acid solutions. The case of an acid solution
is obtained by taking $\bar{c}_{\mathrm{b0}}=0$ and
$\bar{c}_{\mathrm{a0}}=c_{\mathrm{a0}}$, to obtain

$\bar{c}_{\mathrm{b}}=c_{\mathrm{b0}}-\left(\frac{c_{\mathrm{b0}}}{c_{\mathrm{a0}}}\right)\bar{c}_{\mathrm{a}},$
(44)

where $\bar{c}_{\mathrm{a}}=c_{\mathrm{a}}$ and
${\bar{c}_{\mathrm{b}}=c_{\mathrm{b}}}$. This is the equation of a straight
line with slope ${-c_{\mathrm{b0}}/c_{\mathrm{a0}}}$ and ordinate intercept
$c_{\mathrm{b0}}$.

Equations (43) and (44) describe the titrations as straight line paths on the
$(\bar{c}_{\mathrm{a}},\bar{c}_{\mathrm{b}})$ plane. The addition of the base,
with concentration $c_{\mathrm{b0}}$, increases $\bar{c}_{\mathrm{b}}$ along a
straight line, from $\bar{c}_{\mathrm{b0}}$ to $c_{\mathrm{b0}}$, meanwhile
decreases $\bar{c}_{\mathrm{a}}$, from $\bar{c}_{\mathrm{a0}}$ to zero. The
use of the contours of constant pH on the
$(\bar{c}_{\mathrm{a}},\bar{c}_{\mathrm{b}})$ plane and the trajectory
described by equation (43), or (44), give the full description of the
titration experiments. However, it is more practical to describe the titration
experiment as function of the pH, instead of the added strong base. For this,
equation (16) is used, which relates $x$, the pH, with the effective
concentrations $\bar{c}_{\mathrm{a}}$ and $\bar{c}_{\mathrm{b}}$. After some
rearrangement, equation (16) gives

$x-\tfrac{1}{x}+\bar{c}_{\mathrm{b}}=\frac{\bar{c}_{\mathrm{a}}k_{1}\left(x+2k_{2}\right)}{x^{2}+k_{1}x+k_{1}k_{2}}.$
(45)

The use of the ratio between the effective dimensionless concentrations,
$\bar{n}={\bar{c}_{\mathrm{b}}}/{\bar{c}_{\mathrm{a}}}$, in equation (45)
gives

$\bar{n}=\frac{k_{1}\left(x+2k_{2}\right)}{x^{2}+k_{1}x+k_{1}k_{2}}+\frac{1-x^{2}}{x\bar{c}_{\mathrm{a}}}.$
(46)

This equation has been reported by Kalka previously [7], and works well in the
case of $\bar{c}_{\mathrm{a}}$ constant. However, in the titration experiment
the concentrations $\bar{c}_{\mathrm{a}}$ and $\bar{c}_{\mathrm{b}}$ are not
constants. The last term of the right hand side of this equation depends on
$\bar{c}_{\mathrm{a}}$. This dependence can be eliminated using equation (43).
For this purpose, equation (43) must be written in terms of $\bar{n}$,

$\bar{n}=\bar{n}_{0}+c_{\mathrm{b0}}\left(\frac{1}{\bar{c}_{\mathrm{a}}}-\frac{1}{\bar{c}_{\mathrm{a0}}}\right),$
(47)

with $\bar{n}_{0}={\bar{c}_{\mathrm{b0}}}/{\bar{c}_{\mathrm{a0}}}$. Algebraic
manipulation of equation (47) gives

$\frac{1}{\bar{c}_{\mathrm{a}}}=\frac{\bar{n}-\bar{n}_{0}}{c_{\mathrm{b0}}}+\frac{1}{\bar{c}_{\mathrm{a0}}},$
(48)

which can be used in the last term of equation (46) to obtain an expression
for $\bar{n}$ in terms of the pH

$\bar{n}=\left(\bar{n}_{0}-\frac{c_{\mathrm{b0}}}{\bar{c}_{\mathrm{a0}}}\right)\frac{P_{\mathrm{a0}}}{P_{\mathrm{b0}}},$
(49)

where $\bar{n}\geq\bar{n}_{0}$ and

$\displaystyle P_{\mathrm{a0}}$
$\displaystyle=P\left(x=10^{7-\mathrm{pH}},\bar{c}_{\mathrm{a}}=\frac{\bar{c}_{\mathrm{a0}}c_{\mathrm{b0}}}{c_{\mathrm{b0}}-\bar{c}_{\mathrm{b0}}},\bar{c}_{\mathrm{b}}=0\right),$
(50) $\displaystyle P_{\mathrm{b0}}$
$\displaystyle=P\left(x=10^{7-\mathrm{pH}},\bar{c}_{\mathrm{a}}=0,\bar{c}_{\mathrm{b}}=c_{\mathrm{b0}}\right),$
(51)

with $P$ the polynomial given by equation (19). The case of the acid titration
is given by considering $\bar{c}_{\mathrm{b0}}=0$ in equations (49)–(51).

The equivalence points for a diprotic acid occur at $\bar{n}=1,2$ in equation
(49). The first equivalence point is given when the acid and base
concentrations are equal, $\bar{n}=1$, the second equivalence point happens
when the base concentration doubles the acid concentration, $\bar{n}=2$. Since
$\bar{n}(\mathrm{pH})$ must be a monotonically growing function of the pH, it
must fulfill the condition $n^{\prime}(\mathrm{pH})>0$.

##  3 Results and discussion

The validity of equations (31)–(36) has been tested by calculating the pH, at
different effective concentrations $\bar{c}_{\mathrm{a}}$ and
$\bar{c}_{\mathrm{b}}$, of diprotic acids with reported values of p$K_{1}$ and
p$K_{2}$ for 180 diprotic acids [10] and comparing them with the numerical
solution. The average absolute error in the pH, $\epsilon_{\mathrm{pH}}$, at
millimolar concentrations $\bar{c}_{\mathrm{a}}$ and $\bar{c}_{\mathrm{b}}$,
is $\epsilon_{\mathrm{pH}}\lesssim 10^{-5}$ pH units with a standard deviation
$\lesssim 10^{-4}$ pH units. The small error in the calculated pH is caused
mainly by the diprotic acids with
$\mathrm{p}K_{2}-\mathrm{p}K_{1}<-\log_{10}4$.

###  3.1 Aqueous dissolution of a diprotic acid

The case of a diprotic acid is given by using the conditions
$\bar{c}_{\mathrm{b}}=0$ and $\bar{c}_{\mathrm{a}}=c_{\mathrm{a}}$ in
equations (31)–(36). The condition $k_{1}\geq 4k_{2}$, _i.e._ ,

$\mathrm{p}K_{2}-\mathrm{p}K_{1}\geq{\mathrm{log}}_{10}4\approx 0.602,$ (52)

makes the discriminant of $P$ a positive number, $\Delta>0$. This is the
condition for a quartic equation with four distinct real roots. An inspection
of the values of $\mathrm{p}K_{1}$ and $\mathrm{p}K_{2}$ of tabulated diprotic
weak acids indicates that the condition (52) is fulfilled by many diprotic
weak acids [10, 9, 11]. Figure 1 displays, on the p$K_{1}$–p$K_{2}$ plane, the
region given by the condition (52) as the light blue region above the blue
line. It is clear in the Figure that most of the diprotic weak acids (open
black circles) fulfill this condition, however a simple visual inspection of
Figure 1 shows that there are weak diprotic acids that fulfill the condition
${\mathrm{p}K_{2}-\mathrm{p}K_{1}\leq-\log_{10}4}$ (light red region). This
condition can be expressed in terms of the acid constants as $k_{2}/k_{1}\geq
4$, _i.e._ $K_{2}/K_{1}\geq 4$. Diprotic acids in the light red region have
$\mathrm{p}K_{1}>\mathrm{p}K_{2}$, examples of these are several nitrogenous
organic compounds as Piperazine, Quinine, Phenylbiguanide, $L$-Nicotine,
$p$-Benzidine, Sulfamethazine, $m$-Phenylenediamine, $p$-Phenylenediamine,
1,2-Propanediamine, 1,3-Propanediamine, 1,4-Butanediamine, 1,6-Hexanediamine,
1,8-Octanediamine, cis-2,5-Dimethylpiperazine, trans-1,2-Cyclohexanediamine,
cis-1,2-Cyclohexanediamine, and the alcohol 1,3-Diamino-2-propanol [10].

Figure 1: p$K_{1}$–p$K_{2}$ plane for a set of diprotic acids in aqueous
solution [10]. The light blue region is given by the condition
$\mathrm{p}K_{2}-\mathrm{p}K_{1}\geq\log_{10}4$, the light red region is given
by the condition $\mathrm{p}K_{2}-\mathrm{p}K_{1}\leq-\log_{10}4$

###  3.2 Concentrations $\ce{[H2B]}$, $\ce{[HB^{-}]}$, and $\ce{[B^{2-}]}$ in
terms of $\ce{[H3O+]}$

Equation (18) can be written to give an expression for $\det{\mathsf{K}}$,

$\det{\mathsf{K}}=\frac{\bar{c}_{\mathrm{a}}k_{1}x\left(x+2k_{2}\right)}{\left(x-\sigma_{1}\right)\left(x-\sigma_{2}\right)}.$
(53)

This equation can be used in equation (14), to obtain

$\displaystyle z_{0}$
$\displaystyle=\frac{x\left(x-\sigma_{1}\right)\left(x-\sigma_{2}\right)}{k_{1}(x+2k_{2})},$
(54) $\displaystyle z_{1}$
$\displaystyle=\frac{\left(x-\sigma_{1}\right)\left(x-\sigma_{2}\right)}{x+2k_{2}},$
(55) $\displaystyle z_{2}$
$\displaystyle=\frac{k_{2}\left(x-\sigma_{1}\right)\left(x-\sigma_{2}\right)}{x\left(x+2k_{2}\right)}.$
(56)

These concentrations are constrained to be positive numbers. Since
$0<\sigma_{1}\leq 1$ and $\sigma_{2}\leq-1$, it is necessary to have
$x>\sigma_{1}$.

It is also possible to have parametric dependence on $\bar{c}_{\mathrm{a}}$,
$\bar{c}_{\mathrm{b}}$ and $k_{2}$ by using equation (12), which gives as
result equations (55), (56), and

$z_{0}=\bar{c}_{\mathrm{a}}-\frac{\left(x+k_{2}\right)\left(x-\sigma_{1}\right)\left(x-\sigma_{2}\right)}{x\left(x+2k_{2}\right)}$
(57)

instead of (54). The case of a dissolution of the diprotic acid gives
$\sigma_{1}=1$ and $\sigma_{2}=-1$, with $\bar{c}_{\mathrm{b}}=0$ and
$\bar{c}_{\mathrm{a}}=c_{\mathrm{a}}$.

It is convenient to employ logarithmic scaling to describe concentrations and
equilibrium constants of highly diluted solutions and weak acids. The
p-function of $\mathrm{Q}$ is defined as

$\begin{split}\mathrm{pQ}&=-\log_{10}a_{\mathrm{Q}}\\\
&=-\log_{10}\frac{\gamma_{\mathrm{Q}}\ce{[Q]}}{C^{\circ}},\end{split}$ (58)

with $a_{\mathrm{Q}}$ and $\gamma_{\mathrm{Q}}$ as the activity and the
activity coefficient of $\mathrm{Q}$ respectively [4]. Since equilibrium
constants are dimensionless, it is possible to define the p-function of $K$ as
$\mathrm{p}K=-\log_{10}K$ [6].

The case of weak acids, and low concentrations, allows to use the ideal
solution approximation, $\gamma_{\mathrm{Q}}\approx 1$, hence the pH is given
by

$\begin{split}\mathrm{pH}&\approx-\log_{10}\frac{\ce{[H3O+]}}{C^{\circ}}\\\
&\approx-\log_{10}\frac{C^{\standardstate}x}{C^{\circ}}\\\ &\approx
7-\log_{10}x.\end{split}$ (59)

The p-functions for $\ce{[H_{2}B]}$, $\ce{[HB^{-}]}$ and $\ce{[B^{2-}]}$ are
given by $\mathrm{pH_{2}B}\approx 7-\log_{10}z_{0}$, $\mathrm{pHB^{-}}\approx
7-\log_{10}z_{1}$, and $\mathrm{pB^{2-}}\approx 7-\log_{10}z_{2}$,
respectively. The p-function $\mathrm{pH_{2}B}$ can be expressed in two ways,
either using $z_{0}$ from equation (54) to obtain
$\mathrm{pH_{2}B}(k_{1},k_{2})$, or using $z_{0}$ from equation (57) to get
$\mathrm{pH_{2}B}(c_{\mathrm{a}},k_{2})$.

Figure 2 displays the behaviour of pH2B, pHB- and pB2- as functions of the pH
for different concentrations $c_{\mathrm{a}}=2^{n}$ ($n=0,1,\dots,23$) of
oxalic acid, $\ce{H2C2O4}$, which has $K_{1}=5.62\times 10^{-2}$ and
$K_{2}=1.54\times 10^{-4}$ [10]. In this Figure the intersections between the
$\mathrm{pH_{2}B}(k_{1},k_{2})$ curve (red) and the
$\mathrm{pH_{2}B}(c_{\mathrm{a}},k_{2})$ curves (pink) are shown as labeled
black points. These intersections give the pH for the different concentrations
$c_{\mathrm{a}}$.

Figure 2: Behaviour of pH2B (red and pink), pHB- (green) and pB2- (blue) as
functions of the pH for oxalic acid, $\ce{H2C2O4}$. Labeled black dots
indicate the intersection between pH2B$(k_{1},k_{2})$ (red) and
pH2B$(c_{\mathrm{a}},k_{2})$ (pink) at different concentrations,
$c_{\mathrm{a}}=2^{n}$, $n=0,1,2,\dots,23$. Recall that that
$C_{\mathrm{a}}=C^{\standardstate}c_{\mathrm{a}}$, hence for $n=0$,
$C_{\mathrm{a}}=10^{-7}\,\mathrm{M}$, and for $n=23$, $C_{\mathrm{a}}\approx
0.84\,\mathrm{M}$.

###  3.3 Use of physical constraints of the system to obtain approximate
expressions for $\ce{[H3O+]}$

For the diprotic acid, the combined use of equation (57) and the condition
$z_{0}>0$, gives the inequality $P_{z_{0}}<0$ with $P_{z_{0}}$ given by the
monic cubic polynomial

$P_{z_{0}}=x^{3}+\left(k_{2}-c_{\mathrm{a}}\right)x^{2}-\left(1+2c_{\mathrm{a}}k_{2}\right)x-k_{2}.$
(60)

Although this polynomial goes to infinity as $x$ goes to infinity, there are
values of $x$ for which the inequality $P_{z_{0}}<0$ is satisfied. This can be
seen by analyzing the 4-tuple of $P_{z_{0}}$ coefficients,

$\begin{split}\operatorname{coef}[P_{z_{0}}]&=\left(a_{3},a_{2},a_{1},a_{0}\right)\\\
&=\left(1,k_{2}-c_{\mathrm{a}},-(1+2c_{\mathrm{a}}k_{2}),-k_{2}\right).\end{split}$
(61)

The signs of (61) are given by

$\begin{split}\operatorname{sgn}[P_{z_{0}}]&=\left(\operatorname{sgn}a_{3},\operatorname{sgn}a_{2},\operatorname{sgn}a_{1},\operatorname{sgn}a_{0}\right)\\\
&=\left(+,\pm,-,-\right),\end{split}$ (62)

Regardless of the value of $\operatorname{sgn}a_{2}$, there is only one change
of sign in (62), from positive to negative, for this case Descartes rule of
signs gives that $P_{z_{0}}$ must have only one positive root. This positive
root is a function of $k_{2}$ and $c_{\mathrm{a}}$, and gives the upper bound
of $x$. Using the method of Caicedo, et al. the upper bound of $x$ is given by

$x_{\mathrm{ub}}=\tfrac{2}{3}\sqrt{(k_{2}-c_{\mathrm{a}})^{2}+6c_{\mathrm{a}}k_{2}+3}\cos{\left({\theta_{z_{0}}}/{3}\right)}-\frac{k_{2}-c_{\mathrm{a}}}{3},$
(63)

with

$\displaystyle p_{z_{0}}$
$\displaystyle=-\tfrac{1}{3}(k_{2}-c_{\mathrm{a}})^{2}-2c_{\mathrm{a}}k_{2}-1,$
(64) $\displaystyle q_{z_{0}}$
$\displaystyle=\tfrac{2}{27}(k_{2}-c_{\mathrm{a}})^{3}+\tfrac{1}{3}(k_{2}-c_{\mathrm{a}})(1+2c_{\mathrm{a}}k_{2})-k_{2},$
(65) $\displaystyle\Delta_{z_{0}}$
$\displaystyle=-4p_{z_{0}}^{3}-27q_{z_{0}}^{2},$ (66)
$\displaystyle\theta_{z_{0}}$
$\displaystyle=\arctan\left(-\frac{q_{z_{0}}}{2},\frac{\sqrt{\Delta_{z_{0}}}}{6\sqrt{3}}\right).$
(67)

The use of Wolfram Mathematica allows to prove that $\Delta_{z_{0}}>0$ for
$c_{\mathrm{a}}>0$ and $k_{2}>0$. Since $\Delta_{z_{0}}>0$, equation (67)
gives the $0<\theta<\pi$, hence $\cos{(\theta/3)}\geq 1/2$. Furthermore, the
use of Wolfram Mathematica allows to find that $\lim_{c_{\mathrm{a}}\to
0}x_{\mathrm{ub}}=1$.

It was shown above that the dissociation constants of many diprotic acids are
constraint by the condition ${k_{1}\geq 4k_{2}}$. The use of equations (6) and
$\eqref{eq:ka2}$ in the inequality ${k_{1}\geq 4k_{2}}$ leads to the
constraint ${z_{1}^{2}\geq 4z_{0}z_{2}}$ between the concentrations of the
acid and its conjugate bases. The use of equations (55), (56) and (57) in the
inequality ${z_{1}^{2}\geq 4z_{0}z_{2}}$ gives the inequality $P_{z}>0$ with

$P_{z}=x^{3}+2k_{2}x^{2}-(1+4c_{\mathrm{a}}k_{2})x-2k_{2}.$ (68)

Figure 3: Upper (blue) and lower (green) bounds of the pH for oxalic acid
(left) and 1,5-Pentanediamine (right) as functions of the molar concentration
$C_{\mathrm{a}}$. The orange curve represents the exact pH.

This polynomial $P_{z}$ is the same polynomial for a monoprotic weak acid with
$k_{\mathrm{a}}=2k_{2}$. Since
$\operatorname{sgn}[P_{z}]=\left(+,+,-,-\right)$, Descartes rule of signs
indicates that $P_{z}=0$ has only one positive root. Using the method of
Caicedo et al., this root is given by

$x_{\mathrm{lb}}=\tfrac{2}{3}\left(\sqrt{4k_{2}^{2}+3c_{\mathrm{a}}k_{2}+3}\cos{(\theta_{z}/3)}-k_{2}\right),$
(69)

with

$\displaystyle p_{z}$
$\displaystyle=-\tfrac{4k_{2}^{2}}{3}-c_{\mathrm{a}}k_{2}-1,$ (70)
$\displaystyle q_{z}$
$\displaystyle=\tfrac{16k_{2}^{3}}{27}+\tfrac{2c_{\mathrm{a}}k_{2}^{2}}{3}-\tfrac{k_{2}}{3},$
(71) $\displaystyle\Delta_{z}$ $\displaystyle=-4p_{z}^{3}-27q_{z}^{2},$ (72)
$\displaystyle\theta_{z}$
$\displaystyle=\arctan\left(-\frac{q_{z}}{2},\frac{\sqrt{\Delta_{z}}}{6\sqrt{3}}\right).$
(73)

The discriminant $\Delta_{z}$ is a positive quantity, in fact, by using
Wolfram Mathematica it is shown that $\Delta_{z}\geq 4$. Furthermore using the
same software it is shown that $\lim_{c_{\mathrm{a}}\to 0,k_{2}\to
0}{x_{\mathrm{lb}}}=1$, and that $\lim_{c_{\mathrm{a}}\to
0,k_{2}\to\infty}{x_{\mathrm{lb}}}=1/\sqrt{2}$.

The lower and upper bounds to the pH are obtained by
${7-\log_{10}{x_{\mathrm{ub}}}}$ and ${7-\log_{10}{x_{\mathrm{lb}}}}$
respectively. Figure 3 displays the lower and upper bounds to the pH as a
function of the molar concentration $C_{\mathrm{a}}$ for oxalic acid (left)
and 1,5-Pentanediamine (right). In the case of oxalic acid, the exact pH and
the lower pH bound (green) are nearly identical for concentrations
$C_{\mathrm{a}}<10^{-2}\,\mathrm{M}$. On the other hand, for the compound
1,5-Pentanediamine the exact pH and the upper pH bound (blue) are nearly
identical for concentrations $C_{\mathrm{a}}<10^{-4}\,\mathrm{M}$. It is
interesting to notice from both panels of Figure 3 that the upper bound pH
curve (blue) overstimate the pH for a constant difference for concentrations
greater than $10^{-2}\,\mathrm{M}$.

###  3.4 Analysis of the dependence of the $\mathrm{pH}$ on p$K_{2}$

Figure 4: Contours of constant pH (red lines) for different concentrations
$C_{\mathrm{a}}$: (a) $0.1\,\mathrm{M}$; (b) $0.01\,\mathrm{M}$; (c)
$10^{-3}\,\mathrm{M}$; (d) $10^{-6}\,\mathrm{M}$. Contours of
$\delta\mathrm{pH}_{2}$ are shown in shading colors from blue to yellow: dark
blue 0.01, grey 0.02, dark orange 0.06, orange 0.1, light orange 0.14, dark
yellow 0.15, and yellow 0.16. The light blue region has
$\delta\mathrm{pH}_{2}<0.01$. The open circle markers are the values of
$\left(\mathrm{p}K_{1},\mathrm{p}K_{2}\right)$ for different weak diprotic
acids [10].

The pH is calculated by $\mathrm{pH}=7-\log_{10}x$ with $x$ given by equation
(31). The partial derivatives

$\displaystyle\delta\mathrm{pH}_{1}$
$\displaystyle=\left(\frac{\partial\mathrm{pH}}{\partial\mathrm{p}K_{1}}\right)_{C_{\mathrm{a}},\mathrm{p}K_{2}},$
(74) $\displaystyle\delta\mathrm{pH}_{2}$
$\displaystyle=\left(\frac{\partial\mathrm{pH}}{\partial\mathrm{p}K_{2}}\right)_{C_{\mathrm{a}},\mathrm{p}K_{1}},$
(75)

measure how much the pH depends on p$K_{1}$ or p$K_{2}$, respectively.

Figure 4 displays contours of constant $\delta\mathrm{pH}_{2}$ as function of
p$K_{1}$ and p$K_{2}$ for different acid concentrations, $C_{\mathrm{a}}$: (a)
$0.1\,\mathrm{M}$; (b) $0.01\,\mathrm{M}$; (c) $10^{-3}\,\mathrm{M}$; (d)
$10^{-6}\,\mathrm{M}$. Only the acids that fulfill the condition
$\mathrm{p}K_{2}\geq\mathrm{p}K_{1}+\log_{10}4$ are shown. This Figure also
shows contours of constant pH as red curves with their respective call-outs
indicating the value of the pH. In all the panels, the contours of the
derivative $\delta\mathrm{pH}_{2}$ are shown with contour shading from blue to
yellow; the dark blue has $0.01<\delta\mathrm{pH}_{2}<0.02$, the grey has
$0.02<\delta\mathrm{pH}_{2}<0.06$, the dark orange
$0.06<\delta\mathrm{pH}_{2}<0.1$, the orange $0.1<\delta\mathrm{pH}_{2}<0.14$,
the light orange $0.14<\delta\mathrm{pH}_{2}<0.15$, the dark yellow
$0.14<\delta\mathrm{pH}_{2}<0.15$, and the yellow
$0.16<\delta\mathrm{pH}_{2}$. The open circle markers in all the panels of
Figure 4 are the values of $\left(\mathrm{p}K_{1},\mathrm{p}K_{2}\right)$ for
different diprotic weak acids [10]. The maximum value of
$\delta\mathrm{pH}_{2}$ is obtained by evaluating $\delta\mathrm{pH}_{2}$
along the line ${\mathrm{p}K_{2}=\mathrm{p}K_{1}+\log_{10}4}$. By doing this a
function $\delta\mathrm{p}K_{2}(\mathrm{p}K_{1},C_{\mathrm{a}})$ is obtained.
The use of the function $\mathsf{NMaximize}$ of Wolfram Mathematica gives that
${\max{\left(\delta\mathrm{p}K_{2}(\mathrm{p}K_{1},C_{\mathrm{a}})\right)}\approx
0.17153}$ regardless the values of p$K_{1}$ and $C_{\mathrm{a}}$.

Panel (a) of Figure 4 shows that, for $C_{\mathrm{a}}=0.1\,\mathrm{M}$,
p$K_{2}$ has weak influence on the pH for p$K_{1}>4$. This is evident by
observing that the contours of constant pH are practically vertical lines for
$\mathrm{pH}>2.5$ and also by the fact that $\delta\mathrm{pH}_{2}<0.01$ for
the same values of the pH. In the same panel can be observed that the
strongest influence of p$K_{2}$ on the pH, _i.e._
$0.1<\delta\mathrm{p}K_{2}\lesssim 0.17153$, is seen in the regions with
orange and yellow contour shading, and pH$<1.5$. The pH contours in this
region are curved instead of straight vertical lines. Panel (a) shows that the
approximation of considering the pH independent p$K_{2}$ is very good for all
the acids at a concentration $C_{\mathrm{a}}=0.1\,\mathrm{M}$. The highest
observed value of $\delta\mathrm{pH}_{2}$ is about 0.17153 units of pH by 1
unit of p$K_{2}$. This value of $\delta\mathrm{pH}_{2}$ indicates that a
change on 0.5 units of p$K_{2}$ would produce at most a change of 0.08 units
on the pH. This change in the pH is sufficiently small to be within the
experimental error, hence, the heuristic approximation of considering the pH
dependent only on p$K_{1}$ is an good approximation at relatively high
concentrations of the acid.

Panels (b) and (c) of Figure 4 show similarities in shape with panel (a). It
can be seen that for these concentrations the pH is insensitive to the value
of p$K_{2}$ for p$K_{1}>5$ and p$K_{1}>6$, for concentrations
$C_{\mathrm{a}}=10^{-2}$ (b) and $C_{\mathrm{a}}=10^{-3}$ (c) respectively.
Regarding the pH contours, they are insensitive to the p$K_{2}$ value for
pH$>3.5$ and pH$>4.5$ for concentrations $C_{\mathrm{a}}=10^{-2}$ (b) and
$C_{\mathrm{a}}=10^{-3}$ (c) respectively. Panel (d),
$C_{\mathrm{a}}=10^{-6}\,\mathrm{M}$, shows evident differences with respect
to panels (a) to (c). At this low acid concentration the region with
$5.73<\mathrm{pH}<6.5$ displays strong sensitivity to the value of p$K_{2}$,
with the strongest effect for pH$\approx 5.9$ and p$K_{1}\approx 6$.

Panels (a) to (c) of Figure 4 show that the pH contour with
$\mathrm{pH}=-\log_{10}C_{\mathrm{a}}$ is a straight line with negative slope.
This line is a boundary between two regions: one with pH contours that are
asymptotically independent of p$K_{2}$ (vertical lines) and other region with
pH contours that are asymptotically independent of p$K_{1}$ (horizontal
lines). Although panel (d) does not display this straight line with negative
slope, it is clear that there are also regions with asymptotic independence on
p$K_{1}$ and p$K_{2}$.

###  3.5 Strong base titration of diprotic acids and its buffer mixtures

Figure 5: Contours of constant pH on the concentrations plane
$(C_{\mathrm{a}},C_{\mathrm{b}})$ for diprotic acids. Left: maleic acid,
p$K_{1}=1.92$, p$K_{2}=6.23$; Right: 1,8-Octanediamine p$K_{1}=11$,
p$K_{2}=10.1$. The red region is given for $\Delta_{\mathrm{dc}}>0$, the green
region for $\xi_{1}<0$ and $\xi_{2}<0$, and the blue $\xi_{1}>0$ and
$\xi_{2}<0$. The titration line is shown in dashed on both figures. The
equivalence points $n=1,2$ are shown as cyan open circles. The half
equivalence points $n=1/2,3/2,5/2$ are shown as magenta open circles.

Figure 5 displays the contours of constant pH for maleic acid (left panel) and
1,8-Octanediamine (right panel). Maleic acid has p$K_{1}=1.92$ and
p$K_{2}=6.23$, that is p$K_{2}-$p$K_{1}\geq\log_{10}4$, meanwhile the compound
1,8-Octanediamine has p$K_{1}=11$ and p$K_{2}=10.1$, and
p$K_{2}-$p$K_{1}\leq-\log_{10}4$. Maleic acid is an example of a pH given
uniquely by equation (32) with $\bar{y}_{1}$ given by the first case of
equation (36), therefore maleic acid displays only one region in the
$(C_{\mathrm{a}},C_{\mathrm{b}})$ plane. On the other hand the compound
1,8-Octanediamine displays three regions on the
$(C_{\mathrm{a}},C_{\mathrm{b}})$ plane: the red region is given by using
equation (32) with $\bar{y}_{1}$ given by the first case of (36), the green
region given by equation (33) with $\bar{y}_{1}$ given by the second case of
equation (36), and the blue region is given by (33) with $\bar{y}_{1}$ given
by the third case of equation (36).

Figure 6: Titration curves for 10$\,\mathrm{ml}$ of diprotic acids at
concentration $C_{\mathrm{a}}=1\,\mathrm{mM}$ using a volume $V_{\mathrm{b}}$
of strong base with concentration $C_{\mathrm{b}}=1\,\mathrm{mM}$. (Left)
Maleic acid, p$K_{1}=1.92$, p$K_{2}=6.23$; (Right) 1,8-Octanediamine
p$K_{1}=11$, p$K_{2}=10.1$. The red, green and blue curves are calculated
using the first, third, and fourth case of equation (36), respectively.

Figure 6 shows the titrations curves obtained by adding a volume
$V_{\mathrm{b}}$ of a strong base with concentration
$C_{\mathrm{b0}}=1\,\mathrm{mM}$ to a volume $V_{\mathrm{a0}}=10\,\mathrm{ml}$
of solution $C_{\mathrm{a0}}=1\,\mathrm{mM}$ of maleic acid (left) and
1,8-Octanediamine (right). These titration curves are given by the pH along
the dashed lines of Figure 5 for the case $C_{\mathrm{a0}}=1\,\mathrm{mM}$ and
$C_{\mathrm{b0}}=1\,\mathrm{mM}$.

The left panel of Figure 6 displays the typical equivalence points for
diprotic acids with two equivalence points. The first equivalence point occurs
at $V_{\mathrm{b}}\approx 10\,\mathrm{ml}$ with $\mathrm{pH}\approx 5$, the
second equivalence point occurs at $V_{\mathrm{b}}\approx 20\,\mathrm{ml}$
with $\mathrm{pH}\approx 8$. In contrast, the right panel of Figure 6 does not
display equivalence points. The titration curve for 1,8-Octanediamine is made
by joining three different titration curves: red, green and blue (from left to
right). The initial solution of 1,8-Octanediamine has pH slightly below 7, as
the base is added the pH grows rapidly reaching a pH above 10. This behaviour
is described by the red curve of Figure 6 (right). As the volume of added base
increases, the pH grows from 10 to approximately 10.5, following the green
curve. Finally, for $V_{\mathrm{b}}>25\,\mathrm{ml}$ the titration experiment
follows the blue curve reaching a final pH slightly above 10.5 at
$V_{\mathrm{b}}=40\,\mathrm{ml}$.

Figure 7: Titration functions $\bar{n}(\mathrm{pH})$, and
$\bar{n}_{\mathrm{K}}$, obtained from equations (49) (green and red) and (46)
(cyan and magenta), respectively. The first and second equivalence points
occur at $\bar{n},\bar{n}_{K}=1,2$, the half-equivalence points occur at
$\bar{n},\bar{n}_{\mathrm{K}}=1/2,3/2,5/2$. Maleic acid (red and magenta)
displays the typical titration curve of a diprotic acid, 1,8-Octanediamine
(green and cyan) does not diaplay equivalence points.

The use of equation (49) allows to obtain the pH at the equivalence points for
acid and buffer solutions. The equivalence points $\bar{n}=1,2$, and the half-
equivalence points ${\bar{n}=1/2,3/2}$, as functions of the pH, are displayed
in Figure 7 for maleic acid, p$K_{1}=1.92$, p$K_{2}=6.23$, and for
1,8-Octanediamine, p$K_{1}=11$, p$K_{2}=10.1$. It is seen in this Figure that
the equivalence point for $\bar{n},\bar{n}_{\mathrm{K}}=3/2$ are the same for
maleic acid, but not for 1,8-Octanediamine. This Figure also shows that the
use of equation (46) gives a wrong pH at the equivalence points $\bar{n}=1,2$.
The first equivalence point is shifted to more acid pHs meanwhile the second
equivalence point is shifted to basic pHs.

###  3.6 pH stability of buffer solutions

In the titration experiment a volume $V_{\mathrm{b}}$ of a base solution, with
concentration $c_{\mathrm{b0}}$, is added to a volume $V_{\mathrm{B0}}$ of a
buffer solution with concentrations $\bar{c}_{\mathrm{a0}}$ and
$\bar{c}_{\mathrm{b0}}$. The concentrations $\bar{c}_{\mathrm{a}}$ and
$\bar{c}_{\mathrm{b}}$ as functions of $V_{\mathrm{b}}$ are given by equations
(39) and (40). Although equation (49) can be used to analyze the stability of
a buffer solution, it is more convenient to use the parametric curve

$\beta(V_{\mathrm{b}})=\left\\{\mathrm{pH}_{\mathrm{acid}}\left(V_{\mathrm{b}}\right),\mathrm{pH}_{\mathrm{buffer}}\left(V_{\mathrm{b}}\right)\right\\},$
(76)

where $\mathrm{pH}_{\mathrm{acid}}\left(V_{\mathrm{b}}\right)$ is the pH of
the acid as function of added base, and
$\mathrm{pH}_{\mathrm{buffer}}\left(V_{\mathrm{b}}\right)$ is the pH of the
buffer as function of added base. Figure 8 displays $\beta(V_{\mathrm{b}})$
for acid and buffer solutions prepared with the same number of moles of the
acid, $C_{\mathrm{a}}V_{\mathrm{a0}}=7.5\times 10^{-3}$ moles, and titrated
with the same strong base, $C_{\mathrm{b0}}=7.5\,\mathrm{mM}$. The buffer
solutions of panel (a) are prepared by adding $C_{\mathrm{10}}V_{10}=2.5\times
10^{-3}$ moles of the monobasic salt only; The buffer solutions of panel (b)
are prepared by adding $C_{\mathrm{20}}V_{20}=2.5\times 10^{-3}$ moles of the
dibasic salt only; The buffer solutions of panel (c) are prepared by adding
$C_{\mathrm{10}}V_{10}=2.5\times 10^{-3}$ moles of the monobasic salt, and
$C_{\mathrm{20}}V_{20}=2.5\times 10^{-3}$ moles of the dibasic salt. The three
panels show four curves for different acids all with p$K_{1}=1$ and p$K_{2}=1$
(red), p$K_{2}=4$ (green), p$K_{2}=6$ (blue), and p$K_{2}=8$ (cyan).

Figure 8: $\beta(V_{\mathrm{b}})$ curves for buffer solutions of different
acids prepared with: (a) only monobasic salt, (b) only dibasic salt, and (c)
both monobasic and dibasic salts. All the acids have p$K_{1}=1$ and:
p$K_{2}=1$ (red), p$K_{2}=4$ (green), p$K_{2}=6$ (blue), and p$K_{2}=8$
(cyan).

The pH stability of a buffer solution is given when the
$\beta(V_{\mathrm{b}})$ curve is horizontal, _i.e._ regardless the change in
the pH of the acid solution, the pH of the buffer remains stable. The red
curve of panels (a) and (c) of Figure 8 display the best pH stability. These
red curves are produced by acids with p$K_{1}=1$ and p$K_{2}=1$. The red
$\beta(V_{\mathrm{b}})$ curves of panels (a) and (c) display buffer stability
at $\mathrm{pH}_{\mathrm{buffer}}\approx 3$ and
$\mathrm{pH}_{\mathrm{acid}}>4$. The green $\beta(V_{\mathrm{b}})$ curves of
panels (a) and (c) display buffer stability at
$\mathrm{pH}_{\mathrm{buffer}}\approx 4.7$ and
$\mathrm{pH}_{\mathrm{acid}}>6$, for and acid with p$K_{1}=1$ and p$K_{2}=4$.
The cyan curves of panels (b) and (c) of Figure 8 show that the dibasic salt
produces basic pH stability for acids with higher p$K_{2}$.

##  4 Declarations

###  4.1 Ethical Approval

This work does not involve studies in humans and/or animals. There was no need
for ethics committee approval.

###  4.2 Competing interests

The authors declare no competing interests.

###  4.3 Authors’ contributions

Juan C. Morales made analytical and numerical calculations. Carlos A. Arango
performed analytical and numerical calculations, wrote the manuscript,
prepared the figures, and performed the analysis of results.

###  4.4 Funding

This work has been financed by the OMICAS program, Project ID:
FP44842-217-2018, and the internal research grants of Universidad Icesi. The
OMICAS program acronym stands for “In-silico Multiscale Optimization of
Sustainable Agricultural Crops”, a member of the Scientific Colombia
Ecosystem, sponsored by the World Bank, and the Colombian Ministries of
Science, Technology and Innovation (Minciencias), Education, Industry and
Tourism, and the ICETEX.

###  4.5 Availability of data and materials

The data and Wolfram Mathematica codes used for this study are available from
the corresponding author on request.

##  5 Appendix

###  5.1 Mathematical solution of $P=0$

####  5.1.1 The resolvent cubic equation

The solution of equation $P=0$ by Ferrari’s method requires to find its
resolvent cubic equation [13, 14]. The standard procedure to obtain the
resolvent cubic of a quartic equation begins by writing $P=0$ in its
equivalent form

$\left(x^{2}+\tfrac{1}{2}c_{3}x\right)^{2}=\left(\tfrac{1}{4}c_{3}^{2}-c_{2}\right)x^{2}-c_{1}x-c_{0}.$
(77)

The addition of a quantity $y/2$ inside the squared term of the left hand
side, and the addition of the compensation terms on the right hand side gives,
after simplification,

$\left(x^{2}+\tfrac{1}{2}c_{3}x+\tfrac{y}{2}\right)^{2}=\left(\tfrac{1}{4}c_{3}^{2}-c_{2}+y\right)x^{2}+(\tfrac{1}{2}c_{3}y-c_{1})x-c_{0}+\tfrac{1}{4}y^{2}.$
(78)

The left hand side of this equation can be written as a complete square, that
is, equation (78) can be written

$\left(tx+\frac{c_{3}y-2c_{1}}{4t}\right)^{2}=\left(\tfrac{1}{4}c_{3}^{2}-c_{2}+y\right)x^{2}+(\tfrac{1}{2}c_{3}y-c_{1})x-c_{0}+\tfrac{1}{4}y^{2},$
(79)

with $t=t(y)\neq 0$, given that $t^{2}=\tfrac{1}{4}c_{3}^{2}-c_{2}+y$, and

$\left(\frac{c_{3}y-2c_{1}}{4t}\right)^{2}=\tfrac{1}{4}y^{2}-c_{0}.$ (80)

The expansion of equation (80) gives, after simplification, the resolvent
cubic $R=0$, with

$R=y^{3}-c_{2}y^{2}+\left(c_{1}c_{3}-4c_{0}\right)y+\left(4c_{0}c_{2}-c_{0}c_{3}^{2}-c_{1}^{2}\right).$
(81)

This equation has three roots $y_{i}$, $i=1,2,3$. The use of one of these
roots in equations (78) and (79) gives

$\left(x^{2}+\tfrac{1}{2}c_{3}x+\tfrac{y_{i}}{2}\right)^{2}=\left(t_{i}x+\frac{c_{3}y_{i}-2c_{1}}{4t_{i}}\right)^{2},$
(82)

with $i=1,2,3$ and $t_{i}=t(y_{i})$. Each of these equations split in two
quadratic equations,

$\displaystyle
x^{2}+\left(\tfrac{1}{2}c_{3}-t_{i}\right)x+\tfrac{1}{2}y_{i}-\frac{c_{3}y_{i}-2c_{1}}{4t_{i}}$
$\displaystyle=0,$ (83) $\displaystyle
x^{2}+\left(\tfrac{1}{2}c_{3}+t_{i}\right)x+\tfrac{1}{2}y_{i}+\frac{c_{3}y_{i}-2c_{1}}{4t_{i}}$
$\displaystyle=0,$ (84)

with $i=1,2,3$. The roots of $P=0$ satisfy these quadratic equations [13, 14].
The discriminants of the quartic equation, $\Delta$, and its resolvent cubic
equation, $\Delta_{\mathrm{rc}}$, are identical [13, 14]. The restriction
$k_{1}\geq 4k_{2}$ gives that $\Delta>0$, hence $\Delta_{\mathrm{rc}}>0$ and
the resolvent cubic equation (81) must have three distinct real roots. For the
case $\Delta_{\mathrm{rc}}<0$ the cubic $R=0$ has one real root and two non-
real complex conjugate roots [13, 14].

####  5.1.2 Solution of the resolvent cubic equation $R=0$

The third order polynomial equation $R=0$ can be solved by Cardano’s method.
The change of variable $y=\bar{y}+\frac{c_{2}}{3}$ gives the depressed cubic
equation $R_{\mathrm{dc}}=0$, with

$R_{\mathrm{dc}}=\bar{y}^{3}+\bar{p}\bar{y}+\bar{q},$ (85)

and

$\displaystyle\bar{p}$ $\displaystyle=c_{1}c_{3}-\frac{c_{2}^{2}}{3}-4c_{0},$
(86) $\displaystyle\bar{q}$
$\displaystyle=\frac{8c_{0}c_{2}}{3}+\frac{c_{1}c_{2}c_{3}}{3}-\frac{2c_{2}^{3}}{27}-c_{1}^{2}-c_{0}c_{3}^{2},$
(87)

and discriminant $\Delta_{\mathrm{dc}}=-4\bar{p}^{3}-27\bar{q}^{2}$, which is
equal to $\Delta$ and $\Delta_{\mathrm{rc}}$ [13, 14].

The use of Vieta’s substitution, $\bar{y}=\bar{z}-\frac{\bar{p}}{3\bar{z}}$,
gives the polynomial equation

$\bar{z}^{3}-\frac{\bar{p}^{3}}{27\bar{z}^{3}}+\bar{q}=0.$ (88)

Multiplication of (88) by $\bar{z}^{3}$ gives

$\bar{z}^{6}+\bar{q}\bar{z}^{3}-\frac{\bar{p}^{3}}{27}=0,$ (89)

which is equivalent to the quadratic equation
$\xi^{2}+\bar{q}\xi-\tfrac{\bar{p}^{3}}{27}=0$, in the variable
$\xi=\bar{z}^{3}$, with roots

$\begin{split}\xi_{1,2}&=-\frac{\bar{q}}{2}\pm\sqrt{\frac{27\bar{q}^{2}+4\bar{p}^{3}}{108}}\\\
&=-\frac{\bar{q}}{2}\pm\frac{1}{2}\sqrt{-\frac{\Delta_{\mathrm{dc}}}{27}}.\end{split}$
(90)

The physical case of diprotic acids with $k_{1}\geq 4k_{2}$ [9, 11] gives
$\Delta_{\mathrm{dc}}>0$ for $\bar{c}_{\mathrm{a}}\geq 0$ and
$\bar{c}_{\mathrm{b}}\geq 0$, therefore $\xi_{1,2}$ are a complex conjugate
pair. On the other hand, for the less common case of diprotic acids with
$k_{1}<4k_{2}$ is possible to have $\Delta_{\mathrm{dc}}<0$, hence $\xi_{1,2}$
are real conjugates, on part of the plane
$\bar{c}_{\mathrm{a}}$-$\bar{c}_{\mathrm{b}}$, with $\xi_{1}>\xi_{2}$.

It is convenient to define $\zeta=\xi_{1}$ and $\zeta^{*}=\xi_{2}$ for the
case $k_{1}\geq 4k_{2}$ and $\xi=\xi_{1}$ and $\bar{\xi}=\xi_{2}$ for the case
$k_{1}<4k_{2}$, with $\zeta^{*}$ and $\bar{\xi}$ as the complex and real
conjugate of $\zeta$ and $\xi$ respectively.

The use the polar representation for $\zeta$ gives
${\zeta=\|\zeta\|e^{i\theta}}$ for the case $k_{1}\geq 4k_{2}$, with

$\displaystyle\|\zeta\|$
$\displaystyle=\frac{1}{2}\sqrt{\bar{q}^{2}+\frac{\Delta_{\mathrm{dc}}}{27}}=\sqrt{\frac{-\bar{p}^{3}}{27}},$
(91) $\displaystyle\theta$
$\displaystyle=\arctan{\left(-\frac{\bar{q}}{2},\frac{\sqrt{\Delta_{\mathrm{dc}}}}{6\sqrt{3}}\right)},$
(92)

with $\theta\in(0,\pi)$ as the angle between $\zeta$ and the positive real
axis on the Argand plane. The angle $\theta$ is related to the trigonometric
solution obtained by Nickalls for the roots of the cubic equation [15]. The
polar representation of $\xi$ and $\bar{\xi}$, case $k_{1}<4k_{2}$, gives
$\xi=|\xi_{1}|e^{i\theta_{1}}$ and $\bar{\xi}=|\xi_{2}|e^{i\theta_{2}}$ with
$\theta_{1,2}=\frac{\pi}{2}(1-\operatorname{sgn}{\xi_{1,2}})$, and
$\xi>\bar{\xi}$.

The roots $\bar{y}$ of the depressed cubic equation $R_{\mathrm{dc}}=0$ are
given by Cardano’s formula

$\bar{y}=\alpha+\beta,$ (93)

where $\alpha=\sqrt[3]{\zeta}$ and $\beta=\sqrt[3]{\zeta^{*}}$ for $k_{1}\geq
4k_{2}$, and $\alpha=\sqrt[3]{\xi}$ and $\beta=\sqrt[3]{\bar{\xi}}$ for
$k_{1}<4k_{2}$. The cubic roots $\alpha$ and $\beta$ have three values each,
$\alpha_{n}$ and $\beta_{n}$, with $n=0,1,2$. The combined use of the three
roots $\alpha$ and the three roots $\beta$ must give the three roots of
$R_{\mathrm{dc}}=0$.

The cubic roots $\alpha$ and $\beta$, for the case $k_{1}\geq 4k_{2}$, are
given by

$\displaystyle\alpha_{n}$
$\displaystyle=\sqrt[3]{\|\zeta\|}\exp{\left(i\left(\frac{\theta}{3}+\frac{2n\pi}{3}\right)\right)},$
(94) $\displaystyle\beta_{n}$
$\displaystyle=\sqrt[3]{\|\zeta\|}\exp{\left(i\left(-\frac{\theta}{3}+\frac{2n\pi}{3}\right)\right)},$
(95)

with $n=0,1,2$, and

$\sqrt[3]{\|\zeta\|}=\tfrac{1}{3}\sqrt{1+k_{1}Q_{1}+k_{1}^{2}Q_{2}}.$ (96)

for which $Q_{1}$ and $Q_{2}$ are

$\displaystyle Q_{1}$
$\displaystyle=-3k_{2}\bar{c}_{\mathrm{b}}^{2}+(6\bar{c}_{\mathrm{a}}k_{2}+1)\bar{c}_{\mathrm{b}}+2\bar{c}_{\mathrm{a}}-14k_{2},$
(97) $\displaystyle Q_{2}$
$\displaystyle=k_{2}^{2}+(4\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}})k_{2}+3+(\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}})^{2}.$
(98)

The case $k_{1}<4k_{2}$ has

$\displaystyle\alpha_{n}$
$\displaystyle=\sqrt[3]{|\xi_{1}|}\exp{\left(i\left(\frac{\theta_{1}}{3}+\frac{2n\pi}{3}\right)\right)},$
(99) $\displaystyle\beta_{n}$
$\displaystyle=\sqrt[3]{|\xi_{2}|}\exp\left(i\left(\frac{\theta_{2}}{3}+\frac{2n\pi}{3}\right)\right),$
(100)

with $n=0,1,2$.

Since for the case $k_{1}\geq 4k_{2}$ the cubic equation $R_{\mathrm{dc}}=0$
has three real roots, the addition of two cubic roots $\alpha_{n}+\beta_{m}$
must give a real number. This is possible only if
$\operatorname{Im}(\alpha_{n})=-\operatorname{Im}(\beta_{m})$. There are only
three possible combinations that fulfill this requirement:

$\displaystyle\alpha_{0}+\beta_{0}$
$\displaystyle=2\sqrt[3]{\|\zeta\|}\cos{\left(\tfrac{\theta}{3}\right)},$
(101) $\displaystyle\alpha_{2}+\beta_{1}$
$\displaystyle=-2\sqrt[3]{\|\zeta\|}\cos{\left(\tfrac{\theta+\pi}{3}\right)},$
(102) $\displaystyle\alpha_{1}+\beta_{2}$
$\displaystyle=-2\sqrt[3]{\|\zeta\|}\sin{\left(\tfrac{\theta+2\pi}{3}\right)},$
(103)

which are the roots of $R_{\mathrm{dc}}=0$: $\bar{y}_{1}$, $\bar{y}_{2}$, and
$\bar{y}_{3}$, respectively, with $\bar{y}_{1}>\bar{y}_{2}>\bar{y}_{3}$.

The case $k_{1}<4k_{2}$ has only one real solution, and a complex conjugate
pair. Since $\xi>\bar{\xi}$, there are three possibilities:

  * •

$\theta_{1}=\theta_{2}=0$: the roots are $\alpha_{i}+\beta_{i}$ with
$i=0,1,2$. The root $\bar{y}_{1}=\alpha_{0}+\beta_{0}$ is the only real
solution, $\bar{y}_{1}=\sqrt[3]{|\xi_{1}|}+\sqrt[3]{|\xi_{2}|}$.

  * •

$\theta_{1}=\theta_{2}=\pi$: the roots are $\alpha_{i}+\beta_{i}$ with
$i=0,1,2$. The root $\bar{y}_{1}=\alpha_{1}+\beta_{1}$ is the only real
solution, $\bar{y}_{1}=-(\sqrt[3]{|\xi_{1}|}+\sqrt[3]{|\xi_{2}|})$.

  * •

$\theta_{1}=0$, $\theta_{2}=\pi$: the roots are $\alpha_{0}+\beta_{1}$,
$\alpha_{1}+\beta_{0}$, and $\alpha_{2}+\beta_{2}$. The root
$\bar{y}_{1}=\alpha_{0}+\beta_{1}$ is the only real solution,
$\bar{y}_{1}=\sqrt[3]{|\xi_{1}|}-\sqrt[3]{|\xi_{2}|}$.

In summary, the solution $\bar{y}_{1}$ of the depressed cubic equation is

$\bar{y}_{1}=\begin{cases}\tfrac{2}{3}\sqrt{1+k_{1}Q_{1}+k_{1}^{2}Q_{2}}\cos{\left(\tfrac{\theta}{3}\right)},&\Delta_{\mathrm{dc}}>0,\\\
\sqrt[3]{|\xi_{1}|}+\sqrt[3]{|\xi_{2}|},&\Delta_{\mathrm{dc}}<0,\;\xi_{1}>0,\;\xi_{2}>0,\\\
-(\sqrt[3]{|\xi_{1}|}+\sqrt[3]{|\xi_{2}|}),&\Delta_{\mathrm{dc}}<0,\;\xi_{1}<0,\;\xi_{2}<0,\\\
\sqrt[3]{|\xi_{1}|}-\sqrt[3]{|\xi_{2}|},&\Delta_{\mathrm{dc}}<0,\;\xi_{1}>0,\;\xi_{2}<0.\end{cases}$
(104)

The root $y_{1}$ of the resultant cubic equation $R=0$ is given by

$y_{1}=\bar{y}_{1}-\tfrac{1+k_{1}\left(\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}}-k_{2}\right)}{3}$
(105)

This root is substituted in the quadratic equations (83) and (84), with
$t_{1}=t(y_{1})$ given by

$t_{1}=\sqrt{1+\tfrac{1}{4}(\bar{c}_{\mathrm{b}}+k_{1})^{2}+k_{1}\left(\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}}-k_{2}\right)+y_{1}}.$
(106)

The four roots of the quartic equation $P=0$ are given by

$\displaystyle x_{1,2}$
$\displaystyle=\tfrac{1}{2}\left(-\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}-t_{1}\right)\pm\sqrt{\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}-t_{1}\right)^{2}-2y_{1}+\tfrac{(\bar{c}_{\mathrm{b}}+k_{1})y_{1}+2k_{1}\left(1+(2\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}})k_{2}\right)}{t_{1}}}\right),$
(107) $\displaystyle x_{3,4}$
$\displaystyle=\tfrac{1}{2}\left(-\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}+t_{1}\right)\pm\sqrt{\left(\tfrac{\bar{c}_{\mathrm{b}}+k_{1}}{2}+t_{1}\right)^{2}-2y_{1}-\tfrac{(\bar{c}_{\mathrm{b}}+k_{1})y_{1}+2k_{1}\left(1+(2\bar{c}_{\mathrm{a}}-\bar{c}_{\mathrm{b}})k_{2}\right)}{t_{1}}}\right).$
(108)

Only the roots $x_{1}$ and $x_{3}$ can have physical significance, and $x$ is
given by

$x=\begin{cases}x_{1},&\Delta_{\mathrm{dc}}>0,\\\
x_{1},&\Delta_{\mathrm{dc}}<0,\;\xi_{1}>0,\;\xi_{2}>0,\\\
x_{3},&\Delta_{\mathrm{dc}}<0,\;\xi_{1}<0,\;\xi_{2}<0,\\\
x_{3},&\Delta_{\mathrm{dc}}<0,\;\xi_{1}>0,\;\xi_{2}<0.\end{cases}$ (109)

## References

  * Scholz and Kahlert [2018] Scholz, F.; Kahlert, H. Acid–base equilibria of amino acids: microscopic and macroscopic acidity constants. _The Textbook Journal of Chemistry_ 2018, _4_ , 1–9
  * Hamm et al. [2015] Hamm, L.; Nakhoul, N.; Hering-Smith, K. Acid-Base Homeostasis. _Clin J Am Soc Nephrol_ 2015, _10_ , 2232–42
  * Doney et al. [2020] Doney, S. C.; Busch, D. S.; Cooley, S. R.; Kroeker, K. J. The Impacts of Ocean Acidification on Marine Ecosystems and Reliant Human Communities. _Annual Review of Environment and Resources_ 2020, _45_ , 83–112
  * Denbigh [1981] Denbigh, K. G. _The Principles of Chemical Equilibrium: With Applications in Chemistry and Chemical Engineering_ , 4th ed.; Cambridge University Press, 1981
  * Burgot [2012] Burgot, J.-L. _Ionic Equilibria in Analytical Chemistry_ ; Springer, 2012
  * Skoog et al. [2022] Skoog, D. A.; West, D. M.; Holler, F. J.; Crouch, S. R. _Fundamentals of Analytical Chemistry_ , 10th ed.; Cengage, 2022
  * Kalka [2021] Kalka, H. Polyprotic Acids and Beyond—An Algebraic Approach. _Chemistry_ 2021, _3_ , 454–508
  * Caicedo et al. [2023] Caicedo, A.; Morales, J. C.; Arango, C. A. Closed-form expressions for monoprotic weak acid aqueous solutions. _ChemTexts_ 2023, _9_
  * Adams [1916] Adams, E. Q. Relations Between the Constants of Dibasic Acids and of Amphoteric Electorolytes. _Journal of the American Chemical Society_ 1916, _38_ , 1503–1510
  * CRC Handbook [2007] CRC Handbook, _CRC Handbook of Chemistry and Physics_ , 88th ed.; CRC Press, 2007
  * Mchedlov-Petrossyan [2019] Mchedlov-Petrossyan, N. POLYPROTIC ACIDS IN SOLUTION: IS THE INVERSION OF THE CONSTANTS OF STEPWISE DISSOCIATION POSSIBLE? _Ukrainian Chemistry Journal_ 2019, _85_ , 3–45
  * Simms [1926] Simms, H. S. Dissociation of Polyvalent Substances II. Relation of Constants to Chemical Structure. _Journal of the American Chemical Society_ 1926, _48_ , 1251–1261
  * Dickson [1914] Dickson, L. E. _Elementary Theory of Equations_ , 1st ed.; John Wiley and Sons, Inc., 1914
  * Dickson [1922] Dickson, L. E. _First Course in the Theory of Equations_ , 1st ed.; John Wiley and Sons, Inc., 1922
  * Nickalls [1993] Nickalls, R. W. D. A new approach to solving the cubic: Cardan’s solution revealed. _The Mathematical Gazette_ 1993, _77_ , 354–359