Reactive transport in the browser

Write a reaction set the way you would write it down, say whether it sits in a stirred batch or in a column with diffusion and advection, and this page builds the equations, differentiates them into an analytic sparse Jacobian and integrates them with the same stiff solvers as the canister radiolysis page. Nothing is uploaded: the whole thing runs here.

The model text

Six sections, of which two are needed: <SPECIES> and <REACTIONS>. Everything after a # is a comment. A rock matrix beside the column is a few more settings; see A rock matrix.

Colour the syntax, on the Model tab, paints comments, section headings, numbers, function names and the words that mean something where they stand — a setting and the words its value may take, fixed and D on a species line, k and r in a reaction, all and matrix in a parameter. It is a second copy of the text drawn underneath the box you are typing in, which is the only way a plain text area can be coloured while keeping its own undo, selection and caret; turn it off for a very large model, or if the colours are in the way. The choice is remembered between visits.

<SETTINGS>

<SPECIES>

H2O    55.56   fixed
H+     1.0E-7  D=9.3E-9  left=1.0E-7
UO2                      # no number: starts at zero, does not move

A name, the concentration it starts at, and then any of: fixed to hold it there for ever, D= its diffusion coefficient in m²/s, and left= / right= the concentration at a Dirichlet face. A species with no D does not move, which is how a surface site or a mineral is written.

fixed is a statement about the equation, not a nudge: that species' row in the right-hand side and in the Jacobian is set to exactly zero, so nothing — reaction, diffusion, advection, a Dirichlet face, an equilibrium — can move it. It still drives everything else, and is still consumed on paper by any reaction that consumes it; it simply is not replenished from anywhere the model has to account for. That is what a large water reservoir, a buffered pH or an imposed dose field is. Because the row is zero, the answer for such a species is the number you wrote, and that number is what the charts and the CSV show: a stiff step leaves round-off in it — about one part in 10¹³ on the built-in model over nine hundred steps — and that is arithmetic rather than chemistry, so it is written back out. If it is ever larger than round-off could explain, the run says so instead of hiding it.

A fixed species inside an <EQUILIBRIUM> is held there too, and the bookkeeping follows: fix the H+ of a water equilibrium and charge stops being conserved, because the buffer supplies whatever the relation calls for. What is left is one relation in one unknown, and OH- goes to K/[H⁺]. Fixing every species of a relation is refused, since then there is nothing left to satisfy it with.

A name may be a chemical formula: e-, H+, OH-, C2O4-2, UO2+2, U_site.

R= (or mass=, the same thing) is the coefficient on the left of the equation rather than the right: R dC/dt = …. A retardation factor for a sorbing species, a porosity, a water content — anything that divides everything happening to that species, transport and reactions alike. R=4 turns diffusion at D into diffusion at D/4 and a decay at k into one at k/4. It may be a plain number or the name of a parameter, in which case it varies from cell to cell. It has to be greater than zero: a zero would make that equation algebraic rather than differential, which is what <EQUILIBRIUM> is for.

<REACTIONS>

A + 2 B => C, k = 1.4E11                       # mass action: k[A][B]²
A <=> B, kf = 1.0, kb = 0.5                     # both ways
S => P, r = vmax*[S]/(km + [S]), vmax = 2, km = 0.1   # a rate law of your own
=> H2O2, k = 1.0E-8                            # nothing on the left is a source
OH- => , k = 0.1                               # nothing on the right is a sink

The equation comes first, then the rate, then any constants the rate uses. With k = you get mass action and the rate is that constant times each reactant raised to its coefficient. With r = you write the rate yourself and a concentration goes in brackets — [S] — so that a name like H+ is never read as arithmetic. A power needs a plain number as its exponent.

The time unit

A model has one time unit and everything in it is in that unit: the rate constants, the diffusivities, TEND, a velocity. TIME_UNIT says which, and it converts nothing — the solver integrates the numbers it is given. What the setting does is let the page say so: label the time axis, read 500 a in the Times box, print a travel time, head a CSV column. Without it everything is called seconds, which is what a model that says nothing means.

