Phreeqc COM coupling with a transport model (PHRRQC/MATLAB/VBA)

[ericluz]

I am having problems in MATLAB coupling phreeqc. The same code has different results in matlab (function), matlab (main), VBA and PHREEQC.

In matlab (main), only in the EQUILIBRIUM_PHASES has the same result. In VBA only with reactive transport presents different results ( TRANSPORT - KINETICS)

PHRREQC

Code: [Select]

SOLUTION 0 
       temp  30      #ºC
pH 4.61
units mg/L
Na 32684
Cl 53584
Ca 475
Alkalinity 300

END


SOLUTION 1-20
pH 8.39
temp  30 #ºC
units mg/L
Ca 337.4882396
Cl 51578.00105
Na 31680.34765
Alkalinity 0.82597427

RATES
Calcite
-start
1   REM   PARM(1) # area superficial em cm^2
2   REM   PARM(2) # exp factor

10  si_cc = SI("Calcite")
20  IF (M <= 0  and si_cc < 0) THEN GOTO 200
30  k1 = 10^(0.198 - 444.0 / TK )
40  k2 = 10^(2.84 - 2177.0 /TK )
50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)
60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )
80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME # convert from mmol to mol
200 SAVE moles
-end
 
KINETICS 1-20
Calcite  ; formula CaCO3
-m0 1.1# moles
-parms 50000 0.667
-tol   1e-8

INCREMENTAL_REACTIONS true

TRANSPORT
        -cells           20
        -lengths         0.003485 # comprimento de cada celula (m)  0.0697 m/20 = 0.003485
        -shifts          1208 # tempo de simulacao = time_step * shifts (1400*60 s/ 69.49 s) =1208
        -time_step       69.5 # (0,0003485m*0,145/ 7.27x10^-6 m/s)= 69.49 s (Lcel*porosidade/velocidade)
        -flow_direction  forward # direcao do fluxo
        -boundary_conditions   flux  flux # condicao de contorno
        -diffusion_coefficient 0.3e-9
        -dispersivities  0.002
        #-correct_disp    true
        -punch_cells      20 # define quais celulas terao o seu resultado no output


       
SELECTED_OUTPUT
        -file            tent.dat
        -reset           false
        -time
        -totals          Ca
        -high_precision true

USER_GRAPH 1 Example 11
  -headings Ca -la("H+") tot("Ca")
  -axis_titles "Tempo (min)" " mg/L" "pH"
  -plot_concentration_vs x
  10 x = total_time  / 60
  20 PLOT_XY x, TOT("Ca"), symbol = Plus, symbol_size = 2
  30 PLOT_XY x, -la("H+"), symbol = Plus, symbol_size = 2, y-axis 2


END

MATLAB (Function)

Code: [Select]

function y=APELO (par)


global  M0

Finjecao=[32684,53584,475,300];
Matriz=[31680.34765,51578.00105,337.4882396,0.82597427];
iphreeqc = actxserver('IPhreeqcCOM.Object'); %
L1 =['-m0 ',num2str(M0)];
L2 =['-parm ', num2str(par),' 0.667'];

iphreeqc.LoadDatabase('llnl.dat');
iphreeqc.ClearAccumulatedLines;

iphreeqc.AccumulateLine ('SOLUTION 0 ');
iphreeqc.AccumulateLine ('-temp 20');
iphreeqc.AccumulateLine ('-pH 4.6');
iphreeqc.AccumulateLine ('-units mg/L');
iphreeqc.AccumulateLine (['Na ' num2str(Finjecao(1))]);
iphreeqc.AccumulateLine (['Cl ' num2str(Finjecao(2))]);
iphreeqc.AccumulateLine (['Ca ' num2str(Finjecao(3))]);
iphreeqc.AccumulateLine (['Alkalinity ' num2str(Finjecao(4))]);

iphreeqc.AccumulateLine ('SOLUTION 1-40 ');
iphreeqc.AccumulateLine ('-pH 8.6');
iphreeqc.AccumulateLine ('-units mg/L');
iphreeqc.AccumulateLine (['Na ' num2str(Matriz(1))]);
iphreeqc.AccumulateLine (['Cl ' num2str(Matriz(2))]);
iphreeqc.AccumulateLine (['Ca ' num2str(Matriz(3))]);
iphreeqc.AccumulateLine (['Alkalinity ' num2str(Matriz(4))]);

