FARF31: radionuclide migration along a stream tube

This page calculates what SKB’s far-field code FARF31 calculates: the release of radionuclides from a stream tube in fractured rock, given the release into it. Water carries the nuclides along the tube with dispersion; they diffuse into the stagnant pore water of the rock matrix beside the flowing water, to a finite depth, and sorb there; and they decay, parents into daughters. The model is that of Norman and Kjellbert (SKB TR 90-01), with the element-specific diffusivities added for SR 97, as summarised in SKB R-04-51.

The page implements the model as published in TR 90-01 and R-04-51, with numerical methods of its own (see how the page differs), and adds sorption on the fracture surfaces as an option. Everything runs in your browser: the files you open or drop are read there and are not sent anywhere.

The model

Assumptions (TR 90-01, chapters 2–3): stationary flow; transversal dispersion neglected; the longitudinal dispersion and the exchange between flowing and stagnant water replaced by averages over the tube’s cross-section, and constant along the tube; diffusion into the matrix one-dimensional, perpendicular to the flow, to a depth x₀; linear, reversible, instantaneous sorption in the matrix, and, as an option of the page, on the fracture surfaces (FARF31 has none there); decay chains without branches. The distance along the tube is measured as accumulated water travel time ζ, 0 ≤ ζ ≤ t_w, so the tube is described by t_w and a Peclet number rather than by lengths and velocities.

In the flowing water

R_f,i ∂c_i/∂t = −∂c_i/∂ζ + (t_w/Pe) ∂²c_i/∂ζ² + a_w D_e,i ∂c_p,i/∂x|_x=0 − λ_iR_f,ic_i + λ_i−1R_f,i−1c_i−1,   R_f,i = 1 + K_a,ia_w

K_a is the surface sorption coefficient of the element on the fracture walls: the amount sorbed per m² of wall over the concentration in the water (m³/m², that is m). The sorbed nuclides decay with the dissolved ones, and their daughters take up their own equilibrium. With K_a = 0, R_f = 1, this is the equation of TR 90-01 and R-04-51.

In the rock matrix, 0 ≤ x ≤ x₀

R_i ∂c_p,i/∂t = D_e,i ∂²c_p,i/∂x² − R_iλ_ic_p,i + R_i−1λ_i−1c_p,i−1,   R_i = ε + K_d,iρ

Conditions

• Nothing in the tube or the rock at t = 0.
• c_p = c at the fracture wall (x = 0) and ∂c_p/∂x = 0 at x = x₀.
• The inlet flux is given: F_in,i(t) = Q(c_i − (t_w/Pe) ∂c_i/∂ζ) at ζ = 0, and c_i → 0 far downstream.
• The output is the same flux at ζ = t_w. The flow Q cancels, and is not an input.

The flow-wetted surface a_w is per volume of flowing water, so t_wa_w is the F-factor (a/m), the quantity that governs the retention (R-04-51, section 2.3). ρ is 2700 kg/m³ in FARF31 unless its settings file gives another (RHOP).

The solution in the Laplace domain

Transformed in time (variable s), a single nuclide gives an ordinary differential equation in ζ whose decaying solution makes the transfer function from inlet flux to outlet flux

H(s) = exp[(Pe/2)(1 − √(1 + 4t_wg/Pe))],   g = R_f(s + λ) + a_w √(D_eR(s + λ)) tanh(x₀√(R(s + λ)/D_e)),

the transform of the unit response: the outlet flux after a unit pulse at the inlet. For a chain, the matrix equation couples each nuclide to its parent, and so to every ancestor: a parent that decays in the matrix produces a daughter there, which diffuses back. Written as vectors, with Λ the decay matrix (λ_i on the diagonal, −λ_i−1 below it) and R, R_f, D_e diagonal,

M = D_e⁻¹(sI + Λ)R,   G = (sI + Λ)R_f + a_wD_e τ(M),   T(s) = H(G),   τ(z) = √z tanh(x₀√z),

where τ(M) and H(G) are functions of lower-triangular matrices. T_ij is the transform of the response of nuclide i to a unit pulse of nuclide j. TR 90-01 (Appendix C, Proposition 3) writes the same solution as sums over exponentials with recursively defined coefficients; its factors 1/(F_i − F_k) and 1/(h_i² − h_l²) are divided differences here. The page evaluates the matrix functions by Parlett’s recurrence when the diagonal is well separated, and otherwise by the sum over paths, with the divided differences of close points (the same element with close decay constants, U-238 and U-234 for one) taken from Taylor series about their centre, so that equal half-lives need no special treatment.