So a geosphere model written in years has half-lives in years, D_e in m²/a and a velocity in m/a, and writes TIME_UNIT = year; the radiolysis models are in seconds throughout and say nothing. Mixing the two inside one model is the one mistake this cannot catch, and never could: a diffusivity in m²/s beside a half-life in years is simply a different model from the one that was meant.

The two ends of the column

Four kinds of face, under the names the transport literature gives them. Which you want depends on what is on the other side: a reservoir, a wall, water arriving, or nothing at all.

The difference between the first and third types is worth knowing, because it is a common way to get an inlet wrong. A Dirichlet face adds a diffusive flux D(c_b − C)/d on top of what the flow brings, so more enters than the water carries, and the column's mass balance does not close: measured on the test column in resources/tests/rtm/, 6.7 % more goes in than u(c_in − c_out) accounts for. With robin and free the same column balances to one part in 10¹⁵, and its profile matches the flux-type closed solution of van Genuchten and Alves, converging on it by a factor of four each time the cells are doubled.

As arithmetic, free and neumann are the same face — the advective outflow applies at every cell whatever the far boundary says, so a Neumann outlet already lets the flow leave. They are separate words because they say different things and are wrong in different ways, and the page can only tell you which is wrong if it knows which you meant.

A rock matrix: dual porosity

Water flows along a fracture; the rock beside it is stagnant, porous, and open to diffusion. A species carried along the fracture leaks sideways into the rock's pore water, sorbs there, and comes back out long after the water that brought it has gone — which is what holds a radionuclide up in a geosphere by orders of magnitude more than the water's own travel time. Set MATRIX_CELLS above zero and every cell of the column gets a chain of that many stagnant cells behind it, reaching MATRIX_DEPTH metres into the rock:

<SETTINGS>
MODE = transport
CELLS = 20                # along the fracture
LENGTH = 235              # m; with VELOCITY = 1 the length is the travel time in the time unit
ADVECTION = 1
VELOCITY = 1
DIFFUSION = 1
PECLET = 10
LEFT = robin              # water arriving at the stated concentration
RIGHT = free
MATRIX_CELLS = 20         # layers of rock behind each cell
MATRIX_DEPTH = 4.5        # m into the rock
MATRIX_FIRST = 1E-4       # m, the layer at the wall; 0 has it worked out
MATRIX_POROSITY = 0.0019
WETTED_SURFACE = 340      # m² of rock wall per m³ of flowing water; = F / t_w

<SPECIES>
Cs  0  left=1  R=1  Dm=8.5E-7  Rm=2.7    # Dm = De in the rock; Rm = eps_m + rho*Kd

Every species has its Dm=, the effective diffusivity in the rock in m²/s, and its Rm=, the rock's capacity for it: pore water plus sorbed, ε_m + ρ_mK_d, the quantity the SKB reports call the capacity factor. With no Rm the species gets the porosity alone, which is a species that does not sorb; with no Dm it never enters the rock. Like R, Rm may name a parameter and so vary from cell to cell.

The layers grow geometrically with depth from MATRIX_FIRST, because the gradient is steepest at the wall and all but flat far in. Left at zero the first layer is worked out so that the ratio between layers is e, which is the FARFCOMP rule, reproduced step for step so that the thicknesses come out bit for bit the same as that code's. The far face of the last layer is closed: the depth is a depth, not a boundary. WETTED_SURFACE is the area of fracture wall a cubic metre of flowing water touches; APERTURE (2b) says the same thing as 2/aperture. Only one of the two is needed.

Concentrations are per unit pore water everywhere, and so are the rate laws. A cubic metre of rock holds only ε_m of water, so a rate r per unit water changes the rock's inventory by ε_mr, and it is the inventory, R_mC, that the rock's equation is written for: R_m dC/dt = ε_mr + diffusion. That is right for chemistry in the water and wrong for one thing: radioactive decay, which does not ask where an atom sits, sorbed or dissolved. Write such a rate with on = inventory and it removes the same share of everything in the cell:

