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 button at the right end of the tabs gives the page the whole window, without the site's header and footer; the kvot mark at the left of the tabs then goes to the home page, and the switch beside the button changes between light and dark. The page remembers the choice.

The model text

Six sections, of which two are needed: <SPECIES> and <REACTIONS>, and any number of named tables, <TABLE name>. Everything after a # is a comment. A rock matrix beside the column is a few more settings; see A rock matrix. Wherever a number goes it may be written as arithmetic on numbers — 365*86400, 1.33*10**12, 2/3 — and is worked out when the model compiles. Open… loads a model text from a file, and so does dropping the file anywhere on the page; either also reads skbrtm's database files, see Opening skbrtm and BRUM databases.

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. With a rock matrix a species also takes De= (or Dm=), Kd= or Rm=, and Kdf=: see A rock matrix below.

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, and groups in round brackets, one level deep: Fe(OH)3, UO2(CO3)3-4. Square brackets are what a rate law puts round a name, so a name cannot contain those.

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's exponent may be anything that does not depend on a concentration: a number, arithmetic such as (2/3), a constant of the line or a parameter.

Fe+2 + 0.25 O2 => Fe(OH)3, r = (kr/kh)*[O2]*[Fe+2]*[OH-]^2, kr = 1.33*10**12, kh = 1.3*10**-3
A => B, r = k*c^(2/3)*(1 - [B]/cmax), k = 5e-13, c = [A], cmax = 1.8e-6

A constant may be arithmetic, and may use the constants named before it on its line — kf = 2, kb = kf/4. A constant written c = [A] is an alias: in the law, c is the concentration of A, exactly as if [A] had been written there. Both are how skbrtm writes its databases, and they mean here what they say. (skbrtm itself pastes each constant into the law as text, so it reads (kr/kh) above as 1.33*10**12/1.3*10**-3 — a million times less.)

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.

How many cells a semi-infinite end needs past it is a matter of how far upstream the dispersion reaches. It dies away by D/(D + uΔx) per cell, which for a Peclet number and no D of a species' own is (2N − Pe)/(2N + Pe), and unless EXTRA_CELLS says otherwise they are the fewest that bring that under a tenth: five for twenty cells at Pe 10, as Kompartment counts them, and none where the grid's own dispersion is all there is. They are the column going on, so a line of <PARAMETERS> or <INITIAL> that reaches the last cell reaches them too, and a rock matrix is behind them as well. The charts leave them out.

A Danckwerts outlet, free, closes the column at its end instead, and lets out 4a/(1 + a)² of what the semi-infinite column does at a frequency where a = √(1 + 4t_wg/Pe). Against FARF31's exact transfer function, the far-field path below is 2.4 % out where most of a release is with a semi-infinite end, and 6.8 % with a free one. More cells bring the first down by four per doubling, to 0.17 % on 160, and make the second worse, 9.9 % on 160: that column closes on Danckwerts' answer, which is not FARF31's.

What crosses an end. On Against time the Cell picker ends with flux through the left end and through the right end: what crosses the plane at 0 or at LENGTH, per m² of the water's cross-section and positive to the right, in concentration times metres per unit of time. At an outlet that is the release — advection alone at a free end, advection and dispersion at a semi-infinite one. It is worked out from the stored run with exactly the terms the equations put through that face, so what the column gains is what crosses its left end less what crosses its right, to round-off.

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. Say anything about the rock and every cell of the column gets a chain of stagnant cells behind it, twelve unless MATRIX_CELLS says otherwise (0 for none), reaching MATRIX_DEPTH metres into the rock. A far-field path is written the way the SKB reports describe one, and here as Kompartment's far-field block starts one:

<SETTINGS>
MODE = transport
CELLS = 20                  # along the path
LENGTH = 100                # m; any length: TRAVEL_TIME sets the flow
TRAVEL_TIME = 100           # the water's travel time t_w, in the time unit
PECLET = 10
LEFT = robin                # water arriving at left= per m3
RIGHT = semi-infinite       # the rock goes on past the release point
MATRIX_DEPTH = 12.5         # m into the rock
MATRIX_POROSITY = 0.0018
MATRIX_DENSITY = 2700       # kg/m3, the default
TRANSPORT_RESISTANCE = 1E5  # F = a_w t_w: 1000 m2 of wall per m3 of water here
TIME_UNIT = year