The singularities of T lie on the negative real axis: branch points at −λ for an infinitely deep matrix, poles of the tanh for a finite one, and the branch point of the square root in H (TR 90-01, Appendix C.2). An infinitely deep matrix and plug flow (Pe = ∞, when every nuclide exchanges with the matrix and the members of a chain share one R_f) are allowed.

Numerical inversion

The saddle-point contour (default)

For each response and time t the page finds the saddle point s* of e^stT_ij(s) on the real axis: the point where −d ln T/ds, the mean time of the exponentially tilted response, equals t. It integrates along the parabola through it, s = s* + iy − κy², κ = ψ″(s*)/(2t), which is the path of steepest descent to second order — for advection and dispersion alone exactly, the integrand then being a Gaussian in y — with the trapezoidal rule, halving the step until two sums agree. The step also stays inside the strip of analyticity set by the nearest singularity s₀. The terms are then of the size of the answer, so the result is accurate relative to itself on the leading edge, where the transform behaves like a delay and fixed contours amplify rounding by many orders of magnitude, and at high Peclet numbers. Against closed forms it gives ten or more significant digits.

Three safeguards keep it so for responses with a sharp front and a long slow tail, and at high Peclet numbers. In the tail the saddle comes close to s₀, a branch point; the parabola's vertex then stays 1.5/t to the right of it and its focus on it, so that the path wraps the branch point as Hankel's contour does, at the scale on which e^st changes. At the branch point of the fracture's own transform, where 1 + 4t_wg/Pe = 0, the transform reaches e^Pe/2, and it stays large near the real axis beyond it; the parabola is flattened until no member's transform grows along it, and flattened further should its terms, once fallen away, grow again; the step follows the fastest turning phase of any member's transform along the path. The real axis is sampled finely enough that the saddle of a broad response is placed within its own width, also when the singularity lies far out. Under plug flow, where a finite matrix puts a ridge beside every pole of tanh, two halvings in a row must agree. Where the path finds no clear way or its sums do not settle within their budget, the sample comes from de Hoog's method instead, when that agrees with itself at twice its terms.

Talbot’s fixed contour

The page’s alternative: s = rθ(cot θ + i), r = 2M/(5t) with M = 28 nodes (Abate and Valkó 2004), scaled out to the saddle point before a front. Six or seven digits in ordinary cases, fewer above a Peclet number of about 100, where the saddle-point contour is the one to use.

De Hoog, Knight and Stokes

The trapezoidal rule on the Bromwich line Re s = γ, accelerated by the quotient-difference continued fraction, with period 2 × 2t (de Hoog et al. 1982). It shares nothing with the contour methods and serves as the check: with the check ticked, a spread of samples of every response is recomputed by it (or by the saddle-point method when de Hoog is the chosen one) and the Summary gives the largest relative difference over values above 10⁻⁶ of each response’s peak.

Responses, convolution and the output

Each response is sampled with its first and second derivatives (the transforms sT and s²T cost nothing extra), on a logarithmic grid refined until the quintic Hermite interpolant predicts every new midpoint to 2·10⁻⁸ (with the default inversion; to what the alternatives deliver otherwise). Its range comes from the saddle-point estimate of ln h(t). The integral of the interpolant is the mass balance: it must equal T_ij(0), the fraction that leaves the tube, or the part of it that has left by the response's last time when the output ends before its tail does (the integral of h to that time, by de Hoog's method). The page checks every response against it; one that misses is computed again from earlier times on a finer grid, and if it still misses, the status line and the Summary say so and name it.

A release series is piecewise linear between its points and zero before the first and after the last; a time given twice is a step. The convolution of such a series with the interpolated response is done exactly: between the merged breakpoints of the two, the integrand is a polynomial of degree six, integrated by four-point Gauss–Legendre. The daughters’ own releases add to what their parents produce.

Unless given, the output runs from where the first response reaches 10⁻¹⁰ of its peak after its release starts to where the last has fallen below 10⁻⁹ of its peak after its release ends, as far as 10¹² a; when that limit cuts a release short, the Summary says so. The output grid is logarithmic (20 points per decade by default), with the input’s breakpoints and those moved by each response’s peak time, refined until linear interpolation between neighbours is within RELINT of every midpoint. Each peak is refined by a golden-section search. The activity is λN_A per mole with FARF31’s constants, N_A = 6.022045·10²³ and a year of 365.2422 days.

A single nuclide or a short chain takes a fraction of a second. The cost grows with the number of responses (a chain of seven with a release of every member has 28) and with their length, since each is inverted at a few hundred times; such a case takes a few seconds, off the page’s main thread.

What the inputs mean

XX in KDR_XX, KA_XX and DE_XX is the first two characters of the nuclide name, upper or lower case; for a one-letter element the first digit is part of it (U238 → U2, I129 → I1). Nuclides that share those characters share their K_d, K_a and D_e, and the Input tab keeps them together.

