Maple and Acid Rain: A look at solution chemistry

William D. Scott

Division of Environmental Sciences
Murdoch University
Perth, Western Australia 6150

The power of Maple is apparent when tackling the solution of multivariate polynomials. Here the equilibrium chemistry of raindrops is considered, in examples with three levels of complexity. Considered first is the simplest case of ionic equilibrium with a single solute species derived from Carbon Dioxide ($\ce{CO2}$). Then, the case of species generated from three interacting sources, including Sulfur Dioxide ($\ce{SO2}$) and Ammonia ($\ce{NH3}$). Lastly, chemical kinetics are included with oxidation to sulfate ion. The results are all applicable to the problem of acid rain, resulting from excess sulfurous emissions and carbon dioxide. With effort from the user, many more interactive species may be included, even organics and solids. Maple can solve these problems, provided the user evolves an appropriate structure.

These three specific cases are important not only in cloud and rain physics and air pollution but also in all our waterways, in our soils and in our foods. The technique is presented as a series of examples. Each uses three types of conditions; ‘equilibrium mass action constants’, ‘mass balances’ and a single electroneutrality equation.

Carbon Dioxide

pH of a solution of Carbon Dioxide ($\ce{CO2}$) in water ($\ce{H2O}$), as may be appropriate near the turn of the century.

Simply put, $\ce{CO2}$ dissolves in water to produce ‘carbonic acid’. This breaks up into hydrogen ions ($\ce{H+}$) and bicarbonate ions ($\ce{HCO3-}$) which, in turn, break up into more hydrogen ions and carbonate ions ($\ce{CO3^{=}}$). Water is an interloper that supplies a small number of ions and effectively links the acid $\ce{H+}$ and the basic $\ce{OH-}$ ion. Here the equilibrium expressions or ‘mass action ratios’, the K values, are expressions of thermodynamic equilibrium. In a steady, equilibrium state

\begin{array}{ll} \ce{CO2(gas) + H2O(lqd) <=> H2CO3} & K_{hc} = \dfrac{\left[ \ce{H2CO3} \right]}{p_{\ce{CO2}}} = 34/1000 \\ \ce{H+ + HCO3- <=> H2CO3} & K_{1c} = \dfrac{\left[ \ce{H2CO3} \right]}{\left[ \ce{H+} \right] \cdot \left[ \ce{HCO3-} \right]} = 2250000 \\ \ce{H+ + CO3^{=} <=> HCO3-} & K_{1c} = \dfrac{\left[ \ce{HCO3-} \right]}{\left[ \ce{H+} \right] \cdot \left[ \ce{CO3^{=}} \right]} = 21400000000 \\ \ce{H+ + OH- <=> H2O} & K_{w} = \dfrac {1}{\left[ \ce{H+} \right] \cdot \left[ \ce{OH-} \right]} = 10^{14} \end{array}

they must apply, together with the electroneutrality equation

\begin{equation*} \ce{[ H+ ] ~=~ [ OH- ] + [ HCO3- ] + 2 * [ CO3^{=} ]} \end{equation*}

The [ ]’s represent ‘concentration of’ in moles per litre. The ratios given are ‘mass action’ ratios, not activities; it is presumed that the concentrations are small and the activity coefficients can be assumed to be unity. A feature here is that all the constants (the K values given for 25°C), are formation constants so that the numbers used are not small fractions. This makes it easier to put into Maple though it is not essential. As Maple is a symbolic manipulator, it might be better to use the symbolic values of the constants. Later, using another language like Fortran, the values are inserted. This does not change the essential structure of the solution but generally will fill pages with symbolic code; it is not particularly enlightening during testing of an idea. Here we only use floating point numbers when they can not be avoided; this makes the work easier for Maple and removes round off error.

The question asked is: If the concentration of $\ce{CO2}$ in the atmosphere is 365 parts per million in the year 2000, what is the expected pH of rainwater? This corresponds to a partial pressure of $\ce{CO2}$ ($p_{\ce{CO2}}$) of $\dfrac{365 \cdot 29/44}{1000000}$ moles of $\ce{CO2}$ per mole of air or about 240/1000000 atmospheres of $\ce{CO2}$. All we do now is combine the above equations, using Maple. It is noteworthy that there are 5 equations and 5 unknowns in this equation set.

The maple codes omit the square brackets and the superscripts:

>pCO2 := 240/10^6:

>H2CO3 := 34/1000*pCO2:

>HCO3 := H2CO3/(2250000*H):

>CO3 := HCO3/(214*10^8*H):

>OH := 1/(10^14*H):

