Am transport in groundwater

[Sen]

Dear Prof. Parkhurst

The purpose of this topic is to ask you for help with using the Transport function in PHREEQC. I also found a Topic: Dual Domain Transport of U and Cl in the forum, http://https://phreeqcusers.org/index.php/topic,1339.0.html. But I can't download any attachments and it always shows 404 - Attachment Not Found.If it is convenient for you, please send me a copy of the attachment in the link above.


In my thesis, there is a section to calculate the migration of Am in groundwater, and the model assumes that there is a repository in deep-underground, which is considered as a source of pollution and defined as Solution 0. The main rock of the repository is granite, so it is considered as column, and the whole model is simplified as possible. The groundwater in zone I is greatly influenced by the repository, which I defined as SOLUTION 1-20 in the model, and the groundwater in zone II is the original groundwater of granite, which is defined as SOLUTION 21-70 in the model, please see the attachment for the specific geochemical conditions, it should be noted that there is no Am in the initial water of zone I and zone II, Am does not exist in nature. Am is subjected to adsorption and ion exchange in granite groundwater.


In addition to the parameters used in the code, I have collected as many hydrogeological parameters as possible. transmissibility of groundwater = 8.673×10-8(m2/s)，hydraulic conductivity =5.04×10-7(m/s)，Effective diffusion coefficient of Am on granite =2×10-13(m2/s)，Pore water velocity=0.82×10-7(m/s)，Porosity= 5.8×10-3，Bulk density =3.0×106 (g/cm3)，Grams of sediment per liter (from porosity and bulk density)=5.17×105 (g/L)，Distribution coefficient of Am on granite Kd  = 0.005 m3/g.



I wonder if these parameters can make the calculation results more accurate. Here I apologize for my broken knowledge and for adding extra work to you，but I've done my best to understand the meaning of these parameters in the model.


The attached code only calculates the migration of Am in region I. I have rewritten the model according to the PHREEQC manual, and the model seems to run fine, but the plot is obviously incorrect. Can you help me to adjust the code more accurate?


Best wishes，Sen

Code: [Select]

DATABASE D:\USGS\Phreeqc Interactive 3.4.0-12927\database\wateq4f.dat
##DATABASE c:\phreeqc\database\wateq4f.dat


SOLUTION_MASTER_SPECIES
         Am           Am+3      0.0        Am        243.0000
         Am(+3)           Am+3      0.0        Am         
SOLUTION_SPECIES
         Am+3 = Am+3
            log_k      0.0
         Am+3 +  H2O - H+  = AmOH+2
        log_k      -7.200
     Am+3 + 2H2O - 2H+  = Am(OH)2+
        log_k      -15.100
     Am+3 + 3H2O - 3H+  = Am(OH)3
        log_k      -26.200
         Am+3 + F-  = AmF+2
        log_k      3.400
           delta_h    2.89 kcal
     Am+3 + 2F-  = AmF2+
        log_k      5.800
           delta_h     10.77 kcal
     Am+3 +  Cl- = AmCl+2
        log_k      0.240
     Am+3 + 2Cl- = AmCl2+
        log_k      -0.810
            delta_h    13.11 kcal
     Am+3 + SO4-2 = AmSO4+
        log_k      3.500
            delta_h    9.56 kcal
     Am+3 + 2SO4-2 = Am(SO4)2-
        log_k       5.000
            delta_h    16.72 kcal
     Am+3 + NO3- = AmNO3+2
        log_k      1.28
            delta_h    0.43 kcal
     Am+3 + 2NO3- = Am(NO3)2+
        log_k      0.88
            delta_h    2.58 kcal
     Am+3 + CO3-2 = AmCO3+
        log_k      8.000
     Am+3 + 2CO3-2 = Am(CO3)2- 
        log_k      12.900
     Am+3 + 3CO3-2 = Am(CO3)3-3 
        log_k      15.000
         Am+3 + HCO3- = AmHCO3+2
        log_k      3.100
         Am+3 +  H3SiO4- = AmSiO(OH)3+2
        log_k      8.13
           delta_h     4.44 kcal

