Multicomponent diffusion and porosity update (solution volume)

[ahockin]

Hello,

I am trying to model a simple diffusion-only case for calcite dissolution using multi-component diffusion. I have a column with calcite, with 20 cells of 1 cm each. I assume the column is half water, half calcite (initial porosity = 0.5)

As the calcite dissolves, the porosity increases and I would like to use GET_POR and CHANGE_POR to calculate the volume of calcite dissolved in each cell for each timestep and then update the porosity. Two questions I have though are:

1. What initial amount of calcite do I have to specify in EQUILIBRIUM_PHASES - do I give the total amount of calcite in the whole column or the amount per cell (1 cm3 of calcite?) Or do I scale the amount to the volume of water (1000 cm3?)

2. What is the total volume in each transport cell? I read in another forum post that each transport cell is scaled to 1L, is that correct? In which case is the total volume 1000 cm3? Or is the volume of water 1000cm3 and therefore the total volume = 1000cm3 + volume_calcite?

Thank you for any help!

I have attached an example of my code here.

Code: [Select]

#DATABASE concrete_3T_V07_02.dat

TITLE Diffusion Transport, Multi-D and porosity experiment
# Simple calcite diffusion case.
# Calcite: Molar volume 36.9 cm3/mol = MW (100.09 g/mol) / rho (2.71 g/cm3)

# Model name: _multi_D_Phreeqc_Forum

SOLUTION 0 Plain water
  -units mg/L
  pH 5
  Cl 1 mg/L charge
  Na 1 mg/L

SAVE SOLUTION 0
END

SOLUTION 1-20  Initial solutions for transport column
  -units mg/L
  pH 5
  Cl 1 mg/L charge
  Na 1 mg/L

USE solution 1

EQUILIBRIUM PHASES 1-20
    Calcite 0 27 # 10 mole = aprox. 1 kg calcite

SAVE equilibrium_phases 1-20
END

PRINT; -reset false; -status false

TRANSPORT
        -cells           20
        -lengths         0.001 #m -> 1 cm sized cells, total length of column = 20*1 cm = 20cm

        -shifts          1000 #total time = shifts x time_step
        -time_step       6.21e-06 year #years
        -flow_direction  diffusion_only
        -boundary_conditions   flux closed
        -punch_cells     1 2 5 10 20
        -punch_frequency 1
        -print_cells     1 2 5 10 20
        -print_frequency 100
  -multi_D true 2.29e-9 0.5 0 1 #diffusion coeff of pure water m2/s
  -implicit true
  -porosities 0.5

SELECTED_OUTPUT 1
        -file            Phases_multi_D_Phreeqc_Forum.sel
        -reset           false
        -step
        -totals          Na Cl Ca
  -equilibrium_phases Calcite
        -high_precision true

USER_PUNCH 2
-headings TOTAL_TIME Cell Porosities calcite Diff-Ca+2
      -start
10 PUNCH TOTAL_TIME
11 PUNCH CELL_NO
20 PUNCH GET_POR(CELL_NO)
21 PUNCH EQUI("Calcite")
22 PUNCH DIFF_C("Ca+2")

100 REM calculate the porosity
140 cal = EQUI("Calcite") #Moles of  calcite in the equilibrium-phase assemblage.
141 ini_por = 0.5 #% #initial porosity, half calcite, half water
162 molmass = 100.0869 #molar mass of calcite, g/mol
164 dens = 2.711 #density of calcite, g/cm3
165 cal_V = (cal * (molmass/dens)) #Volume of calcite in cm3
166 Vtot = 1000 + cal_V #Volume water + volume of calcite, cm3
168 new_por = 1-cal_V/Vtot #new porosity (-)
169 CHANGE_POR(new_por, CELL_NO) #Changing porosity to new porosity
 -end

SELECTED_OUTPUT 2
        -file            porosities_multi_D_Phreeqc_Forum.txt
        -high_precision true
  -user_punch true


USER_GRAPH 1
 -chart_title "Concentration Calcium"
 -headings Ca_cell_1 Ca_cell_2 Ca_cell_5 Ca_cell_10 Ca_cell_20
 -axis_titles      "Time (hours)" "moles per kg water"
# -plot_concentration_vs time
 -start
 10 x = TOTAL_TIME/3600 #Total_time default is in seconds
 15 c = TOT("Ca") #*1000
# 15 c = MOL("H+") #to plot pH
# 16 c = -lOG10(c)
 20 if (cell_no = 1) then PUT(c, 10)
 30 if (cell_no = 2) then PUT(c, 20)
 40 if (cell_no = 5) then PUT(c, 30)
 50 if (cell_no = 10) then PUT(c, 40)
 55 if (cell_no = 20) then PUT(c, 50)
 60 if (cell_no <> 20) then goto 200
 70 PLOT_XY x, get(10), color = Blue , symbol_size = 2
 80 PLOT_XY x, get(20), color = Red, symbol_size = 2
 90 PLOT_XY x, get(30), color = Green, symbol_size = 2
 100 PLOT_XY x, get(40), color = cyan , symbol_size = 2
 110 PLOT_XY x, get(50), color = black , symbol_size = 2
 200 REM end
#-batch Calcium_multi_D_Phreeqc_Forum.png false false
-end


USER_GRAPH 2
 -chart_title "Calcite_Concentration"
 -headings Calcite_cell_1 Calcite_cell_2 Calcite_cell_5 Calcite_cell_10 Calcite_cell_20
 -axis_titles      "Time (hours)" "moles per kg water"
