Jia,
Thank you for taking a look at this simulation. I did convert to intensive units as suggested (revised input file attached) and got the same result as before as one would expect. However, I still believe the numerical model is failing to correctly simulate the physical system based on multiple lines of evidence:
1) Attached PDF page 1 shows results at 25 and 1000 years. Note the extremely high Ca++ concentration at 1000 years at position 1.5 cm (node 1). More generally I see the Ca species concentrations jumping around in this region for no apparent reason. Physically aqueous Ca should be leaving the system by diffusion out the left boundary causing slow dissolution of the minerals.
2) PDF page 2 left side shows results at 25 years when the number of nodes is increased from 10 to 50. Similar profiles are observed around the first few nodes, even through the physical spacing has changed by 5x.
3) PDF page 2 upper right shows a PHREEQC simulation of this problem. A portlandite carbonation+dissolution front has advanced from the left boundary and by 10,000 years consumed most of that mineral. A calcite dissolution front is also advancing from the left boundary at a slower rate.
4) PDF page 2 bottom right: If one knows the aqueous concentration of CO2 at the left boundary and Ca in portlandite regions, the position of the portlandite carbonation+dissolution front can be predicted using a simple analytic model. Using the CO2 (0.03454 mol/kg) and Ca (0.01495 mol/kg) concentrations from the PHREEQC simulation, the analytic solution identically matches the PHREEQC simulation for the elapsed time required to completely consume portlandite in the domain.
My sense is that dissolved Ca is not transporting by diffusion from node 1 (1.5 cm) to node 0 (0.5 cm) as it should. I would appreciate it if you could take another look at this simulation.
Greg
PortlanditeDissCarb_DIF.x1t
GWB_diagnosis.pdf