SOLUTION_SPECIES
Na+ = Na+; -log_k 0;   -gamma 1e7   0
Ca+2 = Ca+2; -log_k 0; -gamma 1e7   0 
Cl- = Cl-; -log_k 0;   -gamma 1e7   0 
H2O + 0.01e- = H2O-0.01; -log_k -9

SURFACE_MASTER_SPECIES      #Surface Complexation of Am on Granite
        Surf    SurfOH
SURFACE_SPECIES
        SurfOH = SurfOH
                log_k   0.0
##Am
        SurfOH + Am+3 = SurfOAm+2 + H+
                log_k   -0.60
        SurfOH + Am+3 + 2H2O = SurfOAm(OH)2 + 3H+
                log_k   -17.20

EXCHANGE_MASTER_SPECIES
        X       X- 
EXCHANGE_SPECIES
        X- = X-
        log_k           0.0
        3NaX + Am+3 = AmX3 + 3Na+
        log_k           -2.2

SOLUTION 1-20   #Repository groundwater   
        temp        25
        pH          9.0
        pe          -3 
        redox       pe   
        units       mol/L
        density     1.0003     
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2   
        S           1.15e-4
        N           3.67e-4 
        C(4)        2.05e-3   CO2(g) -5.5 
        Si          7.05e-5

SURFACE 1-20
 SurfOH  7.35e-6    2.3     0.09
-donnan 1e-10 
-equil 1

EXCHANGE 1-20
X        4.04e-5
-equil 1
End

SOLUTION 0  Palse solution with Am 
        temp        25
        pH          9.0
        pe          -3 
        units       mol/L
        C(4)        0.0000001  CO2(g) -2.5 
        Si          7.05e-5
        Am          3.51e-10

SURFACE 0
SurfOH  7.35e-6   2.3   0.09
-donnan 1e-10 
-equil 0

EXCHANGE 0
X        4.04e-5
-equil 0
End

PRINT;
-reset false;
-status false

SELECTED_OUTPUT
   -file           Am-000001.sel
   -reset          false
   -solution
   -distance       true
   -time           true
   -step           true
   -pH             true
   -totals         Am
   -molalities     AmSiO(OH)3+2    AmOH+2   Am(OH)2+    AmCO3+    Am(CO3)2-   Am(CO3)3-3   SurfOAm+2   SurfOAm(OH)2   AmX3

USER_GRAPH 1
   -headings Distance Am
   -chart_title "Am transport"
   -axis_titles "Distance, in meters" "Millimoles per kilogram water"
   -axis_scale x_axis 0   auto
   -axis_scale y_axis 0   auto
   -initial_solution false
   -plot_concentration_vs x
  -start
10 x = DIST
20 PLOT_XY x, TOT("Am")*1000, symbol = Circle, color = Orange
30 PLOT_XY x, TOT("AmSiO(OH)3+2 ")*1000, symbol = Plus
40 PLOT_XY x, TOT("AmOH+2 ")*1000, symbol = Plus
50 PLOT_XY x, TOT("Am(OH)2+")*1000, symbol = Plus
60 PLOT_XY x, TOT("AmCO3+ ")*1000, symbol = Plus
70 PLOT_XY x, TOT("Am(CO3)2- ")*1000, symbol = Plus
80 PLOT_XY x, TOT("Am(CO3)3-3 ")*1000, symbol = Plus
90 PLOT_XY x, TOT("SurfOAm+2 ")*1000, symbol = Plus
100 PLOT_XY x, TOT("SurfOAm(OH)2 ")*1000, symbol = Plus
110 PLOT_XY x, TOT("AmX3 ")*1000, symbol = Plus
120 if (x > 10 OR SIM_TIME <= 0) THEN END
 -end