Files

Drop files anywhere on the page or use Open…; several at once is fine, and in.dat is read before in.par and in.ts. # starts a comment in all of them.

• in.dat: PRINT ON|OFF|DEBUG, CASENAME name, DIFFUSIVITY SINGLE|ELEMENT_SPECIFIC, then one line per nuclide, NAME THALF IDAUGH ISOURC: half-life in years, IDAUGH = 1 when the next line is its daughter, ISOURC = 1 when in.ts has a series for it.
• in.par: KEYWORD value in any order: TW, PECLET, ASPEC, EPS, DE (with SINGLE), PENDEP, and per element KDR_XX and (with ELEMENT_SPECIFIC) DE_XX; F, the F-factor t_wa_w, may take the place of TW or ASPEC. The page writes the element lines in the order the elements first appear in in.dat. KA_XX, the fracture’s surface sorption coefficient, is the page’s own keyword: an in.par with KA_ lines goes beyond FARF31’s format, so the page writes them only when some K_a is not zero, and reads K_a = 0 where they are absent.
• in.ts: for each source, in in.dat’s order, its name on a line of its own and then time rate lines. A CSV or text table with a header time, name1, name2, … is read as well, and a two-column table goes to the nuclide selected on the Input tab.
• casename31.prm: numerical settings. FARF31 picks its inversion routine from a line holding only BROMEX (its default), TALBOT or STEAMR, reads KEYWORD value lines such as NPMIN, NPMAX (at most 128 output points), RELINT, BQMIN and RHOP (the rock density), and skips lines it does not know. The page reads those — TALBOT as its fixed Talbot contour, BROMEX as de Hoog’s method, STEAMR as its default — and its own METHOD, TSTART and TEND lines, and lists anything else as not used. It writes BROMEX for FARF31, the page’s method on a METHOD line, and NPMIN and NPMAX only within FARF31’s limits.
• out.ts is written in FARF31’s layout: per nuclide, time mol/a Bq/a, each nuclide from where its rate first reaches BQMIN Bq/a. Dropping an out.ts (FARF31’s or the page’s) draws it over the release as markers and lists the differences at its times.
• out.response, FARF31’s PRINT DEBUG output: the unit responses per pair, each followed by its integral.
• A case (.json) holds everything on the page; CSV holds the release table or the responses.

How the page differs from the original program

• Inversion. FARF31 offers three routines (R-04-51, section 3.2): one based on Talbot’s contour (Talbot 1979, as described in TR 90-01 section 4.3, with constant parameters and a shift from the rightmost singularity), the Bromwich-line routine BROMEX (Gustafson 1991), its default, and Dahlquist’s “steamroller” (Dahlquist 1993). The page uses the saddle-point contour by default, a fixed Talbot contour and de Hoog’s method; see above.
• Convolution. FARF31 convolves the two time series with the trapezoidal rule on the union of their points (TR 90-01, section 4.4). The page integrates the interpolated response against the piecewise-linear release exactly.
• Time grids. FARF31 lays out its grids with the PROPER Time Series Manager (TR 90-01, sections 4.3 and 4.4: first upper bounds from the first two moments of the response and Chebyshev’s inequality, then a controlled interpolation error). The page refines its own grids adaptively (see above), so its times are not FARF31’s; compare at the same times, as dropping an out.ts does.
• Chains. The same model, written as matrix functions rather than as TR 90-01’s recursion; the two are the same solution (see the Laplace domain).
• Additions. Sorption on the fracture surfaces (K_a), the F-factor as an input, an infinitely deep matrix, plug flow, releases built from shapes, a second inversion method as a check, and the mass balance of every response.

How it was checked

The model’s tests (resources/tests/farf31) compare it with closed forms — the inverse Gaussian of advection and dispersion without matrix interaction, the classical solution for matrix diffusion into an infinite matrix under plug flow, the Bateman solution for chains whose members move alike, each also with sorption on the fracture surfaces, the transmission T(0) and the mass balance — and with a separate high-precision implementation of TR 90-01’s recursion (mpmath, 40 digits and more; with the fracture term added for a case with K_a) on made-up cases, among them nine tubes with a sharp front and a long slow tail (Pe 30 to 300, a matrix that holds 0.1 to 10 times the water and fills over 250 travel times). The page agrees with the latter to about 10⁻¹⁰ (10⁻⁸ at worst) over values above 10⁻⁶ of each peak, and every response's integral with T(0). Sharper fronts (Pe 1000 to 10 000), where the fracture's transform reaches e^Pe/2, are checked against de Hoog's method. It also reproduces the outputs of the original program for made-up test cases.

References