Cs => , k = 7.3E-10, on = inventory     # decay: of everything, sorbed included
Cs + X => Y, k = 1.2E-4                   # chemistry: in the water, on = water is the default

Told it acts on the water alone, decay in a sorbing rock runs at ε_m/R_m of its real rate; on the reference case below that is the difference between a release fraction of 0.116 and 0.717.

A line in <INITIAL> or <PARAMETERS> names cells of the column. Add matrix after it to mean the rock behind those cells instead, or fracture to mean the flowing water alone. A parameter with neither applies to both — right for a temperature, wrong for a source into the water, which would otherwise also feed every layer of rock behind that cell. In a parameter's expression j is the layer (0 in the fracture) and xm the depth to its centre in metres.

From the SKB parameterisation. A flow path there is a travel time t_w, a transport resistance F, a Péclet number, a penetration depth, and per species D_e, K_d and the rock's porosity and density. Here that is LENGTH = t_w with VELOCITY = 1 (so time is the model's unit, years if the rates are per year), WETTED_SURFACE = F/t_w, PECLET, MATRIX_DEPTH, and per species Dm = D_e, Rm = ε_m + ρ_mK_d,m, and R = 1 + K_d,f·a_w for sorption on the fracture coating. A unit release into the first cell is a source => X, k = S with S 0 1/Δx fracture in <PARAMETERS>; the release out is v times the last cell's concentration.

What was checked. The layer thicknesses and every rate — advection, dispersion, wall exchange, layer to layer — against SKB's own implementation of TR-19-06's formulation, in five cases including twenty layers from 4×10⁻⁸ m, to one part in 10¹⁴; two of its release curves, one with sorption and decay, to 2×10^-9 and 2×10^-10 at a tolerance of 10^-9; Neretnieks' closed solution for a fracture beside a deep matrix, converging at the upwind scheme's first order as the fracture is refined; and a closed box that conserves mass to 10^-14 and shares it between water and rock exactly by capacity.

On the charts, a Layer picker beside the cell chooses the fracture or a rock layer for the time chart and the gradient, and Against distance can be drawn across the fracture at one layer or into the rock at one cell, the fracture at depth zero.

What may go in an expression

The same arithmetic in a rate law, a constant and a parameter: + - * /, brackets, and ^ or ** for a power. Numbers may be written 1.4E11 or 1.4e11. The functions are:

Three of those are barred from a rate law, and the reason is the Jacobian rather than taste. min and max have a derivative that jumps where their two arguments cross, so an analytic Jacobian would be wrong on one side of the crossing; pow(a,b) differentiates to something containing log a, which is undefined at a zero concentration — which is exactly where a rate law gets evaluated. A power with a plain number as its exponent, [A]^2, has neither problem and is what to write. <PARAMETERS> is worked out once when the model compiles and never differentiated, so all three may be used there.

step and ramp are allowed, with the derivative of the kink taken as zero: they are flat everywhere they are differentiable, and a stiff solver steps over the corner rather than through it. A rate that switches hard is still worth avoiding where a smooth one will do.

Which names are in scope depends on where the expression is. A rate law may use the concentrations it writes in brackets, the constants named on its own line, and any parameter; a constant on the reaction line wins if a parameter has the same name. A parameter may use x, the cell centre in metres, i, the cell index, and any parameter named above it.

<INITIAL>

X   0-9   1.0      # the first ten cells
O2  all   2.1E-4
Fe  20    0.5      # one cell

Optional, and only useful with transport: it overrides the concentration in <SPECIES> cell by cell, which is how you start a front, a plume or a layer.

<PARAMETERS>

DOSE  all  0.64                  # the same everywhere
DOSE  all  0.64*exp(-x/3.0E-5)   # x is the cell centre, in metres
POR   0-4  0.15                  # and a later line overrides an earlier one
POR   5    0.3

A named number that any rate law may use and that may differ from cell to cell. The value is an expression in x, the distance to the cell centre in metres, i, the cell's index, and any parameter named above it. It is worked out once, when the model is compiled, so a parameter is constant in time — it varies in space, not with the solution.

