Chapter 9, Problem 21P

Recall from Sec. 8.2 that determining the chemistry of water exposed to atmospheric ${\text{CO}}_2$ can be determined by solving five nonlinear equations $\left(\text{Eqs}\text{.8}\text{.6 through 8}\text{.10}\right)$ for five unknowns: $c_T,\left[{\text{HCO}}_3^{-}\right],\left[{\text{CO}}_3^{2-}\right],\left[{\text{H}}^{+}\right],\text{ and }\left[{\text{OH}}^{-}\right]$. Employing the parameters from Sec. 8.2 and the program developed in Prob. 9.20, solve this system for conditions in 1958 when the partial pressure of ${\text{CO}}_2$ was 315 ppm. Use your results to compute the pH.

To determine

To calculate: The solutions of the system offive nonlinear equations (refer section 8.2) given by,

$\begin{array}{l}{K}_1={10}^6\frac{\left[H^{+}\right]\left[HC{O}_3{}^{-}\right]}{K_H{p}_{C{O}_2}}\\ K_2=\frac{\left[H^{+}\right]\left[C{O}_3^2{}^{-}\right]}{\left[HC{O}_3{}^{-}\right]}\\ c_T=\frac{K_H{p}_{C{O}_2}}{{10}^6}+\left[HC{O}_3{}^{-}\right]+\left[C{O}_3^2{}^{-}\right]\\ K_w=\left[H^{+}\right]\left[O{H}^{-}\right]\\ 0=\left[HC{O}_3{}^{-}\right]+2\left[C{O}_3^2{}^{-}\right]+\left[O{H}^{-}\right]-\left[H^{+}\right]\end{array}$

And, to compute the $pH$ of the rainwater.

Answer to Problem 21P

Solution: The solution to the system of five nonlinear equations is $\left[H^{+}\right]=2.3419\times {10}^{-6}$, $\left[O{H}^{-}\right]=4.2701\times {10}^{-9}$, $c_T=1.3260\times {10}^{-5}$, $\left[HC{O}_3{}^{-}\right]=2.3375\times {10}^{-6}$, and $\left[C{O}_3^2{}^{-}\right]=5.0025\times {10}^{-11}$.

The $pH$ of the rainwater is $5.6304$.

Explanation of Solution

Given Information:

The system offivesimultaneousnonlinear equations,

$\begin{array}{l}{K}_1={10}^6\frac{\left[H^{+}\right]\left[HC{O}_3{}^{-}\right]}{K_H{p}_{C{O}_2}}\\ K_2=\frac{\left[H^{+}\right]\left[C{O}_3^2{}^{-}\right]}{\left[HC{O}_3{}^{-}\right]}\end{array}$

$\begin{array}{l}{c}_T=\frac{K_H{p}_{C{O}_2}}{{10}^6}+\left[HC{O}_3{}^{-}\right]+\left[C{O}_3^2{}^{-}\right]\\ K_w=\left[H^{+}\right]\left[O{H}^{-}\right]\\ 0=\left[HC{O}_3{}^{-}\right]+2\left[C{O}_3^2{}^{-}\right]+\left[O{H}^{-}\right]-\left[H^{+}\right]\end{array}$

Here, $K_H$ is Henry’s constant, $c_T$ is total inorganic carbon, $\left[HC{O}_3{}^{-}\right]$ is bicarbonate, $K_w,K_1,\text{ and }{K}_2$ are equilibrium constants, $\left[H^{+}\right]$ is hydrogen ion, $\left[C{O}_3^2{}^{-}\right]$ is carbonate, and $\left[O{H}^{-}\right]$ is hydroxyl ion.

The values of constants $K_H,K_1,K_2,\text{ and }{K}_w$ are as follows:

$\begin{array}{l}{K}_H={10}^{-1,46}\\ K_1={10}^{-6.3}\\ K_2={10}^{-10.3}\\ K_w={10}^{-14}\end{array}$

Consider, the value of partial pressure of $C{O}_2$ in 1958 as,

$p_{C{O}_2}=315\text{ ppm}$

Formula used:

For system of n simultaneous nonlinear equations formulated as,

$\begin{array}{l}{f}_1\left(x_1,x_2,\mathrm{...},x_n\right)=0\\ f_2\left(x_1,x_2,\mathrm{...},x_n\right)=0\\ .\\ .\\ .\\ f_{n-1}\left(x_1,x_2,\mathrm{...},x_n\right)=0\\ f_n\left(x_1,x_2,\mathrm{...},x_n\right)=0\end{array}$

From Newton-Raphson method, for $m^{th}$ equation,

$f_{m,i+1}=f_{m,i}+\left(x_{1,i+1}-x_{1,i}\right)\frac{\partial f_{m,i}}{\partial x_1}+\left(x_{2,i+1}-x_{2,i}\right)\frac{\partial f_{m,i}}{\partial x_2}+\mathrm{...}+\left(x_{n,i+1}-x_{n,i}\right)\frac{\partial f_{m,i}}{\partial x_n}$

Here, $m$ represents the equation and the second subscript represents the value of the function at current value $i$ or at next value $i+1$.