TRANSPORT 
-cells                20   
-length               50             ## total  length = length *  cells  # total length is 1000m
-shifts               20  1          ## total diffusion time is  shifts*time_step
-time_step            1.5768e6       ## second  3.1536e7 s= 1a   # total time is 1a

-punch_cells          20
-punch_frequency      2
-print_cells          20
-print_frequency      10

-dispersivity         0.078 
-flow_direction       forward         ##diffusion   
-boundary_cond        flux    flux    ##constant  closed   
-correct_disp         true
-diffc                2.0e-13         ##diffusion coefficient
-multi_d              true    1e-9     5.8e-3    0      1
-porosities           20*5.8e-3
END

PRINT; -user_graph false
SURFACE 0
SOLUTION 0
 temp        25
 pH          9.0
 pe          -3 
 units       mol/L
 C(4)        0.0000001    CO2(g) -2.5 
 Si          7.05e-5
 ##Am         1.27e-7
END

PRINT; -user_graph true
TRANSPORT 
-shifts 40 1
END





[dlparkhurst]

Sorry, I do not have the file that was posted.  All I can suggest is to look at example 13 in the Version 3 user's manual, which demonstrates the use of a stagnant layer.

[Sen]

Dear Prof. Parkhurst

apologize for my rush reply, Before I created this topic, I spent about half a month learning about it, but there were so many formulas that I may not have fully grasped the meaning of migration calculation, and I was very unskilled in the use of transport.

The above code is also the result of my rewrite according to Ex11-15, I hope you can take some time to help me to see the possible problems of the code.

Best wishes，
Sen

[dlparkhurst]

I suggest you start with batch or simple advection calculations. Below is a 1-cell model that allows 1000 pore volumes to move through the cell (using phreeqc.dat). If you look at the output file, you will see that all of the Am is sequestered in the exchanger; very little is dissolved; and the amount on the surface is negligible.

So, the first thing to consider is whether your exchange and surface reactions are reasonable. If so, you really need not do to much more modeling. The transport of Am will be extremely limited. If you are not satisfied with the reactions, consider how to modify the exchange and surface reactions to simulate the appropriate chemistry better.

Code: [Select]

SOLUTION_MASTER_SPECIES
         Am           Am+3      0.0        Am        243.0000
         Am(+3)           Am+3      0.0        Am         
SOLUTION_SPECIES
         Am+3 = Am+3
            log_k      0.0
         Am+3 +  H2O - H+  = AmOH+2
        log_k      -7.200
     Am+3 + 2H2O - 2H+  = Am(OH)2+
        log_k      -15.100
     Am+3 + 3H2O - 3H+  = Am(OH)3
        log_k      -26.200
         Am+3 + F-  = AmF+2
        log_k      3.400
           delta_h    2.89 kcal
     Am+3 + 2F-  = AmF2+
        log_k      5.800
           delta_h     10.77 kcal
     Am+3 +  Cl- = AmCl+2
        log_k      0.240
     Am+3 + 2Cl- = AmCl2+
        log_k      -0.810
            delta_h    13.11 kcal
     Am+3 + SO4-2 = AmSO4+
        log_k      3.500
            delta_h    9.56 kcal
     Am+3 + 2SO4-2 = Am(SO4)2-
        log_k       5.000
            delta_h    16.72 kcal
     Am+3 + NO3- = AmNO3+2
        log_k      1.28
            delta_h    0.43 kcal
     Am+3 + 2NO3- = Am(NO3)2+
        log_k      0.88
            delta_h    2.58 kcal
     Am+3 + CO3-2 = AmCO3+
        log_k      8.000
     Am+3 + 2CO3-2 = Am(CO3)2-
        log_k      12.900
     Am+3 + 3CO3-2 = Am(CO3)3-3
        log_k      15.000
         Am+3 + HCO3- = AmHCO3+2
        log_k      3.100
         Am+3 +  H3SiO4- = AmSiO(OH)3+2
        log_k      8.13
           delta_h     4.44 kcal
SURFACE_MASTER_SPECIES      #Surface Complexation of Am on Granite
        Surf    SurfOH
