Multicomponent Diffusion and Porosity in Diffusive Reactive Transport

[Baborluc]

Dear all,

I am trying to use Transport module to model reactive diffusive transport through clay with non-constant source. I am modelling the reaction as an ion exchange with the exchanger previously saturated with a set of cations. I have made a column of cells where the first and last cell do not contain any exchanger to simulate inlet and outlet reservoir of the column. On the first sight the system behaves as expected but I am not sure how the Phreeqc deals with diffusivity, porosity and geometric factor (more or less defined in Phreeqc as a parameter Archie_n). If I define porosities of each cell, does Phreeqc recalculate the diffusivity or is diffusivity supposed to be recalculated for porosity and geometric factor by hand? In that case, I suppose I would need to use different value of diffusivity in the first and last cell but I am not sure if this is possible with the Multi_D or stagnant keyword?
Also, I don't understand if I would like to look at more elements diffusion if it is actually possible to enter values of diffusivities for more elements.
Please, share any ideas how to improve such a model, thank you.

I am attaching the input file but just to show the main idea:

SOLUTION 1               #initial solution in inlet reservoir (cell 1)
   water 0.16
   units mol/l
   pH 6 charge
   pe 4 O2(g) -0.8
   C(4) 1 CO2(g) -3.5
   Na 0.094
   Cl 0.102
   Sr 0.00096
END

SOLUTION 2             #initial solution for clay saturation
   water 0.00106         #volume of 1 cell (slice) with exchanger
   units mol/l
   pH 6 charge
   pe 4 O2(g) -0.8
#   Na 0.1            #initial solution for saturation
   Na 0.168            #initial solution for saturation + amount of excheangeable Na
   Cl 0.1            #initial solution for saturation
   Ca 1.12e-01         #amount of excheangeable Ca in mol/l
   Mg 3.09e-01         #amount of excheangeable Mg in mol/l
   K 5.3e-02            #amount of excheangeable K in mol/l
   C(4) 1 CO2(g) -3.5
SAVE SOLUTION 2
END

EXCHANGE 1
   X 0.0
   -equilibrate 1
SAVE EXCHANGE 1
END

USE SOLUTION 2
EXCHANGE 2
   CaX2 1.12e-04      #amount of exchangeable Ca in mols in 1 cell (slice)
   MgX2 3.09e-04      #amount of exchangeable Mg in mols in 1 cell (slice)
   NaX 6.8e-05      #amount of exchangeable Na in mols in 1 cell (slice)
   KX 5.3e-05         #amount of exchangeable K in mols in 1 cell (slice)

   -equilibrate 2
SAVE SOLUTION 3
END

USE SOLUTION 3
EXCHANGE 3
        CaX2 1.12e-04      #amount of exchangeable Ca in mols in 1 cell (slice)
   MgX2 3.09e-04      #amount of exchangeable Mg in mols in 1 cell (slice)
   NaX 6.8e-05      #amount of exchangeable Na in mols in 1 cell (slice)
   KX 5.3e-05         #amount of exchangeable K in mols in 1 cell (slice)
   -equilibrate 3
SAVE SOLUTION 4
END

USE SOLUTION 12
USE SOLUTION 2
MIX 12            #160 ml of outlet solution without Sr, mixing with 1 volume of 1 cell/slice
   12 0.00106         
   2 0.16
SAVE SOLUTION 13      #solution for outlet reservoir (12th cell)
END

TRANSPORT                      #Diffusive transport of Sr in 10 cells of clay + 1st cell for inlet and 12th cell for outlet
   -cells 12
   -lengths  1*0.16 10*0.0015 1*0.16      #10 cells of length 1.5 mm = 15 mm in total + 1st cell for inlet and 12th cell for outlet
   -shifts 42                  
   -time_step 1 day                       #reaction time of one pore volume in one cell
   -flow_direction diffusion_only      
   -boundary_conditions closed closed      #boundary conditions on the first and last cell
   #-stagnant 1 1.0 1.0 1.0                             #I am not sure if this keyword would be of any use
   -multi_D true 1.58e-09 0.56 0.01 0.17         
   -diffusion_coefficient 1.58e-09                     #or 1.5e-10 with porosity and geometric factor applied
   -porosities 1*1.0 10*0.56 1*1.0                
   -dispersivities 1*0.0 10*0.0 1*0.0
   -correct_disp true
   -initial_time 0.0
   -punch_cells 1-12
   -punch_frequency 7
   -print_cells 1-12
   -print_frequency 7

[dlparkhurst]

If you are using diffusion only, then the dispersivity does not matter.

I think that using -multi_d, the -diffusion_coefficient setting does not matter. The diffusion coefficient in the -multi_d line applies to any aqueous species that does not have -dw defined. Definitions of -dw are in the databases phreeqc.dat, amm.dat, and pitzer.dat. It is possible to add new elements and aqueous species with SOLUTION_MASTER_SPECIES and SOLUTION_SPECIES data blocks, including -dw for the species.

I think that it should work by defining porosity.

If you want more control, I think it is possible to add MIX definitions that include the form factors, but I am not the one to help. Tony Appelo knows this part of the code. You can check the examples and his web site (hydrochemistry.eu).

[Baborluc]

Thank you for the answer.

Regardless if I enter diffusivity as a species description or in Multi_D keyword, it seems to me, that porosity and the geometric factor are not accounted for (e.g. the diffusivity is not appropriatelly lowered and therefore must be calculated by hand or obtained by fitting). I will check if this is possible to fix by MIX definitions.

[md.muniruzzaman]

Hi,

Another idea to solve you issue can be, you can assign -dw values after the correction of porosity and geometrical factors for all the species.

[Baborluc]

Just to update this post, regarding the porosities, I have understood that because along with the change in porosity there is a change in bulk density and therefore the concentration of exchangeable sites must be recalculated. Otherwise, for non-reactive tracers the porosities specified in Multi_D keyword works fine and for the inlet and outlet cell, which do not contain the exchanger, it is possible to change the porosities to 1.0 by basic function change_por.
But here arises my other question, whether it is possible to change Archie_n factor for certain cells in the same manner as porosity? Because otherwise the diffusivity in the cells without exchanger is too low and the model is inappropriate.

[dlparkhurst]

No, there is no function to update Archie's n. You should probably contact Tony Appelo, who may be able to help you.

[Baborluc]

I was trying to use -stagnant keyword to model decreasing (non-constant) concentration of a tracer in the inlet cell with the closed/closed diffusion_only conditions and I do not understand why Phreeqc keeps giving me this ERROR: Could not find immobile cell solution in TRANSPORT.

What am I missing? It is the same when I use MIX datablock instead of mixing factor defined in -stagnant.

SOLUTION 1-20                
   water 3.534e-04 l            
   Na 0.1               
   Cl 0.1               

SOLUTION 21
   water 1 l                     
   Na 0.1               
   Cl 0.1               

SOLUTION 22
   water 1 l               
   Li 0.001               
   Na 0.1               
   Cl 0.1               

EXCHANGE 1-20
   X 2.81e-04                #ion exchange sites concentration in mol (based on CEC)
   -equilibrate 1
END
   
TRANSPORT                      
   -cells 21         
   -lengths  20*0.0005 1*1e-10         
   -shifts 10                  
   -time_step 1 day               
   -flow_direction diffusion_only      
   -boundary_conditions closed closed      
   -multi_D true 1.58e-09 0.556 0.01 3.92   
   -stagnant 1 1.97e-03 0.556 1.0

[dlparkhurst]

If you use -stagnant with 21 cells, you will have to define SOLUTION 22-43.