Anamolous concentration with TRANSPORT in cell with flux boundary condition

[florisWQ]

Dear PhreeqcUsers,

Thank you for the helpful resource that this forum provides! Here is my first post, so let me know if there are conventions I should follow better next time (:

I've ran into an issue with the TRANSPORT function in conjunction with EQUILIBRIUM_PHASES in the script below that I hope you could help me understand. I've simplified the script I use to a column with just reacting calcite and a real but arbitrarily chosen solution (SOLUTION 1-10) composition from the literature, so this model does not reflect an actual problem. Through this defined column I inject a different solution (SOLUTION 0) that during transport causes dissolution of calcite.

SOLUTION 0 flushes the column from the first cell (0) and the new solution composition throughout the column (1-10 reaches steady state concentrations and calcite dissolution after about 20 shifts. However, the last cell (11), in which there is a flux boundary condition, the concentrations are very different from the rest of the column, for example for the Ca concentration. For Ca, the USER_GRAPH in the script below produces a graph for the last shift (as the attachment do not work at the moment, here is a an imgur link: https://imgur.com/a/21GJ3EB)

I am now wondering why the last cell behaves so differently from the rest of the column when I do have a flux BC defined. I expected that it would get the same concentrations as the rest of main column cells as these solutions would flush through the last cell. Does the last cell function as a 'bucket' that catches the outflowing fluids, or am I making a mistake in my input script for definition of transport behaviour and the BCs?

Many thanks in advance for helping out!

Code: [Select]

SOLUTION 0
temp  8
pH    4
O(0)  0.6
Na    4
S(6)  0.3
N(5)  2.5
Cl    0.75
Alkalinity 0.15
END

SOLUTION 1-10
temp  8
pH  5.8
O(0) 0.6
Na  0.7
S(6) 0.1
N(5) 0 
Cl  0.4

EQUILIBRIUM_PHASES 1-10
goethite 0 100
gibbsite -1 100
Calcite 0 100

SAVE solution 1-10

END

#set up results in graphs
PRINT
 -reset  false; -user_pr true; -status false

USER_GRAPH 1
    -headings          Ca   
    -axis_titles            "Distance" "Concentration, mmol/kgw" ""
    -chart_title            "Concentrations at shift no 1000"
-init false
-axis_scale x_axis auto
-axis_scale y_axis auto
-plot_concentration_vs x
  -start
 10 if (step_no <>1000) then goto 100
 20 x = DIST
80 PLOT_XY x, tot("Ca")*1000
100 REM end
  -end

#reactive transport step
EQUILIBRIUM_PHASES 1-10
Calcite 0 100

TRANSPORT
 -cells  10
 -shifts 1000
 -time_step 1
 -diffusion_coefficient  0.3e-9 #default value of PHREEQC
 -flow_direction  forward
 -boundary_cond   flux    flux
 -correct_disp    true
 -lengths 0.1 #1 is total length flowpath
 -dispersivities 5.37e-3

selected_output
 -file               test.xls
 -equilibrium_phases calcite goethite gibbsite
 -saturation_indices calcite CO2(g)
 -totals             N N(5) S S(6) S(-2) C(4) C(-4) Ca Mg Na K Al
 -pH
 -pe


[dlparkhurst]

I also thought that cell 11 should approach the values of cell 10. However, I wrote to Tony Appelo, and he said that cell 11 is used mostly for electro-diffusion and should be ignored in this case. You can remove it from graphing and punching with -punch_cells 1-10.

If you want the value at the end of the column, you should use the value of cell 10 with a correction of half a time step to reach the end of the column.

[florisWQ]

Many thanks, I'll use the approach that you suggest!