<SPECIES>
Cs  0  left=1  De=1E-4  Kd=0.01    # De in m2/a, Kd in m3/kg

Each species has its De= (or Dm=, the same thing), the effective diffusivity in the rock, and a capacity there: Kd=, the partition coefficient, makes it ε_m + ρ_mK_d with MATRIX_POROSITY and MATRIX_DENSITY, and Rm= gives it outright and may name a parameter, so that it varies from cell to cell. With neither the species has the porosity alone, which is one that does not sorb; with no De it never enters the rock. Kdf= is sorption on the fracture coating, which retards the species in the water by R = 1 + K_d,fa_w; R= says it outright.

The flow-wetted surface a_w, the area of fracture wall a cubic metre of flowing water touches, is given one of three ways: WETTED_SURFACE itself, in m²/m³; APERTURE (2b) for a fracture of two walls, as 2/aperture; or TRANSPORT_RESISTANCE, the F of the SKB reports, as F/t_w. Two of them given have to agree. TRAVEL_TIME gives the flow as the water's travel time; without it t_w is LENGTH over the pore-water speed. A unit release into the path is a robin left end with left at one over the pore-water speed, or a source => X, k = S with S 0 1/Δx fracture in <PARAMETERS>; the release out is what crosses the right end.

The layers grow geometrically with depth, because the gradient is steepest at the wall and all but flat far in, and MATRIX_GRID says how:

MATRIX_FIRST sets the first layer instead, whichever the grid. The far face of the last layer is closed: the depth is a depth, not a boundary. Layers that grow by more than 2.5 are coarse at depth, and the page then says how many would keep the growth to 2. On the U-238 chain among the examples, Po-210's half-life of 0.38 years sizes the first layer at 5 µm: twelve layers to 4.5 m grow by 3.4 and put the releases 4 % out, twenty grow by 2.0 and are within 0.2 %.

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, and makes as much of what is on the right: a daughter born of a parent that sorbs hard is born as many atoms as the parent lost, and then takes up its own share between water and rock, as Kompartment and FARF31 have it. Such a reaction has one species on the left, or none.

Cs => , k = 7.3E-10, on = inventory     # decay: of everything, sorbed included
U234 => Th230, k = 2.8234E-6, on = inventory   # a chain: Th-230 from all of the U-234
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.

How fine. With matched layers, twelve do as well as twenty: what is left is the cells along the path, and they converge on the exact answer at second order. On the path above, Ra-226 (De=1E-4 Kd=4.5E-4, decaying with its half-life of 1600 years, as in the far-field example) arrives at 2.2×10⁻⁵ of a steady release on 20 cells, where FARF31's exact answer is 1.72×10⁻⁵; 40 cells give 1.84×10⁻⁵ and 80 give 1.75×10⁻⁵. A strongly held-back nuclide wants more cells than a tracer does.

What was checked. Kompartment's far-field path, from its own code, on five paths: its matched layers bit for bit, its count of cells past the release point, every entry of its transport matrix and of its chain's ingrowth to one part in 10¹⁵, and its release, read at the same plane with the same weights. The model's transfer function against FARF31's exact one, from far below the path's own frequencies to where it has fallen to e⁻¹⁵. With MATRIX_GRID = reference, the layer thicknesses and every rate 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¹⁴, and 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; a closed box that conserves mass to 10^-14 and shares it between water and rock exactly by capacity; and a decay chain that keeps its atoms in water and rock alike.

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.

A release history

A source that follows the clock is a reaction with nothing on the left whose rate reads a table at the time, t: what a near field releases into a far-field path, year by year, as another code worked it out or as a report tabulates it.

<TABLE release>
t           U238      U234          # the first row names the columns; t in years
0           0         0
1000        2.0E-3    1.5E-3
1.0E4       1.1E-3    9.0E-4
1.0E5       0         0

