Modelling first-order kinetic reaction correctly

[Ellis]

Dear David,

I am trying to model calcite dissolution by hydrochloric acid, given by the first-order rate expression:
 r = k_(H+)*γ_(H+)*c_(H+)

My idea of solving this, is to calculate the concentration, c_(H+) at the fluid-solid interface at every time step. This works without Phreeqc; however, my aim is to solve the problem using PhreeqRM coupled with my transport solver.

The result I get so far with RM seems a bit odd, and I worry that the definitions I have set out in my input script is erroneous (for example, I don’t understand how to include c_(H+) in the RATES block since it is supposed to be an unknown). Thus, I was wondering if you could have a look at my script as shown below. Have I set out the RATES and KINETICS data blocks correctly wrt the rate expression above? If so, then how do I extract c_(H+) from the result?

Code: [Select]

TITLE Reactive Transport
SOLUTION 0  HCl
        units              mol/L
        temp             25.0
        pH                2.0     charge
        Cl                  0.01
SOLUTION 1-100  Initial solution
        units              mol/L
        temp              25.0
        pH                12.5     charge
        pe                  4
        density           1
        water             1  #kg
END
KINETICS 1-100 
        Calcite
         m0              8.5094e-06     
        parms           0.0628         #parms = surface area [cm2]
RATES
        Calcite
        -start
  10  si_cc = SI("Calcite")
  20  if (M <= 0 and si_cc < 0) then goto 200   
  30    k1   = 8.912509e-5                                 # [mol/cm2/s]
  40    area = PARM(1)
  50    rf  = k1*ACT("H+")
  60    rate  = area * rf * (1 - 10^(2/3*si_cc))
  70    moles = rate * TIME
  200 SAVE moles
        -end
END
SELECTED_OUTPUT
        file outputFile.out
        reset false
        high_precision true
USER_PUNCH
    -headings     Calcite pH
    -start
   5  PUNCH KIN("Calcite")   
  10 PUNCH -LA("H+")
    -end

Thanks for your help, comments, or suggestions.

[dlparkhurst]

γ_(H+)*c_(H+) = ACT("H+), so your rate expression seems right. If you want c_(H+), then use the MOL function, which has units of mol/kgw.

You are using the bulk solution activity of H+ in your rate equation, but you refer to the activity at the "fluid-solid interface". If you intend different activities at the interface than the bulk solution, then I am not sure how you intend to model it.

[Ellis]

Thank you for your response.

Now, after Phreeqc has performed the reaction at each time step, I understand that I will have to extract the concentration from each cell by using the code:
C = phreeqc_rm.GetConcentrations(); 

I ask:
(a) is C above equal to rate/[ k_(H+)* γ_(H+) ] ?
(b) if not, what exactly is C,
(c) if I include a “SAVE rate” in the basic program (rate obtained as in line 60 in my script above), how will I extract the saved 'rate'?
(d) What variable do I include in the MOL function to get c_(H+)? I tried several things including: Punch MOL(“c(H+)”), but that gave me error, possibly because the variable name is not recognized by Phreeqc?

PS: regarding your last statement, let me get the above cleared first, then I'll think about a way forward.

[dlparkhurst]

(a-b) The values returned from GetConcentrations are concentrations of H, O, charge, and other elements in the system. The values are not concentrations of individual aqueous species (like H+), but total concentrations of elements. The units of the concentrations are specified by the method SetSolutionUnits with options for mg/L, mol/L, or mass fraction.

It is possible to transport each aqueous species if you are modeling diffusion with methods related to SetSpeciesSaveOn and GetSpeciesConcentrations. However, for advection, you should stick with GetConcentrations, in which case, the concentration for H will be the total number of moles of H in solution--for a liter of water it is about 110 mol/L.

(c) The strategy for returning values that do not have explicit methods (explicit methods include GetConcentrations and other Get... methods) is to use PHREEQC data blocks SELECTED_OUTPUT n and USER_PUNCH n. You can retrieve values written with these data blocks with GetSelectedOutput (or other methods depending on the language you are using).

RATES are a bit tricky because they evolve and are evaluated several times over the course of a time step. I think you probably want the average rate for the time step. The following code in a USER_PUNCH data block would provide the activity of H+ (unitless), the molality of H+ (mol/kgw), H+ concentration in mol/L (for phreeqc.dat, Amm.dat, and pitzer.dat), and the average rate of the calcite reaction (mol/s):

Code: [Select]

-heading act_H+ mol_H+ c_H+ Rate_calcite
10 PUNCH ACT("H+")
20 PUNCH MOL("H+")
30 PUNCH MOL("H+") * TOT("water") / SOLN_VOL
40 PUNCH KIN_DELTA("Calcite")/KIN_TIME

(d) The Basic functions are documented in the manual in the section The Basic Interpreter. MOL("H+") provides the concentration of the aqueous species H+ in mol/kgw.