mirror of
https://git.gfz-potsdam.de/naaice/iphreeqc.git
synced 2025-12-15 16:18:22 +01:00
3d759f6997
1777f2bc Merge branch 'mix_error' of github.com:dlparkhurst/phreeqc3-1 into mix_error da8a610a Fixed examples to have fewer warnings. Fixed headers of some database files. Changed some output from warnings to comments. git-subtree-dir: phreeqc3-examples git-subtree-split: 1777f2bc99bd9b655e838908e0fdfb81bfdd1e3f
363 lines
18 KiB
Plaintext
363 lines
18 KiB
Plaintext
#DATABASE ../database/phreeqc.dat
|
|
TITLE Diffusion through Opalinus Clay in a radial diffusion cell, Appelo, Van Loon and Wersin, 2010, GCA 74, 1201
|
|
# NEW: viscosity effects in solution and Donnan EDL, (and, possibly correct co-ion in Donnan layer to the DLVO values)
|
|
|
|
KNOBS; -tol 1e-16; -diagonal_scale true
|
|
|
|
SOLUTION_MASTER_SPECIES
|
|
# element species alk gfw_formula element_gfw
|
|
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
|
|
# start with finding tortuosity from HTO
|
|
Hto = Hto; log_k 0; -gamma 1e5 0; -dw 2.3e-9 0 0 0 0 0 0.5 # diffusion coefficient is multiplied by (viscos_0 /viscos)^0.5, the viscosity of the DDL is calculated.
|
|
# estimate f_free and f_DL_charge, increase tortuosity
|
|
Cl_tr- = Cl_tr-; log_k 0; -gamma 3.5 0.015; -dw 1.35e-9 0 0 0 0 0 0.5 # increase tortuosity for anions: 2.03e-9 / 1.35e-9 = 1.5
|
|
# use erm_ddl to fit Na
|
|
Na_tr+ = Na_tr+; log_k 0; -gamma 4.0 0.075; -dw 1.33e-9 0 0 0 0 0 0.5 ; -erm_ddl 1.3
|
|
# use interlayer diffusion to fit Cs
|
|
Cs+ = Cs+; log_k 0; -gamma 3.5 0.015; -dw 2.07e-9 0 0 0 0 0 0.5 ; -erm_ddl 1.3
|
|
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
|
|
|
|
|
|
|
|
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
|
|
# uncomment tracer concentrations and kg water 1 by 1... (the experimental water volumes are different)
|
|
Hto 1.14e-6; -water 0.2
|
|
# Cl_tr 2.505e-2; -water 0.502
|
|
# Cs 1; Na_tr 1.87e-7; -water 1.02
|
|
SELECTED_OUTPUT
|
|
-file radial; -reset false
|
|
USER_PUNCH
|
|
# Define symbols and pi...
|
|
1 nl$ = EOL$ # newline
|
|
2 x$ = CHR$(35) # cross '#'
|
|
3 sc$ = CHR$(59) # semicolon ';'
|
|
4 pi = 2 * ARCTAN(1e10) # 3.14159...
|
|
|
|
# Define experimental parameters...
|
|
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)
|
|
# adapted for the harmonic mean calc's in version 3.4.2
|
|
140 CEC = 0.12 * rho_b_eps # CEC / (eq/L porewater)
|
|
150 A_por = 37e3 * rho_b_eps # pore surface area / (m²/L porewater)
|
|
151 correct_$ = ' false'
|
|
# 152 correct_$ = ' true' # if 'true' correct the co-ion concentrations in the Donnan volume
|
|
|
|
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)
|
|
# experimental times (seconds) for HTO, 36Cl, 22Na and Cs, respectively,
|
|
# in order of increasing times...
|
|
200 DATA 86400 * 20, 86400 * 40, 86400 * 45, 86400 * 1000
|
|
210 READ exp_time(1), exp_time(2), exp_time(3), exp_time(4)
|
|
# scale y1-axis (flux) (not used)...
|
|
230 DATA '1', '1', '1', '1'
|
|
240 READ scale_y1$(1), scale_y1$(2), scale_y1$(3), scale_y1$(4)
|
|
# scale y2-axis (mass) (not used)...
|
|
260 DATA '1', '1', '1', '1'
|
|
270 READ scale_y2$(1), scale_y2$(2), scale_y2$(3), scale_y2$(4)
|
|
# scale max of the profile y axes...
|
|
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)
|
|
|
|
# Define model parameters...
|
|
350 Dw = 2.5e-9 # default tracer diffusion coefficient / (m²/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.11 # fraction of free pore water (0.01 - 1)
|
|
400 f_DL_charge = 0.48 # fraction of CEC charge in electrical double layer
|
|
# 400 f_free = 0.2 : f_DL_charge = 0.5 # higher f_free ===> higher f_DL_charge, found from Cl- and Na+
|
|
410 tort_n = -1.00 # exponent in Archie's law, found from HTO
|
|
420 G_clay = por_clay^tort_n # geometrical factor
|
|
430 interlayer_D$ = 'true' # 'true' or 'false' for interlayer diffusion
|
|
440 G_IL = 1300 # geometrical factor for clay interlayers... the initial rise of Cs suggests stagnant water, see Appelo et al for the calculation
|
|
450 punch_time = 60 * 60 * 6 # punch time / seconds
|
|
460 profile$ = 'false' # 'true' or 'false' for c/x profile visualization
|
|
470 IF nfilt1 = 0 THEN thickn_filter1 = 0
|
|
480 IF nfilt2 = 0 THEN thickn_filter2 = 0
|
|
|
|
# See which tracer is present...
|
|
490 IF tot("Hto") > 1e-10 THEN tracer = 1 ELSE \
|
|
IF tot("Cl_tr") > 1e-10 THEN tracer = 2 ELSE tracer = 3
|
|
|
|
# Define clay pore water composition...
|
|
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'
|
|
|
|
# Define phases in which the tracers precipitate...
|
|
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
|
|
|
|
# Write solutions for the cells...
|
|
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)) + ' viscosity calc' + ' correct ' + correct_$
|
|
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'
|
|
|
|
# Write phases in which the tracers precipitate...
|
|
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'
|
|
|
|
# Define mixing factors for the diffusive flux between cells 1 and 2:
|
|
# J_12 = -2 * Dw / (x_1 / g_1 + x_2 / g_2) * (c_2 - c_1)
|
|
# Multiply with dt * A / (V = 1e-3 m³). (Actual volumes are given with SOLUTION; -water)
|
|
# Use harmonic mean: g_1 = por_1 / G_1, g_2 = por_2 / G_2, x_1 = Delta(x_1), etc.
|
|
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
|
|
# Find max time step = 0.5 * V_water * dx * G_factor / (Dw * por * A * fbc)
|
|
# V_water = por * pi * height * ((r + dr)^2 - r^2)
|
|
# A = por * pi * height * r * 2
|
|
# At the inlet of the tracers, fbc = 2...
|
|
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
|
|
# Perhaps the clay has very small cells...
|
|
1540 IF nfilt1 = 0 THEN ff = ff1 ELSE IF ff1 * 2 < ff THEN ff = ff1 * 2
|
|
# Or at the filter1-clay transition, fbc = 1...
|
|
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
|
|
# Perhaps filter2 has very small cells...
|
|
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
|
|
# Check with punch times, set shifts...
|
|
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)
|
|
# Write mixing factors...
|
|
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'
|
|
|
|
# Write TRANSPORT...
|
|
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))
|
|
|
|
# Write USER_GRAPH...
|
|
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 MOL/M2/S" "ACCUMULATED MASS, IN MOL"'
|
|
|
|
|
|
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
|
|
|
|
# Continue with Cs...
|
|
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...
|
|
|
|
# Write TRANSPORT and USER_GRAPH for concentration profile...
|
|
3000 punch nl$ + 'TRANSPORT'
|
|
3010 punch nl$ + ' -shifts 0'
|
|
3020 punch nl$ + ' -punch_fr 2' + sc$ + ' -punch_c 3-' + TRIM(STR$(2 + stag))
|
|
# Write USER_GRAPH...
|
|
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"))'
|
|
# Remove REM if total conc's per kg solid must be plotted (and adapt axis_titles)...
|
|
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
|
|
PRINT
|
|
-selected_out false
|
|
-status false
|
|
INCLUDE$ radial
|
|
END
|