The equations are reformed to help the structure of the solution. In this we simply use the fact that, with each assignment, the variable on the left no longer needs to be considered as a variable. the same could be accomplished by putting the equations as written above into a set and methodically using the subs command. The results should be the same but the above method is best for the Maple novice; this helps with learning as it keeps the expressions simple. Of course, it would also be better to leave the K values in the solution as parameters (though known) so that the values can be easily altered. This would allow one to investigate the effect of error (which can be gross), temperature or whatever – play ‘what if’ games.

The final expression is automatically formed as an indirect substitution into the electroneutrality equation:

>eqn := H-OH-HCO3-2*CO3:


>H := fsolve(eqn,H,0..1):

Collect reveals that there are possibility 3 solutions; only one is realistic. It is interesting that we expect that this ‘pure rainwater’ will have a pH of potatoes or sour milk,

>pH := -log10(H);

which is just slightly acid, despite the large amount of $\ce{CO2}$.

Carbon Dioxide, Sulfur Dioxide and Ammonia

Here we include a few more species, from sulfur dioxide ($\ce{SO2}$) and ammonia ($\ce{NH3}$). Both of these constituents may have natural sources but man contributes a significant quantity of the atmospheric sulfur dioxide with smelting of ores and the burning of sulfurous coals. Ammonia, along with other gases is exuded from swamps and tip sites; animals contribute ammonia through their urine. Sulfur dioxide is an acid material and ammonia is a basic material; it is not clear in any one case whether the rainwater will be acid or basic but the acidity usually wins. The additional equilibrium constants are:

\begin{array}{ll} \ce{SO2(gas) + H2O(lqd) <=> H2SO3} & K_{hs} = \dfrac{\left[ \ce{H2SO3} \right] }{p_{\ce{CO2}}} = 124/100 \\ \ce{H+ + HSO3- <=> H2SO3} & K_{1s} = \dfrac{\left[ \ce{H2SO3} \right] }{\left[ \ce{H+}\right] \cdot \left[ \ce{HSO3-} \right] } = 79 \\ \ce{H+ + SO3^{=} <=> HSO3-} & K_{1s} = \dfrac{\left[ \ce{HSO3-} \right] }{\left[ \ce{H+}\right] \cdot \left[ \ce{SO3^{=}} \right] } = 1600000 \\ \ce{NH3(gas) + H_2O(lqd) <=> NH4OH} & K_{ha} = \dfrac{\left[ \ce{NH_4OH} \right] }{p_{\ce{NH3}}}=57 \\ \ce{NH4+ + OH- <=> NH4OH} & K_{1a} = \dfrac{\left[ \ce{NH_4OH} \right] }{\left[ \ce{NH4+} \right] \left[ \ce{OH-} \right] } = 56470 \end{array}

and the electroneutrality condition contains additional ionic species:

\begin{equation*} \ce{[ H+ ] + [ NH4+ ] ~=~ [ OH- ] + [ HCO3- ] + 2 * [ CO3^{=} ] + [ HSO3- ] + 2 * [ SO3^{=} ]} \end{equation*}

This adds the additional maple statements:

>H := 'H': # Which Redefines H from the last calculation

>pSO2 := 14/10^9: # see Scott and Hobbs[1]

>pNH3 := 7/10^9:  # parts per billion

>H2SO3 := 124/100*pSO2:

>HSO3 := H2SO3/(H*79):

>SO3 := HSO3/(H*16*10**6):

># Primary equation for calculating sulfate production

>NH4OH := 57*pNH3:

>NH4 := NH4OH/(OH*56470):

>eq := H+NH4-OH-HCO3-2*CO3-HSO3-2*SO3:

>HH := fsolve(eq,H,0..1):

>pH := evalf(-log10(HH),3):

And the pH increases in this instance, the ammonia overcoming the weak acidity of the sulfur dioxide.

Kinetics of Sulfate Formation

Now we add an extra whammy and allow oxidation of the $\ce{SO2}$ in the water droplets. This produces both acidity and sulfate ion ($\ce{SO4^{=}}$), effectively sulfuric acid, which is a strong acid; continuity of mass (mass balance) requires that, as the Sulfate is formed, the other S-containing species are depleted. The sulfate ion adds to the electroneutrality equation on the right side. Following Scott and Hobbs[1], the oxidized species is thought to be the sulfite ion ($\ce{SO3^{=}}$). The kinetic equation is the simple first-order D.E. with a right hand side that becomes complex:

