How to find the concentrations of the acid-base equilibrium between sulfuric acid and sulfur dioxide for a known pH using Octave?

Question

The system is an aqueous solution of sulfuric acid and sulfur dioxide. I know how much sulfur there is and I know the pH. The system is in equilibrium. I would like to use Octave to solve the system.

I set up the following equations, ignoring $\ce{SO3^2-}$ because the pH is around 1.5.

By definition: $$ \mathrm{pH} = -\log [\ce{H+}]$$

Conservation of mass: $$[\ce{H2SO4}] + [\ce{H2SO3}] + [\ce{HSO4-}] + [\ce{SO4^2-}] + [\ce{HSO3-}] = M $$

Henderson-Hasselbalch equations for the system: \begin{align} \mathrm{pH} &= \mathrm{p}K_\mathrm{a}(\ce{HSO3-}) + \log \frac{[\ce{HSO3-}]}{[\ce{H2SO3}]}\\ \mathrm{pH} &= \mathrm{p}K_\mathrm{a}(\ce{SO4^2-}) + \log \frac{[\ce{SO4^2-}]}{[\ce{HSO4-}]}\\ \mathrm{pH} &= \mathrm{p}K_\mathrm{a}(\ce{HSO4-}) + \log \frac{[\ce{HSO4-}]}{[\ce{H2SO4}]}\\ \end{align}

I then defined the solution vector to have the following order: $$ x = [ [\ce{H+}], [\ce{HSO4-}], [\ce{SO4^2-}], [\ce{HSO3-}], [\ce{H2SO4}], [\ce{H2SO3}] ]$$

and my code for Octave is attached. For some reason I always get negative values for some of the results.

function y = phCalc(x)
   pH = 1.87;
   pKa1 = -5;   % this is a guess to drive H2SO4 to ionization in this system
   pKa2 = 1.987;% from a textbook
   pKa3 = 1.857;% from a text book
   M = 0.0012 / 32 * 1000 ; % total sulfur moles per volume
   y(1) = 10^(-pH) - x(1);
   y(2) = x(5) * 10^(pH-pKa1) - x(2);
   y(3) = x(2) * 10^(pH-pKa2) - x(3);
   y(4) = x(1) - x(2) - x(3) - x(4);
   y(5) = x(2) + x(3) + x(4) + x(5) + x(6) - M;
   y(6) = x(4) / 10^(pH-pKa3) - x(6);
endfunction

Answer 1

Your code seems to be calculating the differences between a guessed concentration vector x and a calculated concentration vector. These differences in concentrations between guessed and calculated are y. For example, the guessed pH is apparently x(1), and the calculated pH is 10^(-pH). The difference between them is y(1).

Thus, its no surprise that the y values are sometimes negative. To solve the problem you need to feed your function to a numerical solver such as fsolve. That function will take your guessed solution x, other parameters such as pH, the various pKas, and M, and your function phCalc as input, and return the value for x which makes all the y variables zero. That solution is the one you want.

On the chemistry side, you mention sulfur dioxide, but I don't see any equations for the solubility of sulfur dioxide gas in water. Additionally, sulfur dioxide is reduced relative to sulfuric acid, so there would need to be an additional reaction, and parameter (extent of reaction or degree of reduction or something) that describes the relative amounts of reduced sulfur (as sulfur dioxide, bisulfite, or sulfite) vs. oxidized sulfur (as sulfate, bisulfate, sulfuric acid, or sulfur trioxide).

There is also the issue of charge balance: what concentration of inert cations (e.g. sodium or potassium or whatever) are you assuming? Charge balance must hold in bulk aqueous solution. If you don't have any sodium / potassium etc., then you should include the self ionization of water as an equation, hydroxide as a species, and y(7) = x(1) - x(2) - 2*x(3) - x(4) - x(OH) as an equation. If you do have sodium/potassium/whatever at significant concentration, you can neglect hydroxide and just go with the equation y(7) = x(1) - x(2) - 2*x(3) - x(4) + x(Na) where x(Na) is the cation concentration, which I have assumed is sodium.

#

Answer 2

What you want to solve, is the following system: $$ \frac{k_W}{\ce{[H+]}} - \ce{[H+]} + \frac{\ce{[H2SO4]0}\left(\ce{[H+]}~k_{a1}+2~k_{a1}k_{a2}\right)} {\ce{[H+]}^2+\ce{[H+]}~k_{a1}+k_{a1}k_{a2}} + \frac{\ce{[H2SO3]0}\left(\ce{[H+]}~k_{a3}+2~k_{a3}k_{a4}\right)} {\ce{[H+]}^2+\ce{[H+]}~k_{a3}+k_{a3}k_{a4}} = 0 $$ for both concentrations of $\ce{[H2SO4]0}$ and $\ce{[H2SO3]0}$ with the constraint that both concentrations add up to $\mathrm{0.0375~mol/L}$.

I have no idea how to solved this with Octave (see below), but as I am more used to Mathematica, that was my first approach to see where the error came from. The problem that you have with those negative solutions is, that the constrain about the concentrations does not allow two positive concentrations. Have a look at the following picture:

$\hskip3.5cm$

The pH of a pure 0.0375 M sulfuric acid solution would be somewhat around 1.35. Changing the mixture by reducing the amount of sulfuric acid and increasing the amount of sulfurous acid leads to a rising pH. A pure 0.0375 M sulfurous acid solution would then have a pH of about 1.77.

But for concentrations below around $\mathrm{0.025~mol/L}$ the system is solvable, what can be seen in the next image. You could obtain the pH of 1.87 with a mixture of 0.00086859742 M $\ce{H2SO4}$ and 0.024131403 M $\ce{H2SO3}$.

$\hskip3.5cm$

---

Now doing this with Octave

I never used Octave before, so there might be a better solution. But this is my code:

>> function y = Conc(x)
ka1 = 10^3; % wikipedia, sulfuric_acid
ka2 = 10^(-1.99); % wikipedia, sulfuric_acid
ka3 = 10^(-1.857); % wikipedia, sulfurous_acid
ka4 = 10^(-7.172); % wikipedia, sulfurous_acid
ctot = 0.025; % total sulfur concentration in mol/L
ph = 1.87; h = 10^-(ph);
oh = 10^(-14)/h;
xH2SO4 = (h * ka1 + 2 * ka1 * ka2)/(h^2 + h * ka1 + ka1 * ka2);
xH2SO3 = (h * ka3 + 2 * ka3 * ka4)/(h^2 + h * ka3 + ka3 * ka4);
y(1) = x(1) + x(2) - ctot;
y(2) = oh - h + x(1) * xH2SO4 + x(2) * xH2SO3;
endfunction
>> [x, fval, info] = fsolve(@Conc, [1; 2])
x =
  8.6860e-004
  2.4131e-002
fval =
  6.3704e-011  2.8976e-011
info =  1

And as you can see, it gives exactle the same solutions as Mathematica. (fval is the value for y(1) and y(2) that are near zero, so everything is fine. info = 1 means, that everything the solution has converged.)