iphreeqc.AccumulateLine ('RATES');
iphreeqc.AccumulateLine ('Calcite');
iphreeqc.AccumulateLine (' -start');
iphreeqc.AccumulateLine ('1   REM   PARM(1)');
iphreeqc.AccumulateLine ('2   REM   PARM(2)');
iphreeqc.AccumulateLine ('10  si_cc = SI("Calcite")');
iphreeqc.AccumulateLine ('20  IF (M <= 0  and si_cc < 0) THEN GOTO 200');
iphreeqc.AccumulateLine ('30  k1 = 10^(0.198 - 444.0 / TK )');
iphreeqc.AccumulateLine ('40  k2 = 10^(2.84 - 2177.0 /TK )');
iphreeqc.AccumulateLine ('50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)');
iphreeqc.AccumulateLine ('60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )');
iphreeqc.AccumulateLine ('80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M');
iphreeqc.AccumulateLine ('110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))');
iphreeqc.AccumulateLine ('120 rate = rate * (1 - 10^(2/3*si_cc))');
iphreeqc.AccumulateLine ('130 moles = rate * 0.001 * TIME ');
iphreeqc.AccumulateLine ('200 SAVE moles');
iphreeqc.AccumulateLine ('-end');

iphreeqc.AccumulateLine ('KINETICS 1-40');
iphreeqc.AccumulateLine ('Calcite  ; formula CaCO3');
iphreeqc.AccumulateLine (L1); %mols
iphreeqc.AccumulateLine (L2); %cm^2/mol calcite, exp factor
iphreeqc.AccumulateLine ('-tol   1e-6');
iphreeqc.AccumulateLine ('INCREMENTAL_REACTIONS true');

iphreeqc.AccumulateLine ('TRANSPORT');
iphreeqc.AccumulateLine ('-cells           40 ');
iphreeqc.AccumulateLine ('-lengths         0.0017425');
iphreeqc.AccumulateLine ('-shifts          2417');
iphreeqc.AccumulateLine ('-time_step       34.75');
iphreeqc.AccumulateLine ('-flow_direction  forward');
iphreeqc.AccumulateLine ('-boundary_conditions   fluApx  flux');
iphreeqc.AccumulateLine ('-diffusion_coefficient 0.3e-9');
iphreeqc.AccumulateLine ('-dispersivities  0.002');
iphreeqc.AccumulateLine ('-punch_cells      40');

iphreeqc.AccumulateLine ('SELECTED_OUTPUT');
iphreeqc.AccumulateLine ('-file selected.dat');
iphreeqc.AccumulateLine ('-time');
iphreeqc.AccumulateLine ('-reset           false');
iphreeqc.AccumulateLine ('-totals Ca');
iphreeqc.AccumulateLine ('-high_precision true');



try
    iphreeqc.RunAccumulated;
catch
    return;
end

out_PHREEQC = iphreeqc.GetSelectedOutputArray;


tempo=cell2mat(out_PHREEQC(2:end,5));
tempo=tempo/60;
Ca=cell2mat(out_PHREEQC(2:end,9));
y=[tempo,Ca];

end


MATLAB (main)

Code: [Select]

clc
clear all


iphreeqc = actxserver('IPhreeqcCOM.Object');
iphreeqc.LoadDatabase('llnl.dat');

iphreeqc.ClearAccumulatedLines;

iphreeqc.AccumulateLine ('SOLUTION 0 ');
iphreeqc.AccumulateLine ('-temp 20');
iphreeqc.AccumulateLine ('-pH 4.6');
iphreeqc.AccumulateLine ('-units mg/L');
iphreeqc.AccumulateLine ('Na 32684');
iphreeqc.AccumulateLine ('Cl 53584');
iphreeqc.AccumulateLine ('Ca 475');
iphreeqc.AccumulateLine ('Alkalinity 300');

