# acid dissociation constants for surface sites and for water
Ka1 = 1.00e-6
Ka2 = 2.00e-8
Kw = 1.00e-14

# complexation constants for lead
Kalopb = 1.23e6

# mass balances
Caloh = 1e-4
Cpb = 1e-7

# set up master variable
ph = seq(0, 14, 0.01)

# calculate concentrations of hydronium and hydroxide ions
h = 10^-ph
oh = Kw/h

# calculate concentrations of surface sites
aloh2 = Caloh/(1 + Ka1/h + Ka1*Ka2/h^2)
aloh = Ka1*aloh2/h
alo = Ka1*Ka2*aloh2/h^2

# calculate concentrations of lead
pb = Cpb/(1 + Kalopb*alo)
alopb = Kalopb*alo*pb

# plot the results for lead
plot(ph, 100*alopb/Cpb, type = "l", lty = 1, lwd = 2, xlim = c(0, 14), 
     ylim = c(0, 100), ylab = "%Pb or %Surface Sites")

# overlay the results for the surface sites
lines(ph, 100*aloh2/Caloh, col = "lightgray", lty = 2, lwd = 2)
lines(ph, 100*aloh/Caloh, col = "lightgray", lty = 3, lwd = 2)
lines(ph, 100*alo/Caloh, col = "lightgray", lty = 4, lwd = 2)

# add a legend
legend(x = "left", legend = c("AlOH2+", "AlOH", "AlO-", "AlOPb+"), 
       lwd = c(2, 2, 2, 2), 
       col = c("lightgray", "lightgray", "lightgray", "black"), 
       lty = c(2, 3, 4, 1), bty = "n")