The 2nd parameter of -shifts and Mixing factors calculation

[Xiaodong]

I am studying the Appelo's former published paper "Multicomponent diffusion modeling in clay systems with application to the diffusion
of tritium, iodide and sodium in Opalinus Clay". I met two questions that I truly cannot answer by myself.
The fisrt question might be easy for someone, Appelo used two parameters for defining -shifts, for example,
TRANSPORT
-cells 2
-shifts 100 0;

I searched the PHREEQC manual and the whole internet. But there is no clear definition what the 2nd parameters means. Interestingly, it can truly affect the modeling results a lot. Could someone explain clearly how to define the 2nd parameter?

The second question is how to calculate the mixing factors mixf_ij when using -multi_d.
According to the manual, the mixf_ij can be calculated as,

mixf_ij = Dp * delta_t * epsilon_w * A_ij * f_bc / h_ij / V_i

The question is how to calculate A_ij in the equation? Appelo gave an equation as,

A_ij = (w_free/1000 + V_Donna) / delta_x

He mentioned we need to go to the SURFACE datablock to find the V_donna. But normally, we don't define a SURFACE datablock. Could anyone give a hint how to define the mixf_ij or A_ij in a proper way?


 

[dlparkhurst]

It is an undocumented way to set the flow direction: -1 backward, 0 diffusion, and 1 forward. It is equivalent to -direction diffusion_only.


If you are not using surfaces, then V_donna is zero. V_donna is the volume of water attached to the surface that contains the ions necessary to balance the surface charge.

[Xiaodong]

Hi dlparkhurst,

Thanks for the kind and quick reply. You answered my first question perfectly. As for the second question, I found your answer is not in agreement with the published data by Appelo in his paper "Multicomponent diffusion modeling in clay systems with application to the diffusion of tritium, iodide and sodium in Opalinus Clay". In the supporting information of the paper, Appelo listed the table below for calculating mixf_ij.

length/m   Aij / m2     Dw / (m2/s)    ε         θ2       Δt / s      hij / m     VH2O / L    V0 / L
0.5           0.0314       2.24e-9         0.16     6.25    5e5         0.025      0.126         0.8

With the data in the table, I can calculate mixf_ij with the equation,
Dw=2.24e-9
porosity=0.16
n=1
Dp=Dw*porosity**n
dt=5e5
A=0.0314
fbc=2
hij=0.025
Vi=0.001

mixf=Dp*dt*porosity*A*fbc/hij/Vi=0.072.

The calculated mixf is the same as Appelo stated, but I can not dig out how to calculate A_ij according to what you said. Could you explain in more details? Thanks a lot in advance!

[dlparkhurst]

You need to contact Tony. All the MCD work is his.

You have done a better job than I could in sorting this out. I think what I said is true, if there are surfaces, he adds more terms. My guess is that, in the absence of surfaces, all the calculations are done in a method called find_J in transport.cpp. Below are the comments to that method. Beyond this, you'll need to contact Tony directly.

Code: [Select]

find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
/* ---------------------------------------------------------------------- */
{
/* mole transfer of the individual master_species:
* Eqn 1:
* J_ij = DDt * (A_ij / lav) * (-D_i*grad(c) + D_i*z_i*c_i * SUM(D_i*z_i*grad(c)) / SUM(D_i*(z_i)^2*c_i))
* regular column, stagnant FALSE:
*    D_i = temperature-corrected Dw
*    A_ij = A_icell * A_jcell
*    A_icell = (L porewater in i_cell / length_icell) / tort_f_icell /
*       (length_icell / 2)
*    lav = A_icell + A_jcell
*    grad(c) is concentration difference in icell and jcell (dx is in lav),
for activity corrections see Appelo & Wersin, 2007.
** dec. 28, 2015**
included aq_dl in the harmonic mean:
J_ij = - b_i * b_j / (b_i + b_j) * (c_j - c_i) in mol/s (see ex 21 in the manual 3).
b_i = A1 / (G_i * h_i / 2) * Dw for a pore without EDL. A1 = aq1 / h_i (m^2).
with EDL (no aq_il_i in A1, for now):
t_aq1 = aq1 + aq_dl_i. A1 = t_aq1 / h_i. f_free_i = aq1 / t_aq1.
b_i_cat = A1 / (G_i * h_i / 2) * Dw * {f_free + (1 - f_free) * Bm}. Bm  = Boltzmann enrichment in EDL = g_dl.
b_i_ani = A1 / (G_i * h_i / 2) * Dw * {f_free + (1 - f_free) / Bm)}.
22/2/18: now calculates diffusion through EDL's of multiple, differently charged surfaces
*  stagnant TRUE:
*    same eqn for J_ij, but multiplies with 2 * mixf. (times 2, because mixf = A / (G_i * h_i))
*    mixf_ij = mixf / (Dw / init_tort_f) / new_tort_f * new_por / init_por
*    mixf is defined in MIX; Dw is default multicomponent diffusion coefficient;
*    init_tort_f equals multi_Dpor^(-multi_Dn); new_pf = new tortuosity factor.
* Interlayer diffusion (IL) takes the gradient in the equivalent concentrations on X-.
surface area A for IL:
stagnant: ct[icell].mixf_il is mixf * por_il / por.
por_il = interlayer porosity, from -interlayer_D true 'por_il'.
por = free + DL porewater porosity, from -multi_D true 'multi_Dpor'.
in regular column, A is calc'd from (free + DL porewater) and cell-length.
for IL: A * por_il / por.

por_il should be entered for the cell with the maximal cec.
IL water is related to X-, thus the cec (eq/L IL water) is the same for all cells if X is defined.
IL-water = (free + DL porewater) * por_il / por.
for IL: A * aq_il / t_aq.
*/


[Xiaodong]

Dear dlparkhurst,

I studied further Tony's paper, especially checked again how Tony described the samples. I guess I have figured out how he derived all the parameter in the table I showed above. As for the shared surface area A_ij we discussed, it was calculated simply as the cross-section of the studied column (at least this is in agreement with the all calculated numbers he showed in the table), which is reasonable.
Thank you again for your generous answers; it has been incredibly helpful to me!