This is how a dose rate that dies away from a surface is written, or a porosity that changes down a column, or a rate constant that follows a temperature profile. A rate law may use a parameter by name, exactly as it uses a constant declared on its own line; a constant on the reaction line wins if both have the same name.

<EQUILIBRIUM>

      <=> H+ + OH-,    logK = -14
CO2   <=> H+ + HCO3-,  logK = -6.35
HCO3- <=> H+ + CO3-2,  logK = -10.33
NpO2+ <=> NpO2_s,      K = 2.0E4, kf = 1.0E3   # also held during the run

A relation that has to hold rather than a rate that has to be integrated: ∏[C]ν = K, with logK its base-ten logarithm, which is the form a database gives. An empty left-hand side means the reaction has no reactant, as it does for water dissociating.

With EQUILIBRATE = 1 the starting concentrations are put on these relations before the run: a damped Newton iteration in ln C, run separately in every cell, over the species the relations mention and no others — anything they do not touch is left exactly where it was put, zero included. What it may not change is worked out from the relations themselves — the left null space of their stoichiometry, which gives the conserved totals without the page needing a table of elements. So writing H+ at 10^-3 and OH- at 10^-3 gives 10^-7 of each, and the difference between them, which is a charge, is exactly what it was.

Give kf as well and the equilibrium is held during the run too, as a fast reversible pair with kb = kf/K. Without it, it applies to the initial state and nothing else — these are solvers for y' = f(y), and an equilibrium held exactly throughout is an algebraic constraint, a different kind of problem. A kf a few orders above the fastest rate in the model holds the relation to within the solver's tolerance; making it enormous only costs steps.

What it actually solves

For every species s in every cell i:

dC/dt = (what the reactions make and consume in cell i)
      + (diffusion across the two faces of cell i)
      + (advection in from upwind)

Diffusion is a finite-volume difference: the flux between two cells over the distance between their centres, divided by the width of the cell it lands in. On an even grid that is the textbook D/Δx²; on an uneven one it stays conservative, which the textbook form does not. Measured on a closed box, mass is held to about one part in 10¹⁵, and against the analytic erfc profile the error falls by four each time the cells are doubled — second order, as it should be.

Reading the charts

The time chart draws the species ticked in the left-hand panel, at one cell. Against distance draws one species across all the cells, at several times. Which times is up to you: the Times field takes a list — 0, 1 h, 30 d, 500 y, or bare numbers for seconds — and leaving it empty spreads eight over the run.

Gradient puts the whole run in one picture: time across, distance down, concentration as colour. Where Against distance shows a handful of moments, this shows every stored one, so a front's arrival and how far it gets are a shape rather than a set of curves to compare. Distance runs up the axis from zero at the origin, so the left-hand face — the surface the chemistry happens at — is along the bottom.

The bands across it are the cells: with twenty of them the picture has twenty rows and says so. Smooth interpolates between them, which looks better and claims a resolution the run does not have. log concentration plots log₁₀, leaving anything at or below zero blank rather than dropping it to the bottom of the scale, and the note says how many that was. The from time matters more here than elsewhere: a stiff run's first steps are femtoseconds, and without it fifteen decades of an empty axis push the interesting part into the last inch.

A time you ask for is not usually a time the solver stopped at, and nothing is invented to fill the gap: the nearest profile that was actually computed is drawn, the legend names the time it belongs to rather than the time you asked for, and the note under the toolbar says which times moved. Late in a stiff run the steps are wide — tens of thousands of seconds on the built-in model — so ask for 10 hours there and you may well get 10.8.

Taking the run away with you

CSV writes the species ticked in the panel, at the cell the chart is showing, a row per stored step. View in HDF5 Browser writes the whole run as an HDF5 file and opens it in the HDF5 Browser in another tab, without the file ever reaching the disk; HDF5 is the same file as a download. Both need a run behind them, and neither sends anything anywhere: the file is built in this browser and handed to a page of the same origin.