<PARAMETERS>
IN  all  0
IN  0    1/w  fracture              # into the first cell: one per m2 of water's cross-section

<REACTIONS>
=> U238, r = IN*release(t, U238)
=> U234, r = IN*release(t, U234)

IN turns a release per m² of the water's cross-section into one per m³ of the first cell's water; what crosses the right end, on Against time, is then per m² too, and the two compare directly. The time is in the model's own unit, and a table is held at its last row beyond it: end a history with zeros if it ends. The argument may be any arithmetic on t and the line's constants — release(t - 1000) to start it a thousand years late, or release(t/3.156e7) to read a history in years in a model in seconds.

The run is started again at every corner of such a table — every row, moved and scaled with the argument when that is a straight line in t. A stiff solver that has grown its steps over a quiet stretch would otherwise step clean over a short pulse and never see it: a release of ten years in five thousand comes out as nothing at all without the restarts, and whole with them. Kompartment does the same at its switch times. The panel says how many corners there are, and the status line how many the run was started at; a history of a few hundred rows costs a few thousand steps.

Through Kompartment's default far-field path, a history that rises, holds and falls over twenty thousand years comes out as FARF31 convolves it with the path's exact response, to the accuracy of the cells along the path: for Ra-226, which the rock holds back by five orders, the peak is 30 % high on 20 cells, 7 % on 40 and 1.8 % on 80.

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 whose exponent does not depend on a concentration — [A]^2, [A]^(2/3), [A]^n with n a constant or a parameter — 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, any parameter, t, the time in the model's own unit, and a table read at that time — see A release history below; a constant on the reaction line wins if a parameter has the same name, and no parameter may be called t. A parameter, and a starting value in <INITIAL>, may use x, the cell centre in metres, w, the cell's width, xl and xr, its two faces, i, the cell index, the tables, 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
U_site  0     2.1e-4/w/1000       # per m2 of surface, spread through cell 0
B       5-9   0      fixed        # held at 0 in cells 5 to 9, and only there
X       0     1      fixed        # a reservoir: cell 0 held at 1, feeding the rest

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. The value may be an expression in what a parameter can use, worked out in each cell — so an amount per m² of surface goes in as itself divided by w, and stays the same amount per m² whatever the grid.

fixed at the end of the line holds the species in those cells and nowhere else, as fixed on its species line holds it everywhere: its row is zero there, so nothing — reaction, diffusion, a face — can move it, and it still feeds its neighbours and drives its reactions. That is a reservoir cell, or a surface species that must not appear in the water. With no value, U_site 1-9 fixed holds it where it already is. A species held in some cells only cannot be put on the <EQUILIBRIUM> relations by EQUILIBRATE, and the page says so rather than moving it anyway.

<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, w, the cell's width, xl and xr, its faces, i, the cell's index, the tables, 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.

<TABLE name>

<TABLE dose_alpha log>
# distance from the surface (m)   dose rate (Gy/s)
1.0e-7    2.119
1.0e-6    2.02
1.0e-5    0.93
4.0e-5    5.6e-3

<PARAMETERS>
DOSE  all  dose_alpha(x)          # or interp(dose_alpha, x), as FACSIMILE writes it

x increasing and a value, used in a parameter or a starting value as name(x) — at any argument, name(x - 1e-6) included — and in a rate law as name(t), read at the clock. It is interpolated linearly between rows and held at its end values beyond them. After the name, log interpolates the logarithm of the values instead, which follows a profile that dies away exponentially, a dose rate, far better between rows that are far apart, and loglog the logarithm of the values in the logarithm of x: a power law between rows, which is how FARF31 reads a release history. Their values have to be above zero, and loglog's x as well; below the first x the first value is held, so a time of zero reads the first row.

A table may have several columns of values. A first row whose first word starts with a letter names them — t  U238  U234 — and a column is read by its name or its number from 1: release(t, U238), release(t, 2). A table of one column needs neither. Pasting a release history from another code's output is the reason: one table, a column per nuclide.

A rate law reads a table at the clock and never at a concentration: the table's slope would then be part of the Jacobian. A table of place is read in a parameter, as above, and the parameter used in the rate law.