SURFACE_SPECIES
        SurfOH = SurfOH
                log_k   0.0
##Am
        SurfOH + Am+3 = SurfOAm+2 + H+
                log_k   -0.60
        SurfOH + Am+3 + 2H2O = SurfOAm(OH)2 + 3H+
                log_k   -17.20
EXCHANGE_SPECIES
        3NaX + Am+3 = AmX3 + 3Na+
        log_k           -2.2

SOLUTION 1   #Repository groundwater   
        temp        25
        pH          9.0
        pe          -3
        redox       pe   
        units       mol/L
        density     1.0003     
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2   
        S           1.15e-4
        N           3.67e-4
        C(4)        2.05e-3   CO2(g) -5.5
        Si          7.05e-5
END
SURFACE 1
SurfOH  7.35e-6    2.3     0.09
#-donnan 1e-10
-equil 1
END
EXCHANGE 1
X        4.04e-5
-equil 1
END
SOLUTION 0  Palse solution with Am
        temp        25
        pH          9.0
        pe          -3
        units       mol/L
        C(4)        0.0000001  CO2(g) -2.5
        Si          7.05e-5
        Am          3.51e-10
END
ADVECTION
    -cells 1
    -shifts 1000
    -punch_frequency 1000
    -warnings true
END



[Sen]

Dear Prof. Parkhurst

I am honored to have the help and guidance of a geochemist of your historical stature. I should have replied in the first place and with continued discussions with you, but I unfortunately infected with Omicron and my fever didn't go down until a few hours ago.

According to your hint, indeed as you said here the ion exchange reaction is not reasonable and the adsorption reaction is very weak. I have reviewed the literature ion exchange only works at pH < 5. The adsorption of Am by granite in alkaline environment is still not very strong, which is reasonable for the adsorption reaction in the model. I tried to modify the code and the model seems to be working, but the diagrams made by ADVECTION and TRANSPORT are very different, I think the code may not be right, can you take some time to help me?

Best wishes,
Sen

Code: [Select]

SOLUTION_MASTER_SPECIES
         Am           Am+3      0.0        Am        243.0000
         Am(+3)    Am+3      0.0        Am         
SOLUTION_SPECIES
         Am+3 = Am+3
            log_k      0.0
         Am+3 +  H2O - H+  = AmOH+2
        log_k      -7.200
     Am+3 + 2H2O - 2H+  = Am(OH)2+
        log_k      -15.100
     Am+3 + 3H2O - 3H+  = Am(OH)3
        log_k      -26.200
         Am+3 + F-  = AmF+2
        log_k      3.400
           delta_h    2.89 kcal
     Am+3 + 2F-  = AmF2+
        log_k      5.800
           delta_h     10.77 kcal
     Am+3 +  Cl- = AmCl+2
        log_k      0.240
     Am+3 + 2Cl- = AmCl2+
        log_k      -0.810
            delta_h    13.11 kcal
     Am+3 + SO4-2 = AmSO4+
        log_k      3.500
            delta_h    9.56 kcal
     Am+3 + 2SO4-2 = Am(SO4)2-
        log_k       5.000
            delta_h    16.72 kcal
     Am+3 + NO3- = AmNO3+2
        log_k      1.28
            delta_h    0.43 kcal
     Am+3 + 2NO3- = Am(NO3)2+
        log_k      0.88
            delta_h    2.58 kcal
     Am+3 + CO3-2 = AmCO3+
        log_k      8.000
     Am+3 + 2CO3-2 = Am(CO3)2- 
        log_k      12.900
     Am+3 + 3CO3-2 = Am(CO3)3-3 
        log_k      15.000
         Am+3 + HCO3- = AmHCO3+2
        log_k      3.100
         Am+3 +  H3SiO4- = AmSiO(OH)3+2
        log_k      8.13
           delta_h     4.44 kcal