Further, the above equation is rewritten as,

$-f_{m,i}+x_{1,i}\frac{\partial f_{m,i}}{\partial x_1}+x_{2,i}\frac{\partial f_{m,i}}{\partial x_2}+\mathrm{...}+x_{n,i}\frac{\partial f_{m,i}}{\partial x_n}=x_{1,i+1}\frac{\partial f_{m,i}}{\partial x_1}+x_{2,i+1}\frac{\partial f_{m,i}}{\partial x_2}+\mathrm{...}+x_{n,i+1}\frac{\partial f_{m,i}}{\partial x_n}$

Write the matrix equation for partial derivatives.

$\left[Z\right]=\left[\begin{array}{cccc}\frac{\partial f_{1,i}}{\partial x_1}& \frac{\partial f_{1,i}}{\partial x_2}& \mathrm{...}& \frac{\partial f_{1,i}}{\partial x_n}\\ \frac{\partial f_{2,i}}{\partial x_1}& \frac{\partial f_{2,i}}{\partial x_2}& \mathrm{...}& \frac{\partial f_{2,i}}{\partial x_n}\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ \frac{\partial f_{n,i}}{\partial x_1}& \frac{\partial f_{n,i}}{\partial x_2}& \mathrm{...}& \frac{\partial f_{n,i}}{\partial x_n}\end{array}\right]$

Initial and the final values are expressed as,

${\left\{X_i\right\}}^T=\left[\begin{array}{cccc}{x}_{1,i}& x_{2,i}& \mathrm{...}& x_{n,i}\end{array}\right]$

And,

${\left\{X_{i+1}\right\}}^T=\left[\begin{array}{cccc}{x}_{1,i+1}& x_{2,i+1}& \mathrm{...}& x_{n,i+1}\end{array}\right]$

Write the function values at $i$ as,

${\left\{F_i\right\}}^T=\left|\begin{array}{cccc}{f}_{1,i}& f_{2,i}& \mathrm{...}& f_{n,i}\end{array}\right|$

Thus, simplify the partial derivative equation as,

$\left[Z\right]\left\{X_{i+1}\right\}=-\left\{F_i\right\}+\left[Z\right]\left\{X_i\right\}$

The above equation can then be solved iteratively to obtain the solutions to the system of nonlinear equations.

Calculation:

Apply a technique based on Newton-Raphson method to solve the system of nonlinear equations.

Use the following MATLAB code toimplement Newton Raphson method and solve the given system of simultaneous nonlinear equations.

function Code_9_21P()

format short g

x0=[10^-7;10^-7;1e-5;2e-6;5e-11];

KH=10^-1.46;K1=10^-6.3;K2=10^-10.3;Kw=10^-14;pco2=315;

[x,f,e_a,itr]=newtmult(@jfrain,x0,[],[],KH,K1,K2,Kw,pco2);

fprintf('H: %2.4e\n',x(1))

fprintf('OH: %2.4e\n',x(2))

fprintf('cT: %2.4e\n',x(3))

fprintf('HCO3: %2.4e\n',x(4))

fprintf('CO3: %2.4e\n\n',x(5))

fprintf('Maximum relative error = %2.4g percent\n',e_a)

fprintf('Number of iterations = %d',itr)

pH=-log10(x(1));

fprintf('\npH: %2.4f\n',pH)

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[x,f,e_a,itr]=newtmult(func,x0,es,maxit,varargin)

% check for minimum input required

ifnargin<2,error('Provide at least 2 arguments for input'), end

ifnargin<3||isempty(es),es=0.0001;end

ifnargin<4||isempty(maxit),maxit=50;end

% initialize

itr=0; x=x0;

while(1)

[J,f]=func(x,varargin{:});

dx=J\f;

x=x-dx;

% increase counter in every run of the loop

itr=itr+1;

e_a=100*max(abs(dx./x));

% stopping criteria

ifitr>=maxit||e_a<=es, break, end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[J,f]=jfrain(x,varargin)

del=0.00001;

% determine the elements of J matrix

