# Calculating the pH in weak acids#

← Previous | Oldest ↻
August 02, 2022 () by Anders Lervik | Category: chemistry
Tagged: chemistry

## Introduction#

If you have taken a chemistry course, you have probably calculated the pH in a weak acid. A typical problem is: What is the pH in a solution of the weak acid HA (with concentration C) in water? The solution can then be found by solving,

$K_\text{a} = \frac{x^2}{C - x}$

where $$K_\text{a}$$ is the acid dissociation constant, and $$x$$ is the (unknown) concentration of hydrogen ions. We assume here that we do not have to take the autoionization of water into consideration, but when is this true? I will check this with some calculations here.

Note

To give you the answer right now: Solve the equation given above. If the pH is 6 or smaller, this is good enough. If the pH is above 6, solve the more complex equations given below.

## Calculating the pH in a weak acid#

I consider the dissociation of the weak acid HA:

$[\text{HA}] \rightleftharpoons [\text{H}^+] + [\text{A}^-].$

We have the following equations:

• The equilibrium constant for the acid dissociation: $$K_\text{a} = \frac{[\text{H}^+] [\text{A}^-]}{[\text{HA}]}$$.

• The equilibrium constant for autoionization of water: $$K_\text{w} = [\text{H}^+] [\text{OH}^-]$$.

• The mass balance (based on element A): $$[\text{HA}]_0 = [\text{HA}] + [\text{A}^-]$$, where $$[\text{HA}]_0$$ is the initial concentration of HA.

• The electroneutrality: $$[\text{H}^+] = [\text{A}^-] + [\text{OH}^-]$$.

We can rewrite the electroneutrality as $$[\text{A}^-] = [\text{H}^+] - \frac{K_\text{w}}{[\text{H}^+]}$$, combine it with the mass balance and the equilibrium constant to get:

$K_\text{a} = \frac{[\text{H}^+] \left( [\text{H}^+] - \frac{K_\text{w}}{[\text{H}^+]}\right)}{[\text{HA}]_0 - \left([\text{H}^+] - \frac{K_\text{w}}{[\text{H}^+]} \right)},$

which is an equation with one unknown ($$[\text{H}^+]$$) we can solve.

We can also recover the usual General Chemistry 101 formula by assuming that $$[\text{H}^+]^2 \gg K_\text{w}$$:

$K_\text{a} \approx \frac{[\text{H}^+]^2}{[\text{HA}]_0 - [\text{H}^+]}.$

### Solving the approximate equation#

The approximate equation is a quadratic equation, $$[\text{H}^+]^2 + [\text{H}^+] K_\text{a} - K_\text{a} [\text{HA}]_0 = 0.$$ and we can solve it by:

[1]:

import numpy as np
import warnings

[2]:

@np.vectorize
def calculate_ph_approximate(c_ha_zero, ka):
"""Solve the approximate equation"""
poly = [1, ka, -ka * c_ha_zero]
discriminant = poly[1] ** 2 - 4 * poly[0] * poly[2]
if discriminant < 0:
raise ValueError("No positive roots!")
# Calculate the positive root:
x = (-poly[1] + np.sqrt(discriminant)) / (2 * poly[0])
return -np.log10(x)


Let us try this for Example 15.9.2 from Principles of modern chemistry (where the pH is calculated to be 6.65):

[3]:

calculate_ph_approximate(0.001, 4.0e-11)

[3]:

array(6.69901343)


The approximate solution is not spot on here, but it is pretty close.

### Solving the full equation#

The full equation can be rewritten as a cubic equation:

$[\text{H}^+]^3 + [\text{H}^+]^2 K_\text{a} - [\text{H}^+] \left(K_\text{w} + K_\text{a} [\text{HA}]_0\right) - K_\text{w}K_\text{a} = 0.$

According to Descartes’ rule this equation will have exactly one positive root. We can also use the base constant:

$[\text{OH}^-]^3 + [\text{OH}^-]^2 \left(K_\text{b} + [\text{HA}]_0 \right) - [\text{OH}^-] K_\text{w} - K_\text{w}K_\text{b} = 0.$

I will use the first equation and solve for $$[\text{H}^+]$$, but if there are some issues with the solution, I will use the equation for $$[\text{OH}^-]$$. To check if there are some issues, I define two simple checks below:

[4]:

def check_solution_acid(c_ha_zero, ka, kw, x, print_warning=False):
"""Do some checks for the solution."""
c_OH = kw / x
c_A = x - c_OH
w = 0
if c_A < 0:
if print_warning:
warnings.warn(f"Unphysical [A-] = {c_A} for {c_ha_zero}, {ka}")
w += 1
c_HA = c_ha_zero - c_A
if c_HA < 0:
if print_warning:
warnings.warn(f"Unphysical [HA] = {c_HA} for {c_ha_zero}, {ka}")
w += 1
ka_ = x * c_A / c_HA
if not np.isclose(ka_, ka, atol=1e-6):
if print_warning:
warnings.warn(
f"Solution is inconsistent for {c_ha_zero}, {ka_} != {ka}"
)
w += 1
return w