SOLUTION_SPECIES
Na+ = Na+; -log_k 0;   -gamma 1e10   0
Ca+2 = Ca+2; -log_k 0; -gamma 1e10   0 
Cl- = Cl-; -log_k 0;   -gamma 1e10   0 
H2O + 0.01e- = H2O-0.01; -log_k -9

SURFACE_MASTER_SPECIES      #Surface Complexation of Am on Granite
        Surf    SurfOH
SURFACE_SPECIES
        SurfOH = SurfOH
                log_k   0.0
##Am
        SurfOH + Am+3 = SurfOAm+2 + H+
                log_k   -0.60
        SurfOH + Am+3 + 2H2O = SurfOAm(OH)2 + 3H+
                log_k   -17.20



SOLUTION 1-21   #Repository groundwater
        temp        25
        pH          9.0
        pe          -3 
        redox       pe   
        units       mol/L
        density     1.0003     
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2   
        S(6)        1.15e-4
        N(5)        3.67e-4 
        C(4)        2.05e-3   CO2(g) -5.5 
        Si          7.05e-5

SURFACE 1-21
 SurfOH  7.35e-6    2.3      0.09
-donnan 1e-10 
-equil 1
End

SOLUTION 0  Palse solution with Am 
        temp        25
        pH          9.0
        pe          -3 
        units       mol/L
        C(4)        0.0000001  CO2(g) -2.5 
        Am          1.0e-7
        Si          7.05e-5
SURFACE 0
 SurfOH  7.35e-6    2.3      0.09
-donnan 1e-10 
-equil 0
End

ADVECTION
        -cells           20
        -shifts          100
        -punch_cells     20
        -punch_frequency 1
        -print_cells     20
        -print_frequency 10
PRINT; -reset false; -status false

SELECTED_OUTPUT
        -file            Am-adv.sel
        -reset           false
        -step
        -totals          Am

USER_PUNCH
  -heading  Pore_vol
  10 PUNCH (STEP_NO + .5) / 40.
USER_GRAPH 1 Example 11
  -chart_title "Using ADVECTION Data Block"
  -headings Am
  -axis_titles "Pore volumes" "Millimoles per kilogram water"
  -axis_scale x_axis 0 auto
  -axis_scale y_axis 0 auto
  -plot_concentration_vs time
  -start
  10 x = (STEP_NO + 0.5) / cell_no
  20 PLOT_XY x, TOT("Am")*1000, symbol = None
  30 PUT(1, 1)
  -end
COPY cell 101 1-20
END
USER_GRAPH 1
        -detach
END


SOLUTION 0  Palse solution with Am 
        temp        25
        pH          9.0
        pe          -3 
        units       mol/L
        C(4)        0.0000001  CO2(g) -2.5 
        Am          1.0e-7
        Si          7.05e-5
SURFACE 0
 SurfOH  7.35e-6    2.3      0.09
-donnan 1e-10 
-equil 0
End

PRINT
   -reset false
   -echo_input true
   -status false

SELECTED_OUTPUT
   -file           Am-transport.sel
   -reset          false
   -solution
   -distance       true
   -time           true
   -step           true
   -pH             true
   -totals         Am
   -molalities     AmSiO(OH)3+2    AmOH+2   Am(OH)2+    AmCO3+    Am(CO3)2-   Am(CO3)3-3   SurfOAm+2   SurfOAm(OH)2 
USER_PUNCH
   -headings Am_mmol AmSiO(OH)3+2_mmol  AmOH+2_mmol  Am(OH)2+_mmol AmCO3+_mmol  Am(CO3)2-_mmol  Am(CO3)3-3_mmol  SurfOAm+2_mmol  SurfOAm(OH)2_mmol   
  10 PUNCH TOT("Am")*1000, TOT("AmSiO(OH)3+2")*1000,TOT("AmOH+2")*1000,TOT("AmCO3+")*1000,TOT("Am(CO3)2-")*1000,TOT("Am(CO3)3-3")*1000,TOT("SurfOAm+2")*1000,TOT("SurfOAm(OH)2")*1000,