df1dx1=(f1(x(1)+del*x(1),x(2),x(3),x(4),x(5),varargin{:})-...

f1(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(1));

%

df1dx2=(f1(x(1),x(2)+del*x(2),x(3),x(4),x(5),varargin{:})-...

f1(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(2));

%

df1dx3=(f1(x(1),x(2),x(3)+del*x(3),x(4),x(5),varargin{:})-...

f1(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(3));

%

df1dx4=(f1(x(1),x(2),x(3),x(4)+del*x(4),x(5),varargin{:})-...

f1(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(4));

%

df1dx5=(f1(x(1),x(2),x(3),x(4),x(5)+del*x(5),varargin{:})-...

f1(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(5));

%

df2dx1=(f2(x(1)+del*x(1),x(2),x(3),x(4),x(5),varargin{:})-...

f2(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(1));

%

df2dx2=(f2(x(1),x(2)+del*x(2),x(3),x(4),x(5),varargin{:})-...

f2(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(2));

%

df2dx3=(f2(x(1),x(2),x(3)+del*x(3),x(4),x(5),varargin{:})-...

f2(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(3));

%

df2dx4=(f2(x(1),x(2),x(3),x(4)+del*x(4),x(5),varargin{:})-...

f2(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(4));

%

df2dx5=(f2(x(1),x(2),x(3),x(4),x(5)+del*x(5),varargin{:})-...

f2(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(5));

%

df3dx1=(f3(x(1)+del*x(1),x(2),x(3),x(4),x(5),varargin{:})-...

f3(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(1));

%

df3dx2=(f3(x(1),x(2)+del*x(2),x(3),x(4),x(5),varargin{:})-...

f3(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(2));

%

df3dx3=(f3(x(1),x(2),x(3)+del*x(3),x(4),x(5),varargin{:})-...

f3(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(3));

%

df3dx4=(f3(x(1),x(2),x(3),x(4)+del*x(4),x(5),varargin{:})-...

f3(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(4));

%

df3dx5=(f3(x(1),x(2),x(3),x(4),x(5)+del*x(5),varargin{:})-...

f3(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(5));

%

df4dx1=(f4(x(1)+del*x(1),x(2),x(3),x(4),x(5),varargin{:})-...

f4(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(1));

%

df4dx2=(f4(x(1),x(2)+del*x(2),x(3),x(4),x(5),varargin{:})-...

f4(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(2));

%

df4dx3=(f4(x(1),x(2),x(3)+del*x(3),x(4),x(5),varargin{:})-...

f4(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(3));

%

df4dx4=(f4(x(1),x(2),x(3),x(4)+del*x(4),x(5),varargin{:})-...

f4(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(4));

%

df4dx5=(f4(x(1),x(2),x(3),x(4),x(5)+del*x(5),varargin{:})-...

f4(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(5));

%

df5dx1=(f5(x(1)+del*x(1),x(2),x(3),x(4),x(5),varargin{:})-...

f5(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(1));

%

df5dx2=(f5(x(1),x(2)+del*x(2),x(3),x(4),x(5),varargin{:})-...

f5(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(2));

%

df5dx3=(f5(x(1),x(2),x(3)+del*x(3),x(4),x(5),varargin{:})-...

f5(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(3));

%

df5dx4=(f5(x(1),x(2),x(3),x(4)+del*x(4),x(5),varargin{:})-...

f5(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(4));

%

df5dx5=(f5(x(1),x(2),x(3),x(4),x(5)+del*x(5),varargin{:})-...

f5(x(1),x(2),x(3),x(4),x(5),varargin{:}))/(del*x(5));

% define the J matrix

J=[df1dx1 df1dx2 df1dx3 df1dx4 df1dx5

df2dx1 df2dx2 df2dx3 df2dx4 df2dx5

df3dx1 df3dx2 df3dx3 df3dx4 df3dx5

df4dx1 df4dx2 df4dx3 df4dx4 df4dx5

df5dx1 df5dx2 df5dx3 df5dx4 df5dx5];

% determine the elements of f vector.

f11=f1(x(1),x(2),x(3),x(4),x(5),varargin{:});

f22=f2(x(1),x(2),x(3),x(4),x(5),varargin{:});

f33=f3(x(1),x(2),x(3),x(4),x(5),varargin{:});

f44=f4(x(1),x(2),x(3),x(4),x(5),varargin{:});

f55=f5(x(1),x(2),x(3),x(4),x(5),varargin{:});

% define the elements of f vector.

f=[f11;f22;f33;f44;f55];

end

% define the five non-linear equations

% first equation

function f=f1(H,OH,cT,HCO3,CO3,KH,K1,K2,Kw,pco2)

f =1e6*H*HCO3/KH/pco2-K1;

end

% second equation

function f=f2(H,OH,cT,HCO3,CO3,KH,K1,K2,Kw,pco2)

f = H*CO3/HCO3-K2;

end

% third equation

function f=f3(H,OH,cT,HCO3,CO3,KH,K1,K2,Kw,pco2)

f = H*OH-Kw;

end

% fourth equation

function f=f4(H,OH,cT,HCO3,CO3,KH,K1,K2,Kw,pco2)

f = KH*pco2/1e6+HCO3+CO3-cT;

end

% fifth equation

function f=f5(H,OH,cT,HCO3,CO3,KH,K1,K2,Kw,pco2)

f = HCO3+2*CO3+OH-H;

end

Execute the above code to obtain the solutions as,

H: 2.3419e-06

OH: 4.2701e-09

cT: 1.3260e-05

HCO3: 2.3375e-06

CO3: 5.0025e-11

Maximum relative error = 1.265e-12 percent

Number of iterations = 6

pH: 5.6304

Hence, the solution to the system of five nonlinear equations is $\left[H^{+}\right]=2.3419\times {10}^{-6}$, $\left[O{H}^{-}\right]=4.2701\times {10}^{-9},c_T=1.3260\times {10}^{-5},\left[HC{O}_3{}^{-}\right]=2.3375\times {10}^{-6}$, $\left[C{O}_3^2{}^{-}\right]=5.0025\times {10}^{-11}$, and the $pH$ of the rainwater is $5.6304$.