Inside it, /time is the clock in the model's own unit — it converts nothing, any more than the page does — and /Species holds every species where the chart was looking: the batch, or the cell and rock layer the panel had chosen. A column run also has /Cells/cell_007/… for every cell and /Cells/cell_007_rock_03/… for the rock behind it, so a species searched for by name in the browser comes back as one curve per cell. /Grid says where each of those cells is, /Settings what the model was compiled with, and /Model/source keeps the model text itself, a line per value. A dual-porosity model can have more cells than a file should hold datasets, and past four thousand of them the rock layers are left out — the root attribute matrix_written says so when that happens, and the water is always there.

The Jacobian

Every rate law is differentiated symbolically with respect to each species it mentions. Nothing is differenced, so there is no step size to choose and no accuracy to lose. The chemistry is the same in every cell, so the matrix is one block per cell plus the transport terms, and those couple a species only to itself next door — which is why the state is ordered cell by cell rather than species by species: the matrix comes out banded instead of spread across its whole width.

The Jacobian tab draws the pattern of that matrix: the chemistry of each cell in the accent colour, the transport couplings to its two neighbours in the text colour, which is what puts the two thin diagonals either side of the blocks. In batch mode there is one block and it fills the picture. Nothing about the drawing is a guess — the worker folds the real pattern into a cell against itself and against each neighbour, which is the whole of it, and the caption says how much of each block is actually there.

Check Jacobian compares it against a central difference. Read the result with one thing in mind: an entry far smaller than the largest in its column cannot be measured that way at all — the change it makes is lost under the change the big ones make — so those are counted separately rather than called wrong.

Choosing a solver

NDF is the default and is the right first choice. BDF is the same integrator with the κ terms switched off. The Julia ports are independent implementations of other methods; FBDF is the one to try on a large transport problem, and RadauIIA5 when two others disagree. They are described on the canister radiolysis page, which shares them.

Every one of them is available wherever the page is: a run normally happens on a background worker, and where there is no worker to be had — a page opened straight off the disk cannot start one — it happens on the page's own thread instead. The answer is the same; what you lose is the progress bar and the Stop button, because the page cannot repaint until the run is over. Serving the folder rather than opening the file gets both back: python3 -m http.server in it, then the address it prints.

If NDF will not get through a model, try FBDF before anything else. The kind of chemistry this page is for is hard on a corrector: with rate constants up to 10¹⁶ against concentrations of 10^-9, the terms of a rate balance cancel over sixteen digits and the residual cannot be worked out to better than about a part in 10⁷. A Newton correction that has reached that floor stops shrinking, which an ordinary convergence test reads as divergence. This page tells the solver to take such a correction when it is within half the tolerance and the Jacobian is the one for that point; without it, the twenty published spent-fuel cases in resources/tests/rtm/verify-brum.py took millions of steps to cover a few seconds, and with it they run in a few thousand and agree with the reference code.

The examples

The Example picker on the Model tab loads any of them; the text is replaced, and a line under the picker says what the model is and where it comes from. They are grouped by what they are for:

One of them will not run as it stands. HS_case20's own database names UVIO2s+ in a stoichiometry whose rate law reads [UVIO2s+2]; they are different species, and the first appears nowhere else, so it is driven below zero — BRUM's own answer for it is negative too. Turn off keep concentrations non-negative under the solver's advanced settings and this page reproduces that result, negative species and all. The note is at the top of the model text.

Every example is compiled and run by resources/tests/rtm/test-model.js, so the picker cannot offer something broken, and the ones with a published answer are checked against it there.

Where the built-in model came from

Water radiolysis and uranium dissolution at a spent-fuel surface: 127 reactions over 36 species, converted mechanically from the skbrtm example set. Radiolysis enters as a source of G × dose rate, with G in mol/J and the dose rate in Gy/s, so with a litre of water weighing a kilogram the product is mol/L/s and no conversion factor is needed. The dose rate is written as a constant on each radiolysis line, which holds it uniform across the film — the real one is not, it falls away over the range of an alpha particle, about 30 µm. To put that profile in, move it into a parameter and let every radiolysis line use it:

<PARAMETERS>
DOSE  all  0.64*exp(-x/3.0E-5)

then delete the DOSE = 0.64 from the ends of those lines, so that the name resolves to the parameter instead.