USER_GRAPH 1
   -headings Distance Am
   -chart_title "Dual Porosity, First-Order Exchange with Implicit Mixing Factors"
   -axis_titles "Distance, in meters" "Millimoles per kilogram water"
   -axis_scale x_axis 0 auto
   -axis_scale y_axis 0 auto
   -plot_concentration_vs x
   -start
  10 GRAPH_X DIST
  20 GRAPH_Y TOT("Am")*1000
  -end
END


TRANSPORT 
-cells              20
-length             50     # total  length = length *  cells  # total length is 1000m
-shifts             20   1     # total diffusion time is  shifts*time_step
-time_step          3.1536e7   # second   1a   # total time is 10a
-dispersivity       0.078 
-punch_cells        20 
-punch_frequency    1-20 
-print              20   
-print_frequency    10   
-flow_direction     forward  # diffusion
-boundary_cond      flux    flux
-correct_disp       true
-diffc              2.0e-13    #diffusion coefficient
-multi_d            true  1e-9   5.8e-3   0    1
-porosities         20*5.8e-3
END


SOLUTION 0
 temp        25
 pH          9.0
 pe          -3 
 units       mol/L
 C(4)        0.0000001  CO2(g) -2.5 
END
PRINT; -user_graph true
TRANSPORT 
-shifts 60 1
END

[dlparkhurst]

You are not plotting the same quantities. In advection, the X axis is pore volumes (proportional to time)  and in transport the X axis is distance.

You also have not set up the same simulations for advection and transport. You should initialize both the solutions and the surfaces in cells 1-20 between the advection and the transport calculation.

I would eliminate the second SOLUTION_SPECIES data block. You are creating inconsistencies between these species and the other aqueous species (including Am species), and you will have to explain why you made just these changes to the aqueous model.

Code: [Select]

SOLUTION_MASTER_SPECIES
         Am           Am+3   0.0        Am        243.0000
         Am(+3)    Am+3      0.0        Am         
SOLUTION_SPECIES
         Am+3 = Am+3
            log_k      0.0
         Am+3 +  H2O - H+  = AmOH+2
        log_k      -7.200
     Am+3 + 2H2O - 2H+  = Am(OH)2+
        log_k      -15.100
     Am+3 + 3H2O - 3H+  = Am(OH)3
        log_k      -26.200
         Am+3 + F-  = AmF+2
        log_k      3.400
           delta_h    2.89 kcal
     Am+3 + 2F-  = AmF2+
        log_k      5.800
           delta_h     10.77 kcal
     Am+3 +  Cl- = AmCl+2
        log_k      0.240
     Am+3 + 2Cl- = AmCl2+
        log_k      -0.810
            delta_h    13.11 kcal
     Am+3 + SO4-2 = AmSO4+
        log_k      3.500
            delta_h    9.56 kcal
     Am+3 + 2SO4-2 = Am(SO4)2-
        log_k       5.000
            delta_h    16.72 kcal
     Am+3 + NO3- = AmNO3+2
        log_k      1.28
            delta_h    0.43 kcal
     Am+3 + 2NO3- = Am(NO3)2+
        log_k      0.88
            delta_h    2.58 kcal
     Am+3 + CO3-2 = AmCO3+
        log_k      8.000
     Am+3 + 2CO3-2 = Am(CO3)2-
        log_k      12.900
     Am+3 + 3CO3-2 = Am(CO3)3-3
        log_k      15.000
         Am+3 + HCO3- = AmHCO3+2
        log_k      3.100
         Am+3 +  H3SiO4- = AmSiO(OH)3+2
        log_k      8.13
           delta_h     4.44 kcal

SURFACE_MASTER_SPECIES      #Surface Complexation of Am on Granite
        Surf    SurfOH
SURFACE_SPECIES
        SurfOH = SurfOH
                log_k   0.0
##Am
        SurfOH + Am+3 = SurfOAm+2 + H+
                log_k   -0.60
        SurfOH + Am+3 + 2H2O = SurfOAm(OH)2 + 3H+
                log_k   -17.20

