Input file: ..\examples\ex21 Output file: ex21.out Database file: ..\database\phreeqc.dat ------------------ Reading data base. ------------------ SOLUTION_MASTER_SPECIES SOLUTION_SPECIES PHASES EXCHANGE_MASTER_SPECIES EXCHANGE_SPECIES SURFACE_MASTER_SPECIES SURFACE_SPECIES RATES END ------------------------------------ Reading input data for simulation 1. ------------------------------------ TITLE Diffusion through Opalinus Clay in a radial diffusion cell, Appelo and others, 2010, GCA, v. 74, p. 1201-1219. SOLUTION_MASTER_SPECIES Hto Hto 0.0 20 20 Na_tr Na_tr+ 0.0 22 22 Cl_tr Cl_tr- 0.0 36 36 Cs Cs+ 0.0 132.905 132.905 SOLUTION_SPECIES Hto = Hto log_k 0 gamma 1e6 0 dw 2.236e-9 Na_tr+ = Na_tr+ log_k 0 gamma 4.0 0.075 dw 1.33e-9 erm_ddl 1.23 Cl_tr- = Cl_tr- log_k 0 gamma 3.5 0.015 dw 1.31e-9 # dw = dw(water) / 1.55 = 2.03e-9 / 1.55 Cs+ = Cs+ log_k 0 gamma 3.5 0.015 dw 2.07e-9 erm_ddl 1.23 SURFACE_MASTER_SPECIES Su_fes Su_fes- # Frayed Edge Sites Su_ii Su_ii- # Type II sites of intermediate strength Su_ Su_- # Double layer, planar sites are modeled with EXCHANGE SURFACE_SPECIES Su_fes- = Su_fes- log_k 0 Na+ + Su_fes- = NaSu_fes log_k 10 Na_tr+ + Su_fes- = Na_trSu_fes log_k 10 K+ + Su_fes- = KSu_fes log_k 12.4 Cs+ + Su_fes- = CsSu_fes log_k 17.14 Su_ii- = Su_ii- log_k 0 Na+ + Su_ii- = NaSu_ii log_k 10 Na_tr+ + Su_ii- = Na_trSu_ii log_k 10 K+ + Su_ii- = KSu_ii log_k 12.1 Cs+ + Su_ii- = CsSu_ii log_k 14.6 Su_- = Su_- log_k 0 EXCHANGE_SPECIES Na_tr+ + X- = Na_trX log_k 0.0 gamma 4.0 0.075 Cs+ + X- = CsX log_k 2.04 gamma 3.5 0.015 SOLUTION 0-2 column with only cell 1, two boundary solutions 0 and 2. Na 1 Cl 1 END ----- TITLE ----- Diffusion through Opalinus Clay in a radial diffusion cell, Appelo and others, 2010, GCA, v. 74, p. 1201-1219. ------------------------------------------- Beginning of initial solution calculations. ------------------------------------------- Initial solution 0. column with only cell 1, two boundary solutions 0 and 2. -----------------------------Solution composition------------------------------ Elements Molality Moles Cl 1.000e-003 1.000e-003 Na 1.000e-003 1.000e-003 ----------------------------Description of solution---------------------------- pH = 7.000 pe = 4.000 Specific Conductance (uS/cm, 25 oC) = 123 Density (g/cm3) = 0.99708 Volume (L) = 1.00298 Activity of water = 1.000 Ionic strength = 1.000e-003 Mass of water (kg) = 1.000e+000 Total alkalinity (eq/kg) = 1.451e-009 Total carbon (mol/kg) = 0.000e+000 Total CO2 (mol/kg) = 0.000e+000 Temperature (deg C) = 25.00 Electrical balance (eq) = -1.451e-009 Percent error, 100*(Cat-|An|)/(Cat+|An|) = -0.00 Iterations = 3 Total H = 1.110124e+002 Total O = 5.550622e+001 ----------------------------Distribution of species---------------------------- Log Log Log mole V Species Molality Activity Molality Activity Gamma cm3/mol OH- 1.049e-007 1.012e-007 -6.979 -6.995 -0.016 -4.11 H+ 1.035e-007 1.000e-007 -6.985 -7.000 -0.015 0.00 H2O 5.551e+001 1.000e+000 1.744 -0.000 0.000 18.07 Cl 1.000e-003 Cl- 1.000e-003 9.649e-004 -3.000 -3.016 -0.016 18.07 H(0) 1.416e-025 H2 7.078e-026 7.079e-026 -25.150 -25.150 0.000 28.61 Na 1.000e-003 Na+ 1.000e-003 9.652e-004 -3.000 -3.015 -0.015 -1.39 NaOH 9.767e-021 9.769e-021 -20.010 -20.010 0.000 (0) O(0) 0.000e+000 O2 0.000e+000 0.000e+000 -42.080 -42.080 0.000 30.40 ------------------------------Saturation indices------------------------------- Phase SI log IAP log K(298 K, 1 atm) H2(g) -22.05 -25.15 -3.10 H2 H2O(g) -1.50 -0.00 1.50 H2O Halite -7.60 -6.03 1.57 NaCl O2(g) -39.19 -42.08 -2.89 O2 ------------------ End of simulation. ------------------ ------------------------------------ Reading input data for simulation 2. ------------------------------------ KNOBS diagonal_scale true # -tolerance 1e-20 # because of low concentrations SOLUTION 3 tracer solution pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 Hto 1.14e-6 water 0.2 SELECTED_OUTPUT file radial reset false USER_PUNCH 1 nl$ = EOL$ # newline 2 x$ = CHR$(35) # cross '#' 3 sc$ = CHR$(59) # semicolon ';' 4 pi = 2 * ARCTAN(1e10) # 3.14159... 10 height = 0.052 # length of the clay cylinder / m 20 r_int = 6.58e-3 # inner radius of clay cylinder / m 30 r_ext = 25.4e-3 # outer radius 40 thickn_filter1 = 1.8e-3 # tracer-in filter thickness / m 50 thickn_filter2 = 1.6e-3 # tracer-out filter thickness / m 60 por_filter1 = 0.418 # porosity 70 por_filter2 = 0.367 80 G_filter1 = 4.18 # geometrical factor. (for filters, G = por / 10) 90 G_filter2 = 3.67 100 V_end = 0.2 # volume of the tracer-out solution / L 110 thickn_clay = r_ext - r_int # clay thickness / m 120 por_clay = 0.159 130 rho_b_eps = 2.7 * (1 - por_clay) / por_clay # clay bulk density / porosity / (kg/L) 140 CEC = 0.12 * rho_b_eps # CEC / (eq/L porewater) 150 A_por = 37e3 * rho_b_eps # pore surface area / (m2/L porewater) 160 DIM tracer$(4), exp_time(4), scale_y1$(4), scale_y2$(4), profile_y1$(4), profile_y2$(4) 170 DATA 'Hto', 'Cl_tr', 'Na_tr', 'Cs' 180 READ tracer$(1), tracer$(2), tracer$(3), tracer$(4) 200 DATA 86400 * 20, 86400 * 40, 86400 * 45, 86400 * 1000 210 READ exp_time(1), exp_time(2), exp_time(3), exp_time(4) 230 DATA '1', '1', '1', '1' 240 READ scale_y1$(1), scale_y1$(2), scale_y1$(3), scale_y1$(4) 260 DATA '1', '1', '1', '1' 270 READ scale_y2$(1), scale_y2$(2), scale_y2$(3), scale_y2$(4) 280 DATA '0 1.2e-9', '0 2.5e-5', '0 2e-10', '0 auto' 290 READ profile_y1$(1), profile_y1$(2), profile_y1$(3), profile_y1$(4) 300 DATA '0 1.2e-9', '0 2.5e-5', '0 6e-10', '0 auto' 310 READ profile_y2$(1), profile_y2$(2), profile_y2$(3), profile_y2$(4) 350 Dw = 2.5e-9 # default tracer diffusion coefficient / (m2/s) 360 nfilt1 = 1 # number of cells in filter 1 370 nfilt2 = 1 # number of cells in filter 2 380 nclay = 11 # number of clay cells 390 f_free = 0.117 # fraction of free pore water (0.01 - 1) 400 f_DL_charge = 0.45 # fraction of CEC charge in electrical double layer 410 tort_n = -0.99 # exponent in Archie's law, -1.045 without filters 420 G_clay = por_clay^tort_n # geometrical factor 430 interlayer_D$ = 'false' # 'true' or 'false' for interlayer diffusion 440 G_IL = 700 # geometrical factor for clay interlayers 450 punch_time = 60 * 60 * 6 # punch time / seconds 460 profile$ = 'true' # 'true' or 'false' for c/x profile visualization 470 IF nfilt1 = 0 THEN thickn_filter1 = 0 480 IF nfilt2 = 0 THEN thickn_filter2 = 0 490 IF tot("Hto") > 1e-10 THEN tracer = 1 ELSE IF tot("Cl_tr") > 1e-10 THEN tracer = 2 ELSE tracer = 3 520 sol$ = nl$ + ' pH 7.6' + sc$ +' pe 14 O2(g) -1.0' + sc$ +' temp 23' 530 sol$ = sol$ + nl$ + ' Na 240' + sc$ +' K 1.61' + sc$ +' Mg 16.9' + sc$ +' Ca 25.8' + sc$ +' Sr 0.505' 540 sol$ = sol$ + nl$ + ' Cl 300' + sc$ +' S(6) 14.1' + sc$ +' Fe(2) 0.0' + sc$ +' Alkalinity 0.476' 550 tracer_phases$ = nl$ + 'PHASES ' 560 tracer_phases$ = tracer_phases$ + nl$ + ' A_Hto' + nl$ + ' Hto = Hto' + sc$ +' log_k -15' 570 tracer_phases$ = tracer_phases$ + nl$ + ' A_Na_tr' + nl$ + ' Na_trCl = Na_tr+ + Cl-' + sc$ + ' log_k -14' 580 tracer_phases$ = tracer_phases$ + nl$ + ' A_Cl_tr' + nl$ + ' NaCl_tr = Na+ + Cl_tr-' + sc$ +' log_k -14' 590 tracer_phases$ = tracer_phases$ + nl$ + ' A_Cs' + nl$ + ' CsCl = Cs+ + Cl-' + sc$ + ' log_k -13' 600 DIM tracer_equi$(4) 610 FOR i = 1 TO 4 620 tracer_equi$(i) = nl$ + 'A_' + tracer$(i) + ' 0 0' 630 NEXT i 650 punch nl$ + 'PRINT ' + sc$ + ' -reset false' + sc$ + ' -echo_input true' + sc$ + ' -user_print true' 660 IF nfilt1 = 0 THEN GOTO 800 670 punch nl$ + x$ + ' filter cells at tracer-in side...' 680 r1 = r_int - thickn_filter1 690 xf1 = thickn_filter1 / nfilt1 700 FOR i = 1 TO nfilt1 710 num$ = TRIM(STR$(i + 3)) + sc$ 720 V_water = 1e3 * height * por_filter1 * pi * (SQR(r1 + xf1) - SQR(r1)) 730 punch nl$ + 'SOLUTION ' + num$ + ' -water ' + STR$(V_water) 740 punch sol$ + nl$ 750 r1 = r1 + xf1 760 NEXT i 800 punch nl$ + nl$ + x$ + ' cells in Opalinus Clay...' 810 r1 = r_int 820 x = thickn_clay / nclay 830 FOR i = 1 TO nclay 840 num$ = TRIM(STR$(i + 3 + nfilt1)) + sc$ 850 V_water = 1e3 * height * por_clay * pi * (SQR(r1 + x) - SQR(r1)) 860 punch nl$ + 'SOLUTION ' + num$ + ' -water ' + STR$(V_water * f_free) 870 punch sol$ 880 IF f_free = 1 and tracer = 1 THEN GOTO 960 890 punch nl$ + 'SURFACE ' + num$ + ' -equil ' + num$ 900 punch nl$ + ' Su_ ' + TRIM(STR$(f_DL_charge * CEC * V_water)) + STR$(A_por) + ' ' + STR$(V_water) 910 punch nl$ + ' Su_ii ' + TRIM(STR$(7.88e-4 * rho_b_eps * V_water)) 920 punch nl$ + ' Su_fes ' + TRIM(STR$(7.4e-5 * rho_b_eps * V_water)) 930 IF f_free < 1 THEN punch nl$ + ' -Donnan ' + TRIM(STR$((1 - f_free) * 1e-3 / A_por)) 940 punch nl$ + 'EXCHANGE ' + num$ + ' -equil ' + num$ 950 punch nl$ + ' X ' + TRIM(STR$((1 - f_DL_charge) * CEC * V_water)) + nl$ 960 r1 = r1 + x 970 NEXT i 1000 IF nfilt2 = 0 THEN GOTO 1200 1010 punch nl$ + nl$ + x$ + ' tracer-out filter cells...' 1020 r1 = r_ext 1030 xf2 = thickn_filter2 / nfilt2 1040 FOR i = 1 TO nfilt2 1050 num$ = TRIM(STR$(i + 3 + nfilt1 + nclay)) + sc$ 1060 V_water = 1e3 * height * por_filter2 * pi * (SQR(r1 + xf2) - SQR(r1)) 1070 punch nl$ + 'SOLUTION ' + num$ + ' -water ' + STR$(V_water) 1080 punch sol$ + nl$ 1090 r1 = r1 + xf2 1100 NEXT i 1200 punch nl$ + x$ + ' outside solution...' 1210 num$ = TRIM(STR$(4 + nfilt1 + nclay + nfilt2)) + sc$ 1220 punch nl$ + 'SOLUTION ' + num$ + ' -water ' + STR$(V_end) 1230 punch sol$ 1240 punch nl$ + 'END' 1300 punch nl$ + tracer_phases$ 1310 punch nl$ + 'EQUILIBRIUM_PHASES ' + num$ + tracer_equi$(tracer) 1312 If tracer = 3 THEN punch nl$ + tracer_equi$(tracer + 1) 1320 punch nl$ + 'END' 1400 IF nfilt1 > 0 THEN gf1 = por_filter1 / G_filter1 1410 IF nfilt2 > 0 THEN gf2 = por_filter2 / G_filter2 1420 g = por_clay / G_clay 1500 IF nfilt1 = 0 THEN GOTO 1530 1510 r1 = r_int - thickn_filter1 1520 ff = (SQR(r1 + xf1) - SQR(r1)) * xf1 * G_filter1 / (r1 * 2) / 2 1530 ff1 = (SQR(r_int + x) - SQR(r_int)) * x * G_clay / (r_int * 2) / 2 1540 IF nfilt1 = 0 THEN ff = ff1 ELSE IF ff1 * 2 < ff THEN ff = ff1 * 2 1550 IF nfilt1 > 0 THEN ff1 = (SQR(r_int + x) - SQR(r_int)) * (xf1 / gf1 + x / g) / (2 * r_int * 2) 1560 IF nfilt1 > 0 AND ff1 < ff THEN ff = ff1 1570 IF nfilt2 > 0 THEN ff1 = (SQR(r_ext + xf2) - SQR(r_ext)) * xf2 * G_filter2 / (r_ext * 2) 1580 IF nfilt2 > 0 AND ff1 < ff THEN ff = ff1 1590 dt_max = 0.5 * ff / Dw 1610 IF punch_time < dt_max THEN dt = punch_time ELSE dt = dt_max 1620 punch_fr = 1 1630 IF dt < punch_time THEN punch_fr = ceil(punch_time / dt) 1640 dt = punch_time / punch_fr 1650 shifts = ceil(exp_time(tracer) / dt) 1700 punch nl$ + nl$ + x$ + ' mixing factors...' 1710 r1 = r_int 1720 IF nfilt1 > 0 THEN r1 = r_int - thickn_filter1 1730 A = height * 2 * pi 1740 FOR i = 0 TO nfilt1 + nclay + nfilt2 1750 IF i = 0 OR i = nfilt1 + nclay + nfilt2 THEN fbc = 2 ELSE fbc = 1 1760 IF i > nfilt1 OR nfilt1 = 0 THEN GOTO 1810 1770 IF i < nfilt1 THEN mixf = Dw * fbc / (xf1 / gf1) * dt * A * r1 / 1e-3 1780 IF i = nfilt1 THEN mixf = 2 * Dw / (xf1 / gf1 + x / g) * dt * A * r1 / 1e-3 1790 IF i < nfilt1 THEN r1 = r1 + xf1 ELSE r1 = r1 + x 1800 GOTO 1880 1810 IF i > nfilt1 + nclay THEN GOTO 1860 1820 mixf = Dw * fbc / (x / g) * dt * A * r1 / 1e-3 1830 IF i = nfilt1 + nclay AND nfilt2 > 0 THEN mixf = 2 * Dw / (xf2 / gf2 + x / g) * dt * A * r1 / 1e-3 1840 IF i < nfilt1 + nclay THEN r1 = r1 + x ELSE r1 = r1 + xf2 1850 GOTO 1880 1860 mixf = Dw * fbc / (xf2 / gf2) * dt * A * r1 / 1e-3 1870 r1 = r1 + xf2 1880 punch nl$ + 'MIX ' + TRIM(STR$(i + 3)) + sc$ + STR$(i + 4) + STR$(mixf) 1890 NEXT i 1900 punch nl$ + 'END' 2000 punch nl$ + 'TRANSPORT' 2010 stag = 2 + nfilt1 + nclay + nfilt2 2020 punch nl$ + ' -warnings true' 2030 punch nl$ + ' -shifts ' + TRIM(STR$(shifts)) 2040 punch nl$ + ' -flow diff' + sc$ + ' -cells 1' + sc$ + ' -bcon 1 2' + sc$ + ' -stag ' + TRIM(STR$(stag)) 2050 punch nl$ + ' -time ' + STR$(dt) 2060 punch nl$ + ' -multi_D true ' + STR$(Dw) + STR$(por_clay) + ' 0.0 ' + TRIM(STR$(-tort_n)) 2070 punch nl$ + ' -interlayer_D ' + interlayer_D$ + ' 0.001 0.0 ' + TRIM(STR$(G_IL)) 2080 punch nl$ + ' -punch_fr ' + TRIM(STR$(punch_fr)) + sc$ + ' -punch_c ' + TRIM(STR$(2 + stag)) 2180 FOR i = 0 to 1 2190 punch nl$ + 'USER_GRAPH ' + TRIM(STR$(tracer + i)) + ' Example 21' + nl$ 2200 punch nl$ + ' -chart_title " ' + tracer$(tracer + i) + ' Diffusion to Outer Cell"' 2210 punch nl$ + ' -plot_tsv_file ex21_' + tracer$(tracer + i) + '_rad.tsv' 2220 punch nl$ + ' -axis_scale x_axis 0 ' + TRIM(STR$(exp_time(tracer + i) / (3600 * 24))) 2230 punch nl$ + ' -axis_titles "Time, in days" "Flux, in moles per square meter per second" "Accumulated mass, in moles"' 2240 punch nl$ + ' -plot_concentration_vs time' 2250 punch nl$ + ' 10 days = total_time / (3600 * 24)' 2260 punch nl$ + ' 20 a = equi("A_' + tracer$(tracer + i) + '")' 2270 punch nl$ + ' 30 IF get(1) = 0 AND total_time > 0 THEN put(total_time, 1)' 2280 punch nl$ + ' 40 dt = get(1)' 2290 A = 2 * pi * r_ext * height 2300 i$ = TRIM(STR$(2 + i)) 2310 punch nl$ + ' 50 plot_xy days - dt / (2 * 3600 * 24), (a - get(' + i$ + ')) / dt /' + STR$(A) + ', color = Green, symbol = None' 2320 punch nl$ + ' 60 put(a, ' + i$ + ')' 2330 punch nl$ + ' 70 plot_xy days, equi("A_' + tracer$(tracer + i) + '"), y_axis = 2, color = Red, symbol = None' 2340 IF tracer < 3 THEN GOTO 2360 2350 NEXT i 2360 punch nl$ + 'END' 2400 IF profile$ = 'true' THEN GOSUB 3000 2410 IF tracer < 3 THEN END # finished for Hto and Cl 2420 IF profile$ = 'false' THEN punch nl$ + 'USER_GRAPH ' + TRIM(STR$(tracer)) + sc$ + ' -detach' ELSE punch nl$ + 'USER_GRAPH ' + TRIM(STR$(tracer + 4)) + sc$ + ' -detach' 2440 tracer = tracer + 1 2450 punch nl$ + 'TRANSPORT' 2460 shifts = ceil((exp_time(tracer) - exp_time(tracer - 1))/ dt) 2480 punch nl$ + ' -shifts ' + TRIM(STR$(shifts)) 2490 punch nl$ + ' -punch_fr ' + TRIM(STR$(punch_fr)) + sc$ + ' -punch_c ' + TRIM(STR$(2 + stag)) 2500 punch nl$ + 'END' 2510 IF profile$ = 'true' THEN GOSUB 3000 2520 END # finished... 3000 punch nl$ + 'TRANSPORT' 3010 punch nl$ + ' -shifts 0' 3020 punch nl$ + ' -punch_fr 2' + sc$ + ' -punch_c 3-' + TRIM(STR$(2 + stag)) 3030 punch nl$ + 'USER_GRAPH ' + TRIM(STR$(tracer)) + sc$ + ' -detach' 3040 punch nl$ + 'USER_GRAPH ' + TRIM(STR$(tracer + 4)) + ' Example 21' + nl$ 3050 punch nl$ + ' -chart_title "' + tracer$(tracer) + ' Concentration Profile: Filter1 | Clay | Filter2"' 3060 REM punch nl$ + ' -plot_tsv_file + tracer$(tracer) + '_prof.tsv' 3070 punch nl$ + ' -axis_scale x_axis 0 ' + TRIM(STR$((thickn_filter1 + thickn_clay + thickn_filter2) * 1e3)) 3080 punch nl$ + ' -axis_scale y_axis ' + profile_y1$(tracer) 3090 punch nl$ + ' -axis_scale sy_axis ' + profile_y2$(tracer) 3100 punch nl$ + ' -axis_titles ' + '"Distance, in millimeters" "Free pore-water molality" "Total molality"' 3110 punch nl$ + ' -headings ' + tracer$(tracer) + '_free ' + tracer$(tracer) + '_tot' 3120 punch nl$ + ' -plot_concentration_vs x' 3130 punch nl$ + ' -initial_solutions true' 3140 punch nl$ + ' 10 IF cell_no = 3 THEN xval = 0 ELSE xval = get(14)' 3150 punch nl$ + ' 20 IF (' + TRIM(STR$(nfilt1)) + ' = 0 OR cell_no > ' + TRIM(STR$(nfilt1 + 3)) + ') THEN GOTO 60' 3160 punch nl$ + ' 30 IF (cell_no = 4) THEN xval = xval + 0.5 * ' + TRIM(STR$(xf1)) 3170 punch nl$ + ' 40 IF (cell_no > 4 AND cell_no < ' + TRIM(STR$(nfilt1 + 4)) + ') THEN xval = xval + ' + TRIM(STR$(xf1)) 3180 punch nl$ + ' 50 GOTO 200' 3190 punch nl$ + ' 60 IF (cell_no = ' + TRIM(STR$(4 + nfilt1)) + ') THEN xval = xval + 0.5 * ' + TRIM(STR$(xf1)) + ' + 0.5 * ' + TRIM(STR$(x)) 3200 punch nl$ + ' 70 IF (cell_no > ' + TRIM(STR$(4 + nfilt1)) + ' AND cell_no < ' + TRIM(STR$(4 + nfilt1 + nclay)) + ') THEN xval = xval + ' + TRIM(STR$(x)) + ' ELSE GOTO 90' 3210 punch nl$ + ' 80 GOTO 200' 3220 punch nl$ + ' 90 IF (cell_no = ' + TRIM(STR$(4 + nfilt1 + nclay)) + ') THEN xval = xval + 0.5 * ' + TRIM(STR$(x)) + ' + 0.5 * ' + TRIM(STR$(xf2)) 3230 punch nl$ + ' 100 IF (cell_no > ' + TRIM(STR$(4 + nfilt1 + nclay)) + ' AND cell_no <= ' + TRIM(STR$(3 + nfilt1 + nclay + nfilt2)) + ') THEN xval = xval + ' + TRIM(STR$(xf2)) 3240 punch nl$ + ' 110 IF (cell_no = ' + TRIM(STR$(4 + nfilt1 + nclay + nfilt2)) + ') THEN xval = xval + 0.5 * ' + TRIM(STR$(xf2)) 3250 punch nl$ + ' 200 y1 = TOT("' + tracer$(tracer) + '")' 3260 punch nl$ + ' 210 plot_xy xval * 1e3, y1, color = Blue, symbol = Plus' 3270 punch nl$ + ' 220 IF cell_no = 3 THEN put(y1, 15)' 3280 punch nl$ + ' 230 IF (cell_no < ' + TRIM(STR$(4 + nfilt1)) + ' OR cell_no > ' + TRIM(STR$(3 + nfilt1 + nclay)) + ') THEN GOTO 400' 3290 punch nl$ + ' 240 y2 = SYS("' + tracer$(tracer) + '") / (tot("water") + edl("water"))' 3310 punch nl$ + ' 250 REM y2 = y2 / ' + TRIM(STR$(rho_b_eps)) + x$ + ' conc / kg solid' 3320 punch nl$ + ' 260 plot_xy xval * 1e3, y2, symbol = Circle, y_axis = 2' 3330 punch nl$ + ' 270 IF (cell_no > ' + TRIM(STR$(5 + nfilt1)) + ') THEN GOTO 400' 3340 punch nl$ + ' 280 IF ' + TRIM(STR$(nfilt1)) + ' THEN plot_xy ' + TRIM(STR$(thickn_filter1 * 1e3)) + ', get(15), color = Black, symbol = None' 3350 punch nl$ + ' 290 IF ' + TRIM(STR$(nfilt2)) + ' THEN plot_xy ' + TRIM(STR$((thickn_filter1 + thickn_clay) * 1e3)) + ', get(15), color = Black, symbol = None' 3360 punch nl$ + ' 300 put(0, 15)' 3370 punch nl$ + ' 400 put(xval, 14)' 3380 punch nl$ + 'END' 3390 RETURN END ------------------------------------------- Beginning of initial solution calculations. ------------------------------------------- Initial solution 3. tracer solution WARNING: USER_PUNCH: Headings count doesn't match number of calls to PUNCH. -----------------------------Solution composition------------------------------ Elements Molality Moles Alkalinity 4.760e-004 9.520e-005 Ca 2.580e-002 5.160e-003 Cl 3.000e-001 6.000e-002 Hto 1.140e-009 2.280e-010 K 1.610e-003 3.220e-004 Mg 1.690e-002 3.380e-003 Na 2.400e-001 4.800e-002 S(6) 1.410e-002 2.820e-003 Sr 5.050e-004 1.010e-004 ----------------------------Description of solution---------------------------- pH = 7.600 pe = 13.120 Equilibrium with O2(g) Specific Conductance (uS/cm, 23 oC) = 29877 Density (g/cm3) = 1.01169 Volume (L) = 0.20165 Activity of water = 0.990 Ionic strength = 3.657e-001 Mass of water (kg) = 2.000e-001 Total carbon (mol/kg) = 4.793e-004 Total CO2 (mol/kg) = 4.793e-004 Temperature (deg C) = 23.00 Electrical balance (eq) = -1.312e-004 Percent error, 100*(Cat-|An|)/(Cat+|An|) = -0.10 Iterations = 6 Total H = 2.220258e+001 Total O = 1.111286e+001 ----------------------------Distribution of species---------------------------- Log Log Log mole V Species Molality Activity Molality Activity Gamma cm3/mol OH- 5.192e-007 3.419e-007 -6.285 -6.466 -0.181 -3.26 H+ 3.238e-008 2.512e-008 -7.490 -7.600 -0.110 0.00 H2O 5.551e+001 9.899e-001 1.744 -0.004 0.000 18.06 C(4) 4.793e-004 HCO3- 3.229e-004 2.295e-004 -3.491 -3.639 -0.148 26.37 NaHCO3 9.167e-005 9.972e-005 -4.038 -4.001 0.037 19.41 CaHCO3+ 2.585e-005 1.871e-005 -4.588 -4.728 -0.140 9.84 MgHCO3+ 1.824e-005 1.230e-005 -4.739 -4.910 -0.171 5.70 CO2 1.235e-005 1.343e-005 -4.908 -4.872 0.037 30.80 CaCO3 4.029e-006 4.383e-006 -5.395 -5.358 0.037 -14.61 MgCO3 1.615e-006 1.757e-006 -5.792 -5.755 0.037 -17.09 CO3-2 1.612e-006 4.110e-007 -5.793 -6.386 -0.593 -1.89 SrHCO3+ 5.910e-007 4.200e-007 -6.228 -6.377 -0.148 (0) NaCO3- 4.415e-007 3.137e-007 -6.355 -6.503 -0.148 0.60 SrCO3 2.913e-008 3.169e-008 -7.536 -7.499 0.037 -14.14 Ca 2.580e-002 Ca+2 2.378e-002 6.598e-003 -1.624 -2.181 -0.557 -17.02 CaSO4 1.991e-003 2.166e-003 -2.701 -2.664 0.037 7.42 CaHCO3+ 2.585e-005 1.871e-005 -4.588 -4.728 -0.140 9.84 CaCO3 4.029e-006 4.383e-006 -5.395 -5.358 0.037 -14.61 CaOH+ 5.899e-008 4.315e-008 -7.229 -7.365 -0.136 (0) CaHSO4+ 4.754e-010 3.477e-010 -9.323 -9.459 -0.136 (0) Cl 3.000e-001 Cl- 3.000e-001 2.017e-001 -0.523 -0.695 -0.172 18.53 H(0) 0.000e+000 H2 0.000e+000 0.000e+000 -44.617 -44.580 0.037 28.61 Hto 1.140e-009 Hto 1.140e-009 1.140e-009 -8.943 -8.943 0.000 (0) K 1.610e-003 K+ 1.591e-003 1.061e-003 -2.798 -2.974 -0.176 9.40 KSO4- 1.899e-005 1.349e-005 -4.722 -4.870 -0.148 (0) Mg 1.690e-002 Mg+2 1.511e-002 4.618e-003 -1.821 -2.336 -0.515 -20.64 MgSO4 1.771e-003 1.927e-003 -2.752 -2.715 0.037 5.76 MgHCO3+ 1.824e-005 1.230e-005 -4.739 -4.910 -0.171 5.70 MgCO3 1.615e-006 1.757e-006 -5.792 -5.755 0.037 -17.09 MgOH+ 7.505e-007 5.510e-007 -6.125 -6.259 -0.134 (0) Na 2.400e-001 Na+ 2.377e-001 1.722e-001 -0.624 -0.764 -0.140 -0.90 NaSO4- 2.248e-003 1.597e-003 -2.648 -2.797 -0.148 21.10 NaHCO3 9.167e-005 9.972e-005 -4.038 -4.001 0.037 19.41 NaCO3- 4.415e-007 3.137e-007 -6.355 -6.503 -0.148 0.60 NaOH 5.413e-018 5.889e-018 -17.267 -17.230 0.037 (0) O(0) 2.437e-004 O2 1.218e-004 1.325e-004 -3.914 -3.878 0.037 30.24 S(6) 1.410e-002 SO4-2 8.029e-003 1.874e-003 -2.095 -2.727 -0.632 15.89 NaSO4- 2.248e-003 1.597e-003 -2.648 -2.797 -0.148 21.10 CaSO4 1.991e-003 2.166e-003 -2.701 -2.664 0.037 7.42 MgSO4 1.771e-003 1.927e-003 -2.752 -2.715 0.037 5.76 SrSO4 4.204e-005 4.573e-005 -4.376 -4.340 0.037 24.16 KSO4- 1.899e-005 1.349e-005 -4.722 -4.870 -0.148 (0) HSO4- 5.993e-009 4.384e-009 -8.222 -8.358 -0.136 40.64 CaHSO4+ 4.754e-010 3.477e-010 -9.323 -9.459 -0.136 (0) Sr 5.050e-004 Sr+2 4.623e-004 1.281e-004 -3.335 -3.892 -0.557 -16.73 SrSO4 4.204e-005 4.573e-005 -4.376 -4.340 0.037 24.16 SrHCO3+ 5.910e-007 4.200e-007 -6.228 -6.377 -0.148 (0) SrCO3 2.913e-008 3.169e-008 -7.536 -7.499 0.037 -14.14 SrOH+ 3.694e-010 2.590e-010 -9.432 -9.587 -0.154 (0) ------------------------------Saturation indices------------------------------- Phase SI log IAP log K(296 K, 1 atm) Anhydrite -0.65 -4.91 -4.26 CaSO4 Aragonite -0.24 -8.57 -8.32 CaCO3 Calcite -0.10 -8.57 -8.47 CaCO3 Celestite 0.03 -6.62 -6.65 SrSO4 CO2(g) -3.43 -4.87 -1.44 CO2 Dolomite -0.25 -17.29 -17.04 CaMg(CO3)2 Gypsum -0.34 -4.92 -4.58 CaSO4:2H2O H2(g) -41.48 -44.58 -3.10 H2 H2O(g) -1.56 -0.00 1.55 H2O Halite -3.03 -1.46 1.57 NaCl O2(g) -1.00 -3.88 -2.88 O2 Pressure 0.1 atm, phi 1.000. Strontianite -1.01 -10.28 -9.27 SrCO3 Sylvite -4.55 -3.67 0.88 KCl ------------------ End of simulation. ------------------ ------------------------------------ Reading input data for simulation 3. ------------------------------------ PRINT selected_output false status false PRINT reset false user_print true SOLUTION 4 water 1.3963e-003 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SOLUTION 5 water 7.7322e-005 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 5 equilibrate 5 Su_ 5.0965e-004 5.2840e+005 6.6087e-004 Su_ii 7.4371e-006 Su_fes 6.9841e-007 donnan 1.6711e-009 EXCHANGE 5 equilibrate 5 X 6.2290e-004 SOLUTION 6 water 9.5113e-005 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 6 equilibrate 6 Su_ 6.2692e-004 5.2840e+005 8.1293e-004 Su_ii 9.1484e-006 Su_fes 8.5911e-007 donnan 1.6711e-009 EXCHANGE 6 equilibrate 6 X 7.6624e-004 SOLUTION 7 water 1.1291e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 7 equilibrate 7 Su_ 7.4419e-004 5.2840e+005 9.6500e-004 Su_ii 1.0860e-005 Su_fes 1.0198e-006 donnan 1.6711e-009 EXCHANGE 7 equilibrate 7 X 9.0957e-004 SOLUTION 8 water 1.3070e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 8 equilibrate 8 Su_ 8.6146e-004 5.2840e+005 1.1171e-003 Su_ii 1.2571e-005 Su_fes 1.1805e-006 donnan 1.6711e-009 EXCHANGE 8 equilibrate 8 X 1.0529e-003 SOLUTION 9 water 1.4849e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 9 equilibrate 9 Su_ 9.7873e-004 5.2840e+005 1.2691e-003 Su_ii 1.4282e-005 Su_fes 1.3412e-006 donnan 1.6711e-009 EXCHANGE 9 equilibrate 9 X 1.1962e-003 SOLUTION 10 water 1.6628e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 10 equilibrate 10 Su_ 1.0960e-003 5.2840e+005 1.4212e-003 Su_ii 1.5994e-005 Su_fes 1.5019e-006 donnan 1.6711e-009 EXCHANGE 10 equilibrate 10 X 1.3396e-003 SOLUTION 11 water 1.8407e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 11 equilibrate 11 Su_ 1.2133e-003 5.2840e+005 1.5733e-003 Su_ii 1.7705e-005 Su_fes 1.6626e-006 donnan 1.6711e-009 EXCHANGE 11 equilibrate 11 X 1.4829e-003 SOLUTION 12 water 2.0186e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 12 equilibrate 12 Su_ 1.3305e-003 5.2840e+005 1.7253e-003 Su_ii 1.9416e-005 Su_fes 1.8233e-006 donnan 1.6711e-009 EXCHANGE 12 equilibrate 12 X 1.6262e-003 SOLUTION 13 water 2.1966e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 13 equilibrate 13 Su_ 1.4478e-003 5.2840e+005 1.8774e-003 Su_ii 2.1127e-005 Su_fes 1.9840e-006 donnan 1.6711e-009 EXCHANGE 13 equilibrate 13 X 1.7696e-003 SOLUTION 14 water 2.3745e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 14 equilibrate 14 Su_ 1.5651e-003 5.2840e+005 2.0295e-003 Su_ii 2.2839e-005 Su_fes 2.1448e-006 donnan 1.6711e-009 EXCHANGE 14 equilibrate 14 X 1.9129e-003 SOLUTION 15 water 2.5524e-004 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SURFACE 15 equilibrate 15 Su_ 1.6824e-003 5.2840e+005 2.1815e-003 Su_ii 2.4550e-005 Su_fes 2.3055e-006 donnan 1.6711e-009 EXCHANGE 15 equilibrate 15 X 2.0562e-003 SOLUTION 16 water 5.0266e-003 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 SOLUTION 17 water 2.0000e-001 pH 7.6 pe 14 O2(g) -1.0 temp 23 Na 240 K 1.61 Mg 16.9 Ca 25.8 Sr 0.505 Cl 300 S(6) 14.1 Fe(2) 0.0 Alkalinity 0.476 END PHASES A_Hto Hto = Hto log_k -15 A_Na_tr Na_trCl = Na_tr+ + Cl- log_k -14 A_Cl_tr NaCl_tr = Na+ + Cl_tr- log_k -14 A_Cs CsCl = Cs+ + Cl- log_k -13 EQUILIBRIUM_PHASES 17 A_Hto 0 0 END MIX 3 4 6.6932e-004 MIX 4 5 1.9640e-004 MIX 5 6 1.5725e-004 MIX 6 7 1.8971e-004 MIX 7 8 2.2216e-004 MIX 8 9 2.5461e-004 MIX 9 10 2.8706e-004 MIX 10 11 3.1951e-004 MIX 11 12 3.5196e-004 MIX 12 13 3.8441e-004 MIX 13 14 4.1686e-004 MIX 14 15 4.4931e-004 MIX 15 16 7.7653e-004 MIX 16 17 4.2533e-003 END TRANSPORT warnings true shifts 1120 flow_direction diff cells 1 bcond 1 2 stagnant 15 timest 1.5429e+003 multi_d true 2.5000e-009 1.5900e-001 0.0 9.9000e-001 interlayer_d false 0.001 0.0 700 punch_frequency 14 punch_cells 17 USER_GRAPH 1 Example 21 WARNING: No cell-lengths were read; length = 1 m assumed. WARNING: No dispersivities were read; disp = 0 assumed. -chart_title " Hto Diffusion to Outer Cell" -plot_tsv_file ex21_Hto_rad.tsv -axis_scale x_axis 0 20 -axis_titles "Time, in days" "Flux, in moles per square meter per second" "Accumulated mass, in moles" -plot_concentration_vs time 10 days = total_time / (3600 * 24) 20 a = equi("A_Hto") 30 IF get(1) = 0 AND total_time > 0 THEN put(total_time, 1) 40 dt = get(1) 50 plot_xy days - dt / (2 * 3600 * 24), (a - get(2)) / dt / 8.2988e-003, color = Green, symbol = None 60 put(a, 2) 70 plot_xy days, equi("A_Hto"), y_axis = 2, color = Red, symbol = None END WARNING: Maximum iterations exceeded, 100 WARNING: Numerical method failed with this set of convergence parameters. WARNING: Trying smaller step size, pe step size 10, 5 ... WARNING: Maximum iterations exceeded, 200 WARNING: Numerical method failed with this set of convergence parameters. WARNING: Trying reduced tolerance 1e-016 ... TRANSPORT shifts 0 punch_frequency 2 punch_cells 3-17 USER_GRAPH 1 -detach USER_GRAPH 5 Example 21 -chart_title "Hto Concentration Profile: Filter1 | Clay | Filter2" -axis_scale x_axis 0 2.2220e+001 -axis_scale y_axis 0 1.2e-9 -axis_scale sy_axis 0 1.2e-9 -axis_titles "Distance, in millimeters" "Free pore-water molality" "Total molality" -headings Hto_free Hto_tot -plot_concentration_vs x -initial_solutions true 10 IF cell_no = 3 THEN xval = 0 ELSE xval = get(14) 20 IF (1 = 0 OR cell_no > 4) THEN GOTO 60 30 IF (cell_no = 4) THEN xval = xval + 0.5 * 1.8000e-003 40 IF (cell_no > 4 AND cell_no < 5) THEN xval = xval + 1.8000e-003 50 GOTO 200 60 IF (cell_no = 5) THEN xval = xval + 0.5 * 1.8000e-003 + 0.5 * 1.7109e-003 70 IF (cell_no > 5 AND cell_no < 16) THEN xval = xval + 1.7109e-003 ELSE GOTO 90 80 GOTO 200 90 IF (cell_no = 16) THEN xval = xval + 0.5 * 1.7109e-003 + 0.5 * 1.6000e-003 100 IF (cell_no > 16 AND cell_no <= 16) THEN xval = xval + 1.6000e-003 110 IF (cell_no = 17) THEN xval = xval + 0.5 * 1.6000e-003 200 y1 = TOT("Hto") 210 plot_xy xval * 1e3, y1, color = Blue, symbol = Plus 220 IF cell_no = 3 THEN put(y1, 15) 230 IF (cell_no < 5 OR cell_no > 15) THEN GOTO 400 240 y2 = SYS("Hto") / (tot("water") + edl("water")) 250 REM y2 = y2 / 1.4281e+001 260 plot_xy xval * 1e3, y2, symbol = Circle, y_axis = 2 270 IF (cell_no > 6) THEN GOTO 400 280 IF 1 THEN plot_xy 1.8000e+000, get(15), color = Black, symbol = None 290 IF 1 THEN plot_xy 2.0620e+001, get(15), color = Black, symbol = None 300 put(0, 15) 400 put(xval, 14) END END -------------------------------- End of Run after 22.075 Seconds. --------------------------------