This is how a tabulated dose-rate profile goes in as it is, where skbrtm fits a function to it first — and a fit is only as good as where it was started: skbrtm's beta and gamma fits come out as straight lines, 30 % low at the surface.

A surface

A dissolving surface — sites, a coating, a reacting film — is written as species that live in the first cell, stated per m² of surface and spread through that cell: U_site 0 2.1e-4/w/1000 in <INITIAL>, with no D so that they stay put. The amount per m² is then right on any grid.

What the chemistry does there is another matter. A rate law is written in concentrations, so a reaction between two surface species goes as the square of something divided by the cell's width, and when that width comes from the grid, refining the grid moves the answer: in skbrtm's fuel-dissolution example the dissolved uranium per m² halves each time the cells are halved. SURFACE_LAYER pins the first cell's thickness, making it a property of the model — a reaction zone, as the batch model's 40 µm alpha layer is — and lays the rest of the column out behind it. Refine the grid then and the surface does not change.

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.

Opening skbrtm and BRUM databases

Open… takes skbrtm's database files as well as a model text: select any of reaction.in, solutions.in, diffusion.in, the sourceterm and doserate files, and the Python script that runs them — or drop them on the page, or drop an example's whole folder, script and databases together; what is in it besides those is passed over. The script says which files are read, in what order, and for how long; without it every database selected is read and TEND is one second, to be set by hand. The files are read as skbrtm reads them — a SOLUTION sets every species in its cells, a species the DIFFUSION block does not name diffuses at 10^-9 m²/s, names are the same whatever their case — and become one model text, with rate constants per second as skbrtm integrates them.

Where skbrtm would read the text differently from what it says, the model follows the text and the difference is written at the top of it. Every rate law is read both ways, and one whose arguments skbrtm pastes in so that the value changes is named, with the factor. A Dirichlet end holds the face here, where skbrtm holds the first cell; a REGRESSION becomes a <TABLE> of its own data, where skbrtm fits a function to it; a constant species becomes fixed, in the cells its SOLUTION names; and a lone cell_id means that cell, where skbrtm puts it in cell 0.

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 list above it, at one cell, or what crosses one end of the column (see The two ends of the column). 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 in the model's own unit — and leaving it empty spreads eight over the run.

Gradient puts the whole run in one picture: time across, distance up, 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 list, 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, /Ends/left/… and /Ends/right/… what crosses each end of the column, /Settings what the model was compiled with, and /Model/source keeps the model text itself, a line per value. The cells past a semi-infinite end are under /Cells too, after the column's own, and say so. 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. Its BDF formulas switch, under Advanced settings, runs the same integrator with the κ terms off. The menu is Kompartment's, under each method's own name, and the canister radiolysis page shares it.

Auto starts on an explicit method — Tsit5, or Vern7 below a relative tolerance of 10⁻⁵, at a tenth of the tolerances — watches every step for stiffness as DifferentialEquations.jl's default does, and hands the rest of the run to the NDF once it finds it. A model that never turns stiff, A → B → C or the Brusselator, stays explicit and costs less than on the NDF; a transport column turns stiff in its first steps and costs what the NDF does. Rosenbrock 2-3 and Dormand–Prince 4-5 are Kompartment's own one-step methods: the first forms its Jacobian and factorises it every step — once a run, on a model whose Jacobian cannot change, such as a tracer column or a decay chain — and the second is explicit, for a model that is not stiff. The Julia ports are independent implementations of other methods: FBDF is the one to try on a large transport problem, RadauIIA5 when two others disagree, DefaultODEAlgorithm to compare with Julia's own default, and Tsit5 and Vern7 for a model that is not stiff.

Two settings come from Kompartment too. First step sets the step the run starts with, where the derivative at the start does not say what is about to happen; 0 lets the solver choose. Norm control, the NDF's, judges the error against the size of the whole solution rather than each species in each cell against its own, as MATLAB's NormControl does — looser wherever the species differ by orders of magnitude, so it is off.

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, a reaction set of this kind can take millions of steps to cover a few seconds, and with it runs in a few thousand.

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:

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.