MATLAB: Model of mineral precipitation

chemistrygeochemistryiterationmineralsrate of reaction

Looking to model the precipitation of a mineral, which occurs at a rate, r (mM Kg-1 hour-1), which is dependent on the ion activity product (IAP). As the mineral precipitates, the concentration of species in solution decreases, and so the IAP also reduces and the rate slows down. I'm new to Matlab so just looking to get some basic pointers on iterative calculations ([Fe] depends on the rate of mineral precipitation and the rate depends on [Fe]) as well as modelling through time using discrete time steps (e.g. [Fe] vs time).
Code I have so far is basic chemistry:
cla;
clf;
%Basic rate calculation
Fe=0.01*0.26/1000; %mM, activity coefficient, convert to M
Si=1.8/1000; %Input mM, convert to M
pH=8;
H=10^(-pH);
k=1.11*10^(-9);
b=1.81;
Ksp=4.68*10^(22);
IAP=(Fe^(3)*Si^(2))/H^(6);
logIAP=log10(IAP); 
logKsp=log10(Ksp); 
rate=k*exp(b*(logIAP-logKsp)); %mm Fe Kg-1 hour-1
%Fe calculation
time=0.01;
Fet=Fe-(rate*time);
Thanks!

Best Answer

Hello Rosalie,
I'm not a chemist so I will take your equations as given although log10 and exp don't seem compatible somehow. Here are a couple of ways to go about this. First a rate function is defined and is at the end of the code, since when running Matlab scripts all functions have to be at the end. The rate function is shown as a function of both time and the amount of Fe. In your case the rate does not explicitly depend on time, only on Fe, but the syntax of the second method requires the time input to be there even if you don't need it. Note that the rate is negative.
I assumed that the listed value for Fe was the starting point at t=0. The first method is simple and just steps through time using the rate calculated from the previous time step, similar to what you had written. The second method uses one of the Matlab differential equation solvers which calculates its own time steps and is quite a bit more accurate. But if you zoom in on the plot you can see that the results are pretty close.
% method 1
dt = .01;                 % time step
% preallocate 12 hours worth of results, 1201 points
Fe = zeros(1,1201);       % row vector

t = zeros(1,1201);        % row vector
% initial condition 
t(1) = 0;
Fe(1) = 0.01*0.26/1000;           % mM, activity coefficient, convert to M
for n = 2:1201
  % step forward, calculate rate based on previous value and increment Fe
  t(n) = t(n-1) + dt;
  Fe(n) = Fe(n-1) + Ferate(t(n-1),Fe(n-1))*dt;     
end
% method 2
[t2 Fe2] = ode45(@(t,y) Ferate(t,y), [0 12], 0.01*0.26/1000); 
figure(1);
plot(t,Fe,t2,Fe2)
legend('method 1','method2')
function rate = Ferate(t,Fe)
Si=1.8/1000; %Input mM, convert to M
pH=8;
H=10^(-pH);
k=1.11*10^(-9);
b=1.81;
Ksp=4.68*10^(22);
IAP=(Fe.^(3)*Si^(2))/H^(6);
logIAP=log10(IAP); 
logKsp=log10(Ksp); 
rate=-k*exp(b*(logIAP-logKsp)); %mm Fe Kg-1 hour-1
end
Related Question