iphreeqc.AccumulateLine ('SOLUTION 1-20 '); %Matriz
iphreeqc.AccumulateLine ('-temp 20');
iphreeqc.AccumulateLine ('-pH 8.6');
iphreeqc.AccumulateLine ('-units mg/L');
iphreeqc.AccumulateLine ('Na 31680.34765');
iphreeqc.AccumulateLine ('Cl 51578.00105');
iphreeqc.AccumulateLine ('Ca 337.4882396' );
iphreeqc.AccumulateLine ('Alkalinity337.4882396');

iphreeqc.AccumulateLine ('RATES');
iphreeqc.AccumulateLine ('Calcite');
iphreeqc.AccumulateLine (' -start');
iphreeqc.AccumulateLine ('1   REM   PARM(1)');
iphreeqc.AccumulateLine ('2   REM   PARM(2)');
iphreeqc.AccumulateLine ('10  si_cc = SI("Calcite")');
iphreeqc.AccumulateLine ('20  IF (M <= 0  and si_cc < 0) THEN GOTO 200');
iphreeqc.AccumulateLine ('30  k1 = 10^(0.198 - 444.0 / TK )');
iphreeqc.AccumulateLine ('40  k2 = 10^(2.84 - 2177.0 /TK )');
iphreeqc.AccumulateLine ('50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)');
iphreeqc.AccumulateLine ('60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )');
iphreeqc.AccumulateLine ('80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M');
iphreeqc.AccumulateLine ('110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))');
iphreeqc.AccumulateLine ('120 rate = rate * (1 - 10^(2/3*si_cc))');
iphreeqc.AccumulateLine ('130 moles = rate * 0.001 * TIME ');
iphreeqc.AccumulateLine ('200 SAVE moles');
iphreeqc.AccumulateLine ('-end');

iphreeqc.AccumulateLine ('KINETICS 1-20');
iphreeqc.AccumulateLine ('Calcite  ; formula CaCO3');
iphreeqc.AccumulateLine ('-m0 1.1'); %mols
iphreeqc.AccumulateLine ('-parms 50000 0.667'); %cm^2/mol calcite, exp factor
iphreeqc.AccumulateLine ('-tol   1e-8');
iphreeqc.AccumulateLine ('INCREMENTAL_REACTIONS true');

iphreeqc.AccumulateLine ('TRANSPORT');
iphreeqc.AccumulateLine ('-cells           20 ');
iphreeqc.AccumulateLine ('-lengths         0.003485');
iphreeqc.AccumulateLine ('-shifts          1208');
iphreeqc.AccumulateLine ('-time_step       69.49');
iphreeqc.AccumulateLine ('-flow_direction  forward');
iphreeqc.AccumulateLine ('-boundary_conditions   flux  flux');
iphreeqc.AccumulateLine ('-diffusion_coefficient 0.3e-9');
iphreeqc.AccumulateLine ('-dispersivities  0.002');
iphreeqc.AccumulateLine ('punch_cells      20');

iphreeqc.AccumulateLine ('SELECTED_OUTPUT');
iphreeqc.AccumulateLine ('-file selected.dat');
iphreeqc.AccumulateLine ('-time');
iphreeqc.AccumulateLine ('-totals Ca');
iphreeqc.AccumulateLine ('-high_precision true');

iphreeqc.RunAccumulated;

out_PHREEQC = iphreeqc.GetSelectedOutputArray;

tempo=cell2mat(out_PHREEQC(2:end,5));
Ca=cell2mat(out_PHREEQC(2:end,9));
y=[tempo,Ca];
plot (tempo,Ca)

I do not understand why to run phreeqc through the function the "iphreeqc.RunAccumulated " I need to be in the "try catch", in main the code runs without problem

[dlparkhurst]

You don't say what differs I would need to see the output of two runs that differ (SetOutputFileOn and related methods). Please make as simple as possible. Make sure you are using the same database in the PHREEQC calculation.

I'm not sure I understand the problem with AccumlateLine and try/catch. You need not use AccumulateLine, it is just a convenience that ensures the end-of-line characters are present. RunString or RunFile are alternatives.

[ericluz]

both methods are using the same database, the outputs have differences of values

[dlparkhurst]

IPhreeqc and PHREEQC use the same code for the same version number. There are small differences depending on the processor that is used, and there can be small differences if different convergence tolerances are used (add KNOBS; -conv 1e-12 to ensure the same tolerance is used). Again, I need to see the entire output for a simple example to be able to explain he differences.