# -plot_concentration_vs time
 -start
300 x = TOTAL_TIME/3600
310 c = EQUI("Calcite") #*1000
320 if (cell_no = 1) then PUT(c, 300)
340 if (cell_no = 2) then PUT(c, 320)
350 if (cell_no = 5) then PUT(c, 340)
360 if (cell_no = 10) then PUT(c, 350)
370 if (cell_no = 20) then PUT(c, 360)
380 if (cell_no <> 20) then goto 440
390 PLOT_XY x, get(300), color = Blue , symbol_size = 2
400 PLOT_XY x, get(320), color = Red, symbol_size = 2
410 PLOT_XY x, get(340), color = Green, symbol_size = 2
420 PLOT_XY x, get(350), color = cyan , symbol_size = 2
430 PLOT_XY x, get(360), color = black , symbol_size = 2
440 REM end
#-batch Calcite_multi_D_Phreeqc_Forum.png false false
 -end


USER_GRAPH 3
 -chart_title "Porosity over time"
 -headings Porosity_cell_1 Porosity_cell_2 Porosity_cell_5 Porosity_cell_10 Porosity_cell_20
 -axis_titles      "Time (hours)" "Porosity"
# -plot_concentration_vs time
 -start
 1000 x = TOTAL_TIME/3600
 1500 c = GET_POR(CELL_NO)
 2000 if (cell_no = 1) then PUT(c, 1000)
 3000 if (cell_no = 2) then PUT(c, 2000)
 4000 if (cell_no = 5) then PUT(c, 3000)
 5000 if (cell_no = 10) then PUT(c, 4000)
 5500 if (cell_no = 20) then PUT(c, 5000)
 6000 if (cell_no <> 20) then goto 20000
 7000 PLOT_XY x, get(1000), color = Blue , symbol_size = 2
 8000 PLOT_XY x, get(2000), color = Red, symbol_size = 2
 9000 PLOT_XY x, get(3000), color = Green, symbol_size = 2
 10000 PLOT_XY x, get(4000), color = cyan , symbol_size = 2
 11000 PLOT_XY x, get(5000), color = black , symbol_size = 2
 20000 REM end
#-batch Porosity_multi_D_Phreeqc_Forum.png false false
 -end
END



[dlparkhurst]

I have delayed answering your question, because I am not sure of the answer. I am confident that the surface area of the cell for diffusion is calculated by taking the volume of solution and dividing by the cell length. So, if you have 1 liter (0.001 m^3) of water and a cell length of 0.01 m, the surface area is 0.1 m^2.

I am pretty sure that the water volume does not change as porosity is changed with CHANGE_POR, but the diffusion coefficient should change by por^n (as defined by -multi_D). I did some tests and could not convince myself that this was happening.

Anyway, here is a file from Tony Appelo that demonstrates the use of CHANGE_POR with further discussion in the paper https://www.hydrochemistry.eu/pub/appt_wapo.pdf. And the following comments:

"There is still something that should be considered in the example: the solutions
must be distributed over neighboring cells when the volume changes. In the input
file, the diffusion coefficient is adapted by the por^n factor, but PHREEQC
calculates the surface area for diffusion from the water volume / cell-length,
and thus the water volume should be adapted as well. It can be done with
SOLUTION_MIX after each step."

Implementation is left as an exercise for the user ;-).

[ahockin]

Thank you very much for the reply!

For updating the solution after each step using SOLUTION_MIX - I think this works if I were to know in advance how much the porosity will change, then it is possible to automatically update the volume by a fixed/known amount. However, as far as I can tell in our case, the change in the volume over time will also change and it is not possible to use a BASIC command to pass the volume change to the SOLUTION_MIX keyword (found by calculating the change in the volume of calcite). I think that would require something similar to CHANGE_POR, where the variable can be updated with each time step in TRANSPORT. I think a round-about way would be to update the solution volume after a set amount of transport steps, but that quickly becomes cumbersome and is also less accurate I would imagine.

I'm still struggling with the initial amount of calcite in each cell. There is a relationship between the volume of water in the solution [m3], the porosity which I specify in TRANSPORT and the cell length [m]. From what I understand, if my cell lengths increase or decrease, I also have to take this into account when specifying the initial amount of calcite in EQUILIBRIUM_PHASES, correct?

Are the moles of calcite given in EQUILIBRIUM_PHASES then really moles per meter squared, since transport is considered only in 1-dimension?

Something like:

Mole_Calcite (mole/m^2) = (1-porosity)*Soln_Vol [m^3] / Phase_Vol_Calcite [m^3/mole] / {Soln_Vol/Cell_length [m]}

Therefore, for a porosity of 0.5, Soln_Vol = 0.001 m^3, Phase_Vol_Calcite = 37 cm^3/mol (3.7e-5 m^3/mol) and a cell length of 0.001 m (1mm) -> Mole_Calcite = 13.5 moles/m^2

Again, thank you for your help!

[dlparkhurst]

If you have 1 L of water in a cell, you would need 1 liter of calcite to have porosity of 0.5. If your cell length is 0.001 m, then the surface area for the water is 0.001 m^3/0.001 m = 1 m^2. The total volume of the cell would be 2 L and the total surface area 1 m^2.

Alternatively, you could use MIX or -water to define 0.5 L of solution and use 13.5 mol of calcite in each cell.

It is possible to write out input and then run it (see example 8), although it is tedious and trying to do it sequentially would be hard. If you really want to update after each step, it would be best to use IPhreeqc and either a scripting or programming language.