#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