\begin{equation*} \dfrac{d}{dt} \left[ \ce{SO4^{=}} \right] ~=~ k \cdot \left[ \ce{SO3^{=}} \right] \end{equation*}

The solution follows the lead above with an additional term in the electroneutrality equation:

>eqn := eq - 2*so4:

but solve still copes well with the task. It produces three solutions to the cubic expression, only one of which is realistic:

>HH := solve(eqn,H):

>H := HH[1];

The expression for the kinetic constant (turnover constant) is put in terms of the fractional turnover per day. The rate of depletion of $\ce{SO3^{=}}$ ion or the rate of increase of $\ce{SO4^{=}}$ ion is given by a complex expression with squared terms, cubic terms, square roots and cubic roots:

> k := 1/10*60*24:

> rate := k*SO3;

Continuing, we plot the rate of reaction and find that it is an ever-decreasing, hyperbolic function of the amount of $\ce{SO4^{=}}$ ion. The presence of the ammonia is, effectively, catalytic to the reaction. The rate of production eventually goes to zero, presumably because the rise in acidity decreases the $\ce{SO3^{=}}$ concentration nearly to zero.

> opt1 := title=`Rate of Sulfate Production`:

> plot(rate,SO4,0..1/1000,opt1);
> limit(rate,SO4=infinity);

Maple has some difficulty integrating the formula with all the mixed roots and powers. Probably the best approach is to use dsolve with the numeric option or the Runge-Kutta integration routine. However, such a smooth function lends itself to the use of a simple fit. Hence the solution was completed using a quadratic fit to the inverse of the rate over the range of interest. It is a more than adequate fit:

>irate := 1/rate:

>dat := [y,x]:

>for i from 0 to 10 do

> i/10000:

> evalf(evalc(subs(SO4=",irate))):

> dat := dat, [","]:


>dat := convert([dat],matrix):


>rr := regression(dat, y = a + b*x + c*x**2):


>irate1 := a + b*SO4 + c*SO4^2;


>opt1 := title = `Fitting of the inverse rates`:


The lines lie one over the other. The last operation is to integrate the form and obtain the time at which a certain amount of $\ce{SO4^{=}}$ is present in the drops. For plotting, the arguments are reversed by simply calculating a few points and reversing order. Compare the final curve with that of Scott and Hobbs[1]. If one uses precisely the environment of $\ce{CO2}$, $\ce{SO2}$ and $\ce{NH3}$ used by these authors, the results overlay the curves in the paper.

># Integrate 

>t := int(irate1,SO4=0..s);

># Calculate data points to allow reversal of plot 

>dat := NULL:

>for i from 0 to 100 do

>   i/100000:

>   evalf(subs(s=",t)):

>   dat := dat, [",96000*""]:  # conversion to micrograms/cm3 


>dat := [dat]:

>plot(dat,`Reaction time in days`=0..1,SO4=0..100);

A last calculation gives the pH variation with time. This should give a rough indication of the time required to create acid rain in the given polluted, wet environment. For this an additional term in the data preparation loop allows the $\ce{SO4^{=}}$ to become a parameter in the calculation.

>dat := NULL:

>for i from 0 to 100 do

>   i/100000:

>   evalf(subs(s=",t)):

>   -log10(evalf(subs(SO4="",H))):

>   dat := dat, ["","]:


>dat := [dat]:

>opt1 := title = ` Increasing Acidity of Rainwater`:

>plot(dat,`Reaction time in days`=0..1,pH=5.4..6.2);

The pH is now approaching the value for beer. Note that, by investigating the limiting form of the solution it is easily shown that no limit exists. That is, without some additional feedback mechanism, the pH will decrease, ever more slowly, but without a limit. Scott and Hobbs[1] thought that their model predicted a limiting pH.

This article has presented some examples of the mathematics of practical solution chemistry, as can be done with Maple. Though these examples are relatively simple, complicated structures are tractable with Maple. For Example, a recent article considers the solution chemistry of kidney stones[2]. This is but an extension of the structures given above.


This use of Maple was conceived at the 1992 Maple Summer Workshop and my outside study leave with Stan Devitt and the Symbolic Manipulation Group, Computer Sciences Department at the University of Waterloo. It became possible with the help of Professor Peter Kloeden of the Department of Computing and Mathematics at Deakin University.


1 W. D. Scott and P. V. Hobbs (1967). The formation of sulfate in water droplets, Journal of the Atmospheric Sciences, 24 (1), 54-7.

2 W. D. Scott, T. J. Wrigley and K. M. Webb (1991). A computer model of struvite solution chemistry, Talanta, 38 (8), 889-985.