def check_solution_base(c_ha_zero, ka, kw, x, print_warning=False):
"""Do some checks for the solution."""
c_H = kw / x
c_A = c_H - x
w = 0
if c_A < 0:
if print_warning:
warnings.warn(f"Unphysical [A-] = {c_A} for {c_ha_zero}, {ka}")
w += 1
c_HA = c_ha_zero - c_A
if c_HA < 0:
if print_warning:
warnings.warn(f"Unphysical [HA] = {c_HA} for {c_ha_zero}, {ka}")
w += 1
kb = kw / ka
kb_ = x * c_HA / c_A
if not np.isclose(kb_, kb, atol=1e-6):
if print_warning:
warnings.warn(
f"Solution is inconsistent for {c_ha_zero}, {kb_} != {kb}"
)
w += 1
return w

[5]:

@np.vectorize
def calculate_ph_full(c_ha_zero, ka, kw=1e-14):
"""Solve the full equation"""
poly_h = [1, ka, -(kw + ka * c_ha_zero), -ka * kw]
roots_h = np.roots(poly_h)
# Select the positive solution:
x_h = roots_h[roots_h > 0][0]
# Do some checks:
w = check_solution_acid(c_ha_zero, ka, kw, x_h, print_warning=False)
if w > 0:  # Fix problems by solving the base equation:
kb = kw / ka
poly_oh = [1, (c_ha_zero + kb), -kw, -kb * kw]
roots_oh = np.roots(poly_oh)
x_oh = roots_oh[roots_oh > 0][0]
x_h = kw / x_oh
check_solution_base(c_ha_zero, ka, kw, x_oh, print_warning=True)
return -np.log10(x_h)


I will also try this method for Example 15.9.2 from Principles of modern chemistry (where the pH is calculated to be 6.65):

[6]:

calculate_ph_full(0.001, 4.0e-11)

[6]:

array(6.65054607)


Ah, this is a more accurate result! I will run some more comparisons below.

## Comparing the solutions#

[7]:

from matplotlib import pyplot as plt
import seaborn as sns
import matplotlib as mpl

%matplotlib inline

mpl.rcParams["ytick.minor.visible"] = True
mpl.rcParams["xtick.minor.visible"] = True
sns.set_theme(style="ticks", context="talk", palette="muted")

[8]:

ka = [1e-2, 1e-4, 1e-8, 1e-10]
concentrations = 10.0 ** np.arange(-12, 3, 1)

fig, axes = plt.subplots(
constrained_layout=True, ncols=4, figsize=(12, 3), sharex=True, sharey=True
)

for kai, axi in zip(ka, axes):
ph_approx = calculate_ph_approximate(concentrations, kai)
ph_full = calculate_ph_full(concentrations, kai)
axi.scatter(-np.log10(concentrations), ph_approx, label="Approximate")
axi.scatter(-np.log10(concentrations), ph_full, label="Full solution")
axi.set(xlabel="pC", title=f"pKa = {-np.log10(kai)}")
axes[0].set(ylabel="pH")
axes[0].legend()
sns.despine(fig=fig)


We do get some errors for low concentrations here. Let us investigate a grid of points:

[9]:

ka = 10.0 ** np.arange(-14, 1, 0.2)
concentrations = 10.0 ** np.arange(-14, 2, 0.2)
ka_grid, conc_grid = np.meshgrid(ka, concentrations)
ph_approx = calculate_ph_approximate(conc_grid, ka_grid)
ph_full = calculate_ph_full(conc_grid, ka_grid)
error_grid = abs(ph_approx - ph_full) / abs(ph_full)

[10]:

fig, ax = plt.subplots(constrained_layout=True, figsize=(10, 8))
im = ax.contourf(
-np.log10(conc_grid),
-np.log10(ka_grid),
error_grid * 100,
cmap="vlag",
levels=40,
)

im2 = ax.contour(
-np.log10(conc_grid),
-np.log10(ka_grid),
error_grid * 100,
levels=[1, 10],
colors="w",
)

ax.clabel(
im2, im2.levels, inline=True, fontsize="small", fmt="Error = %.1f %%"
)

ax.set(xlabel="pC", ylabel="pKa")

ax.scatter(-np.log10(0.001), -np.log10(4e-11), color="white")
ax.annotate(
"Example calculation",
(-np.log10(0.001), -np.log10(4e-11)),
xytext=(-100, -75),
color="white",
textcoords="offset points",
arrowprops={"arrowstyle": "->"},
)
fig.colorbar(im, ax=ax, label="Error (%)")

[10]:

<matplotlib.colorbar.Colorbar at 0x7fb8165fef40>


To get an error below 1%:

• If the concentration is lower than $$10^{-6.7} = 2 \times 10^{-7}$$ we should not use the approximation.
• If the concentration is higher than $$2 \times 10^{-7}$$ , we can use it as long as $$K_\text{a} [\text{HA}]_0 > 3.3 K_\text{w}$$ (that is $$\text{p}K_\text{a} + \text{pC} < 13.4$$).

