WellRTSim

A MATLAB simulator for mineral scaling in geothermal wells. It couples transient one-dimensional two-phase wellbore flow to multi-mineral equilibrium and kinetic geochemistry through PHREEQC, and feeds the growing scale back into the conduit geometry and wall roughness, so the deposit and the flow evolve together.

The physics in one line: a four-equation drift-flux formulation transports mass, momentum and energy; PHREEQC speciates the liquid at each cell and time step; a kinetic rate law converts local supersaturation into mineral volume; the reduced diameter and raised roughness change the flow that produced it.

What it is for

Predicting where and when scale forms along a well, and how much production is lost as it does. Because the deposition rate depends on the flashing depth, which depends on the pressure profile, which depends on the deposit, these cannot be computed separately.

Requirements

• MATLAB — developed on R2026a.
• Parallel Computing Toolbox — PHREEQC is evaluated on a process pool (instances: in each case's chemistry.md). Set instances: 1 to avoid it.
• PhreeqcMatlab — https://github.com/simulkade/PhreeqcMatlab. Supplies IPhreeqc and the PHREEQC thermodynamic databases.

Water and steam properties use an IAPWS-IF97 implementation included here; no additional toolbox is needed for those.

Quick start

run('<path to PhreeqcMatlab>/startup.m');   % puts IPhreeqc and the databases on the path
cd('<path to WellRTSim>');
setenv('WELLRT_SIMDIR','Simulation/Krafla/KJ-9_bestfit');
main

main.m defaults to Simulation/Krafla/KJ-9_bestfit/ if WELLRT_SIMDIR is unset.

Headless, for batch runs:

WELLRT_HEADLESS=1 WELLRT_SIMDIR='Simulation/Tonkin/LO_T' matlab -batch main

To inspect results:

addpath('viewer');
result_view                       % pick a file interactively
result_view('Simulation/Krafla/KJ-9_bestfit')

Environment variables

variable	effect
WELLRT_SIMDIR	case directory to run
WELLRT_HEADLESS=1	disable all live plotting and drawnow
WELLRT_CHEM_AUDIT=1	print the per-step calcium balance; Ca fluid / wall must read 1.0000

Repository layout

main.m                  driver: init -> steady state -> time loop
functions/
  Hydrodynamics/        drift-flux fluxes, Newton iteration, timestep control
  Chemistry/            PHREEQC coupling, transport-reaction step, partitioning
  Wellspec/             casing, deviation and geometry handling
  init/                 parameter and state initialisation, steady solve, HDF5 setup
  IO/                   HDF5 output, wellhead pressure history
  Graphics/             live plots
viewer/                 standalone results viewer for the HDF5 output
Simulation/<field>/<case>/
  params.md             grid, boundary conditions, timestepping
  chemistry.md          feed composition, minerals, kinetics, database
  chemistry.pht         PHREEQC input template
  *.csv                 casing, deviation, feed zones, measured logs
  results/              indexed HDF5 output, created on demand

A case is defined entirely by files in its own directory — no code changes are needed to add one.

Cases included

case	description
Krafla/KJ-9_bestfit	Calcite scaling in Krafla well KJ-9 (Iceland), July 1979 – June 1980, calibrated against a caliper survey and measured discharge chemistry. The reactive case.
Tonkin/LO_T	Low-enthalpy production after Tonkin et al., adiabatic. Hydrodynamics only.
Tonkin/LO_T_heat	As above, with heat loss to the formation.
EastMesa/Mesa6	East Mesa well 6-1, steady-to-transient regression. Run via runMesa6TonkinTransient().

The Krafla calibration

Over 324 days between two cleanings, against six independent measurements:

quantity	model	measured
deposit above 670 m	6.86 m³	6.90 (caliper ID-427)
flash-front depth	644 m	625–730 m
initial flow	43.7 kg/s	43.5–46.1
wellhead Ca	1.86 mg/kg	1.88
wellhead CO₂	989 mg/kg	989
wellhead H₂S	76.3 mg/kg	73.5
flow at day 250	34.8 kg/s	29.8

The last row is the one genuine miss, and it is structural rather than under-tuned: reproducing that late collapse requires a large deposit in the slotted liner, which the measured wellhead calcium forbids. Both constraints are measurements, so whatever caused the collapse is outside the model.

The deposition timescale is not a free knob fitted to one number. It is pinned twice over — the wellhead calcium fixes it at τ ≈ 6×10⁷ kg·s·mol⁻¹, and the deposit volume above 670 m has its maximum at the same value. See Simulation/Krafla/KJ-9_bestfit/README.md for the full configuration, the mesh-convergence evidence, and the residual misfits; chemistry/INDEX.md there maps each published figure to the setup that produced it.

Output

Each run writes results/<prefix>_NNN.h5, index auto-incrementing:

group	contents
/meta	x (cell centres, height above the well bottom), Dp0 (as-built diameter)
/wellhead	time series: flow, pressure, enthalpy, temperature, quality, phase densities and velocities, per-component concentrations
/profiles/<time>	per-cell snapshot: P, T_C, Dp, alpha_g, quality, mass fluxes, element_*, SI_*, scaleFrac_*
/chemistry	per-step chemistry diagnostics
/inputs	verbatim copies of params.md, chemistry.md and chemistry.pht

Two conventions worth knowing when post-processing:

• meta/x is height above the well bottom, not depth. Convert with depth = L - x.
• element_* are liquid-phase mass fractions; ppm = value × 10⁶. They are the basis the saturation indices are computed on, which is not the same as a total-discharge concentration when the fluid has flashed.
• scaleFrac_* is the incremental mineral fraction for each saved interval, so cumulative per-mineral volume must be accumulated over the profile sequence.

Because /inputs embeds the setup files, every result file records exactly the configuration that produced it.

Limitations

• Pure-water carrier fluid. Dissolved solids and non-condensable gases do not affect the two-phase properties. Do not apply above roughly 5 wt % TDS or 5 mol % CO₂ without extending the property model — this is the highest-priority gap.
• Lumped deposition kinetics. One characteristic time per mineral stands in for nucleation, transport to the wall, attachment and erosion. It is an effective, site- and mineral-specific quantity calibrated against field data, not a transferable constant.
• Gas partitioning. With partition: 2 the molecular-CO₂ distribution coefficient is applied to total dissolved carbonate, which over-strips carbon as degassing converts it to bicarbonate. It does not prevent the Krafla case reproducing the measured outlet CO₂, but it is a known defect; see that case's notes.
• Deliverability slope. The Krafla case gives dQ/dP₀ ≈ −0.55 kg/s/bar against −3.2 to −6.1 measured — a long-standing, unresolved discrepancy.
• Flow-regime-dependent kinetics and transport of suspended solids are not modelled.

This repository and the archived release

This repository is the development source. Simulation output is not tracked here — results/*.h5 is git-ignored, because the files are large and are regenerated by running a case.

A separate Zenodo deposit carries a citable snapshot including the HDF5 result files behind the published figures, so the figures can be re-plotted without re-running anything. Use this repository to run and modify the code; cite the Zenodo release.

One difference between the two: in the archived copy the PHREEQC database is given as a bare phreeqc.dat in each chemistry.md, resolved on the MATLAB path. Here it may be an absolute path pointing at a local PhreeqcMatlab checkout — change that line to suit your installation, or to select a different database.

Licence

CC BY-NC-SA 4.0 for academic and research use — see license.md. For commercial licensing contact oleg.melnik@earth.ox.ac.uk or oemelnik@gmail.com.

Citation

Cite the archived release rather than a git revision. The DOI below is reserved on Zenodo and becomes active when the deposit is published.

@software{wellrtsim,
  author  = {Melnik, Oleg},
  title   = {WellRTSim: A MATLAB Wellbore Reactive Transport Simulator},
  year    = {2026},
  version = {2.0},
  doi     = {10.5281/zenodo.21476565},
  url     = {https://github.com/crystalworkshop/WellRTSim}
}