SOLUTION 1-21   #Repository groundwater
        temp        25
        pH          9.0
        pe          -3
        redox       pe   
        units       mol/L
        density     1.0003     
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2   
        S(6)        1.15e-4
        N(5)        3.67e-4
        C(4)        2.05e-3   CO2(g) -5.5
        Si          7.05e-5

SURFACE 1-21
SurfOH  7.35e-6    2.3      0.09
-donnan 1e-10
-equil 1
End

SOLUTION 0  Palse solution with Am
        temp        25
        pH          9.0
        pe          -3
        units       mol/L
        C(4)        0.0000001  CO2(g) -2.5
        Am          1.0e-7
        Si          7.05e-5
End

ADVECTION
        -cells           20
        -shifts          100
        -punch_cells     20
        -punch_frequency 1
        -print_cells     20
        -print_frequency 10

USER_GRAPH 1 Example 11
  -chart_title "Using ADVECTION Data Block"
  -headings Am
  -axis_titles "Pore volumes" "Millimoles per kilogram water"
  -axis_scale x_axis 0 auto
  -axis_scale y_axis 0 auto
  -plot_concentration_vs time
  -start
  10 x = (STEP_NO + 0.5) / cell_no
  20 PLOT_XY x, TOT("Am")*1000, symbol = None
  30 PUT(1, 1)
  -end

END
USER_GRAPH 1
        -detach
END


SOLUTION 0  Palse solution with Am
        temp        25
        pH          9.0
        pe          -3
        units       mol/L
        C(4)        0.0000001  CO2(g) -2.5
        Am          1.0e-7
        Si          7.05e-5
SOLUTION 1-21   #Repository groundwater
        temp        25
        pH          9.0
        pe          -3
        redox       pe   
        units       mol/L
        density     1.0003     
        Na          3.23e-2
        Ca          3.42e-4
        K           1.16e-4
        Mg          1.19e-4
        Cl          1.79e-2   
        S(6)        1.15e-4
        N(5)        3.67e-4
        C(4)        2.05e-3   CO2(g) -5.5
        Si          7.05e-5

SURFACE 1-21
SurfOH  7.35e-6    2.3      0.09
-donnan 1e-10
-equil 1
End
USER_GRAPH 1
   -headings Distance Am
   -chart_title "Dual Porosity, First-Order Exchange with Implicit Mixing Factors"
   -axis_titles "Pore Volumes" "Millimoles per kilogram water"
   -axis_scale x_axis 0 auto
   -axis_scale y_axis 0 auto
   -plot_concentration_vs x
   -start
10 GRAPH_X TOTAL_TIME/(20*3.1536e7)
  20 GRAPH_Y TOT("Am")*1000
  -end
END

TRANSPORT
-cells              20
-length             50     # total  length = length *  cells  # total length is 1000m
-shifts             100   1     # total diffusion time is  shifts*time_step
-time_step          3.1536e7   # second   1a   # total time is 10a
-dispersivity       0.078
-punch_cells        20
-punch_frequency    1-20
-print              20   
-print_frequency    10   
-flow_direction     forward  # diffusion
-boundary_cond      flux    flux
-correct_disp       true
-diffc              2.0e-13    #diffusion coefficient
-multi_d            true  1e-9   5.8e-3   0    1
-porosities         20*5.8e-3
END


SOLUTION 0
temp        25
pH          9.0
pe          -3
units       mol/L
C(4)        0.0000001  CO2(g) -2.5
END
PRINT; -user_graph true
TRANSPORT
-shifts 60 1
END

[Sen]

Dear Prof. Parkhurst
Sorry for my late Reply again. The Omicron pandemic I may have experienced recurrent infections for a short period of time and did not feel slightly better until today. Thank you very much for your selfless help, which has deepened my understanding of the model. Always your work has been very rewarding. Again, thank you for all your help!

Best wishes，
Sen