For the example calculation I did, the error ($$\text{pH} = 6.7$$ vs. $$\text{pH} = 6.65$$) was around 1% (as can also be seen in the figure above). For this case: $$K_\text{a} [\text{HA}]_0 = 4 \times 10^{-14} > 3.3 K_\text{w}$$. Also, from the figure above, the error should be about 10% for a concentration of $$10^{-5}$$:

[11]:

pH1 = calculate_ph_full(1e-5, 4.0e-11, 1e-14)
pH2 = calculate_ph_approximate(1e-5, 4e-11)
error = abs(pH2 - pH1) / pH1
print(f"Error = {error*100:.2f}%")

Error = 10.13%


Let us also check how the error varies with the final pH:

[12]:

fig, ax = plt.subplots(constrained_layout=True, figsize=(8, 6))
ax.scatter(ph_full.ravel(), 100.0 * error_grid.ravel(), s=10, zorder=2)
ax.set(xlabel="pH", ylabel="Error (%)")
ax.axhline(y=1, ls=":", color="k", lw=3, label="Error = 1%")
idx = np.where(error_grid > 0.01)
pH_min = np.min(ph_full[idx])
limits = ax.get_xlim()
ax.axvspan(
xmin=pH_min,
xmax=max(limits),
lw=3,
label=f"pH > {pH_min:.2f}\n(Error > 1%)",
alpha=0.3,
zorder=1,
)

ax.set_xlim(limits)
ax.legend()

sns.despine(fig=fig)


So the text from Wikipedia:

Note that in this example, we are assuming that the acid is not very weak, and that the concentration is not very dilute, so that the concentration of [OH−] ions can be neglected. This is equivalent to the assumption that the final pH will be below about 6 or so.

is a nice summary. If we use the quadratic equation and get a pH below 6, then we do not need to put in the extra work with solving the cubic equation. Here is a visualization of that rule:

[13]:

fig, ax = plt.subplots(constrained_layout=True, figsize=(10, 8))
im = ax.contourf(
-np.log10(conc_grid),
-np.log10(ka_grid),
error_grid * 100,
cmap="vlag",
levels=40,
)

im2 = ax.contour(
-np.log10(conc_grid),
-np.log10(ka_grid),
error_grid * 100,
levels=[1],
colors="w",
)

ax.clabel(
im2, im2.levels, inline=True, fontsize="small", fmt="Error = %.1f %%"
)

im3 = ax.contour(
-np.log10(conc_grid),
-np.log10(ka_grid),
ph_full,
levels=[1, 2, 3, 4, 5, 6],
colors="w",
linestyles="dashed",
)

ax.clabel(im3, im3.levels, inline=True, fontsize="small", fmt="pH = %.1f")

ax.set(xlabel="pC", ylabel="pKa")

fig.colorbar(im, ax=ax, label="Error (%)")

[13]:

<matplotlib.colorbar.Colorbar at 0x7fb81627ac70>


## Discovering the “rules” by a decision tree#

I found the rules stated above with the following decision tree:

[14]:

import pandas as pd
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import export_graphviz
import graphviz

[15]:

data = {
"ka": -np.log10(ka_grid.ravel()),
"conc": -np.log10(conc_grid.ravel()),
"ka*cons/kw": ka_grid.ravel() * conc_grid.ravel() / 1e-14,
"class": np.zeros_like(ka_grid.ravel()),
}
data["class"][error_grid.ravel() > 0.01] = 1
data = pd.DataFrame(data)
data

[15]:

ka conc ka*cons/kw class
0 1.400000e+01 14.0 1.000000e-14 1.0
1 1.380000e+01 14.0 1.584893e-14 1.0
2 1.360000e+01 14.0 2.511886e-14 1.0
3 1.340000e+01 14.0 3.981072e-14 1.0
4 1.320000e+01 14.0 6.309573e-14 1.0
... ... ... ... ...
5995 4.975930e-14 -1.8 6.309573e+15 0.0
5996 -2.000000e-01 -1.8 1.000000e+16 0.0
5997 -4.000000e-01 -1.8 1.584893e+16 0.0
5998 -6.000000e-01 -1.8 2.511886e+16 0.0
5999 -8.000000e-01 -1.8 3.981072e+16 0.0

6000 rows × 4 columns

[16]:

X = data[["ka", "conc", "ka*cons/kw"]].to_numpy()
y = data["class"].to_numpy()
tree = DecisionTreeClassifier(max_depth=2, random_state=2)
tree.fit(X, y)

dot_data = export_graphviz(
tree,
out_file=None,
feature_names=["pKa", "pC", "Ka*cons/Kw"],
class_names=["Approximation OK", "Approximation not ok"],
rounded=True,
filled=True,
)
from IPython.display import display

display(graphviz.Source(dot_data))


## Summary#

The simplest thing to do is to use the quadratic equation. If the pH is 6 or smaller, this is good enough. If the pH is above 6, redo the calculation with the cubic equation.