Chemistry & Chemical Kinetics — Reaction-Diffusion, Combustion and Acid-Base Equilibria

Chemistry as a Dynamical System

Every chemical reaction is a dynamical system: concentrations change over time according to differential equations, equilibria emerge as fixed points, and nonlinear coupling between reactions can produce oscillations, waves, and self-organised spatial patterns. The same mathematical toolkit that describes a pendulum — phase portraits, stability analysis, bifurcation diagrams — applies directly to chemical kinetics.

This makes chemistry uniquely suited to simulation. Unlike a live laboratory experiment, a simulation lets you freeze time, reverse it, zoom into individual reaction fronts, and instantly adjust rate constants to see what happens at the edge of instability. The six simulations in this spotlight span from reactions happening in microseconds to biogeochemical cycles spanning centuries.

Layer 1: Self-Organisation — Turing Patterns

Reaction-Diffusion Simulation

In 1952, Alan Turing published a paper that had nothing to do with computers. He showed that a system of two chemicals — an activator that promotes its own production and an inhibitor that suppresses it — combined with different diffusion rates can spontaneously break spatial symmetry and create stable patterns: spots, stripes, and labyrinths. These "Turing patterns" explain the stripes on zebras, the spots on leopards, and the ridges on fingerprints.

PDEs (continuous):
  ∂U/∂t = D_U · ∇²U  -  U·V²  +  f·(1-U)
  ∂V/∂t = D_V · ∇²V  +  U·V²  -  (f+k)·V

U = activator concentration  [0, 1]
V = inhibitor concentration  [0, 1]
D_U, D_V = diffusion coefficients  (D_U > D_V for patterns)
f = feed rate (replenishment of U)
k = kill rate (removal of V)

Pattern regimes (f, k parameter space):
  f=0.035, k=0.065 → spots (leopard)
  f=0.060, k=0.062 → stripes (zebra)
  f=0.025, k=0.060 → moving spots (coral)
  f=0.014, k=0.054 → labyrinthine (fingerprint)
  f=0.039, k=0.058 → worms + spots

WebGL implementation:
  Ping-pong framebuffers (A↔B each frame)
  Fragment shader: Laplacian via 9-tap stencil
  Texture format: RG32F (U in R, V in G)
  Resolution: 512×512 → 262,143 cells updated per frame

The Reaction-Diffusion simulation implements the Gray-Scott system on a 512×512 WebGL grid updated in real time. Pan through the (f, k) parameter space and watch the pattern morph — from isolated spots to labyrinthine channels to travelling waves — without restarting the simulation. Paint initial seeds with your cursor and watch the pattern grow from your touch.

Layer 2: Rapid Kinetics — Combustion

Combustion Simulation

Combustion is the most energetically consequential chemistry on Earth — it powers every aircraft, ship, and most electricity grids. It is also one of the most complex: real hydrocarbon flames involve hundreds of simultaneous reactions, chain branching, free radicals, and thermal feedback. The simulation focuses on the key phenomenology: ignition, flame front propagation, and the conditions for explosive deflagration versus detonation.

Arrhenius rate law:
  k(T) = A · exp(-E_a / R·T)
  A   = pre-exponential factor [mol⁻¹·m³·s⁻¹]
  E_a = activation energy [J/mol]
  R   = 8.314 J/(mol·K)
  Doubling T from 300K to 600K can increase k by 10⁶×

H₂/O₂ chain branching (simplified):
  H + O₂ → OH + O    (chain branching, E_a=70 kJ/mol)
  O + H₂ → OH + H    (chain branching)
  H + O₂ + M → HO₂ + M  (chain termination at high P)

Explosion limits:
  First (low P):   chain termination at walls > branching → no explosion
  Second (~1 atm): branching > termination  → explosion
  Third (high P):  thermal runaway → explosion

Deflagration vs Detonation:
  Deflagration: flame travels at 0.5–5 m/s, subsonic
  Detonation:   shock-coupled flame at 1500–3000 m/s, supersonic
  Chapman-Jouguet velocity: D_CJ = √(2(γ²-1)·q)  (q = heat release)

The Combustion simulation models a 2D reactive flow using a simplified one-step Arrhenius mechanism. Set the fuel-to-air equivalence ratio, initial temperature, and pressure; then ignite with a spark. Watch the flame front accelerate as heat release from the reaction zone preheats the unburned mixture. Vary the equivalence ratio below the lean limit and see the flame extinguish; go rich and see incomplete combustion and soot precursor formation.

Layer 3: Equilibrium Chemistry — Acid-Base Reactions

Acid-Base Simulation

Acid-base chemistry underpins biochemistry (enzyme active sites, blood buffering), environmental science (ocean acidification, acid rain), and industrial processing (titrations, pH control in reactors). The mathematics of equilibrium is elegant: a handful of equilibrium constants fully determine the pH and speciation of even complex multi-component systems.

Water equilibrium:
  Kw = [H⁺][OH⁻] = 1×10⁻¹⁴  (at 25°C)
  pH + pOH = 14

Weak acid equilibrium:
  HA ⇌ H⁺ + A⁻
  Ka = [H⁺][A⁻] / [HA]       pKa = -log₁₀(Ka)

Henderson-Hasselbalch equation:
  pH = pKa + log₁₀([A⁻]/[HA])
  Buffer range: pKa ± 1  (10:1 to 1:10 ratio)

Titration curve (strong acid vs strong base):
  Before eq. pt:  pH = -log[H⁺]_excess
  At eq. pt:      pH = 7 (pure water)
  After eq. pt:   pH = 14 + log[OH⁻]_excess

Weak acid titration:
  Half-eq. pt:    pH = pKa  (buffer region centre)
  Equivalent pt:  pH > 7  (conjugate base hydrolysis)

Examples:
  Acetic acid:       pKa = 4.76  (vinegar)
  Carbonic acid:     pKa = 6.35  (blood, ocean CO₂)
  Ammonium:          pKa = 9.25  (cleaning products)
  Phosphoric acid:   pKa = 2.15, 7.20, 12.35  (biological buffer)

The Acid-Base simulation renders interactive titration curves for any combination of strong/weak acids and bases. Add virtual drops of titrant with a slider and watch the pH curve jump sharply at the equivalence point. Enable the indicator overlay — litmus, phenolphthalein, bromothymol blue — and see which colour range each covers. The buffer capacity graph shows exactly where the solution resists pH change most effectively.

Layer 4: Phase Transition — Crystal Growth

Crystal Growth Simulation

When a supersaturated solution or an undercooled melt begins to crystallise, the shape of the growing crystal is determined by a competition between surface energy anisotropy and diffusion-limited transport. The same mathematics governs snowflake dendrites, metal grain growth in alloys, and the mineral structures in geological formations.

Diffusion-Limited Aggregation (DLA):
  Random walker released from boundary
  Sticks on first contact with growing cluster
  Fractal dimension D ≈ 1.71 in 2D

Phase-field model (Kobayashi 1993):
  ∂φ/∂t = ε²·∇²φ + φ(1-φ)(φ - 1/2 + m)
  ∂T/∂t = ∇²T + ∂φ/∂t      (latent heat coupling)

  φ ∈ [0,1]: order parameter (liquid=0, solid=1)
  ε = interface thickness parameter
  m = driving force = (T_melt - T) / T_melt

Anisotropy (dendritic arms):
  ε(θ) = ε̄ · (1 + δ·cos(N·θ))
  N=4 → cubic crystal (6-fold → hexagonal)
  δ=0.05: weak anisotropy → compact crystal
  δ=0.15: strong anisotropy → pronounced dendrites

Snowflake criterion:
  N=6, δ≈0.05, D=0.5 → hexagonal snowflake dendrites
  Arm spacing ∝ √(D·t)  (diffusion-limited)

The Crystal Growth simulation lets you choose between DLA (aggregate random walkers) and a phase-field model. Watch DLA produce the iconic fractal cauliflower shape, then switch to phase-field and tune the anisotropy strength to grow snowflake-like dendrites. Changing the cooling rate shifts the morphology from compact grains to long, spindly arms — exactly the behaviour metallurgists control to engineer microstructure in castings.

Layer 5: Drug Dynamics — Pharmacokinetics

Pharmacokinetics Simulation

After swallowing a tablet, a drug follows a precisely quantifiable journey: absorption from the gut into plasma, distribution into body compartments, metabolism (primarily in the liver), and excretion (primarily by kidneys). This ADME process is governed by first-order differential equations that every pharmacologist and drug designer must understand.

One-compartment model (IV bolus):
  dC/dt = -k_e · C
  C(t) = C₀ · e^{-k_e·t}
  Half-life: t₁/₂ = ln(2)/k_e = 0.693/k_e

Two-compartment model:
  dC_p/dt = -(k_e + k_pt)·C_p + k_tp·C_t + F·D(t)/V_p
  dC_t/dt =   k_pt·C_p - k_tp·C_t

  C_p = plasma concentration
  C_t = tissue concentration
  k_pt = distribution to tissue, k_tp = redistribution
  F = bioavailability (oral: F < 1)

Key parameters:
  V_d = volume of distribution = dose / C₀
  CL = clearance = k_e · V_d  [L/h]
  AUC = ∫₀^∞ C(t)dt = dose·F / CL
  Steady-state (multiple doses):
    C_ss = F·dose / (CL·τ)     τ = dosing interval

Therapeutic window:
  MEC (min effective) < C_ss < MTC (min toxic)
  Design dosing regimen to stay within window

The Pharmacokinetics simulation implements the two-compartment model with configurable initial dose, bioavailability, clearance rates, and distribution volume. Select a dosing regimen — single IV bolus, repeated oral doses, or infusion — and watch the plasma concentration curve evolve. The therapeutic window overlay immediately shows whether the chosen regimen keeps the patient in the effective zone or overshoots into toxicity.

Layer 6: Planetary Chemistry — The Carbon Cycle

Carbon Cycle Simulation

The global carbon cycle is chemistry at planetary scale: carbon atoms flow between the atmosphere, ocean, terrestrial biosphere, and geological reservoirs on timescales from seconds (photosynthesis) to millions of years (carbonate rock formation). The mathematics of these fluxes — essentially a giant system of first-order ODEs with nonlinear feedbacks — determines atmospheric CO₂ and thus global temperature.

Four main reservoirs (GtC):
  Atmosphere:  ~870  (growing ~5 GtC/yr)
  Land biosphere: ~2600  (NPP ≈ 120 GtC/yr)
  Ocean surface:  ~900  (exchange ~90 GtC/yr)
  Deep ocean: ~37,000

Ocean chemistry linkage:
  CO₂(g) ⇌ CO₂(aq) → H₂CO₃ → H⁺ + HCO₃⁻ → H⁺ + CO₃²⁻
  pCO₂ = α · [CO₂(aq)]     Henry's law
  pH = -log[H⁺]           Ocean acidification: ΔpH = -0.1 since 1750

Box model (simplified):
  dC_atm/dt = E_fossil + E_land_use - F_land - F_ocean
  dC_ocean/dt = F_ocean - F_carbonate - F_export
  E_fossil ≈ 10 GtC/yr  (2020s)
  Airborne fraction: ~44% of emissions stay in atmosphere

Temperature feedback:
  ΔT = λ · ΔF = λ · 3.7 · ln(CO₂/CO₂_0)  (forcing in W/m²)
  λ ≈ 0.8 K/(W/m²)   (equilibrium climate sensitivity)

The Carbon Cycle simulation animates carbon flows between atmosphere, land, and ocean reservoirs in real time. Increase fossil fuel emissions, trigger deforestation, or stimulate ocean iron fertilisation, and watch CO₂ accumulate and temperature rise in response. The ocean acidification graph runs alongside, showing how carbonate chemistry limits the ocean's buffering capacity over decades.

Complete Chemistry Collection

Cross-Collection Connections

Chemistry bridges almost every other category on the platform. The Arrhenius rate law governing combustion chain branching is the same exponential in the Maxwell-Boltzmann simulation that explains why hot gas molecules react faster: both are Boltzmann factors. The reaction-diffusion pattern-formation mechanism reappears in the Cellular Automata simulation — Gray-Scott and the activator-inhibitor CA rule share the same mathematical structure. Pharmacokinetics is a direct application of the ODE methods covered in Learning #16 (ODEs in Biology), and the carbon cycle box model uses precisely the same compartmental approach. Crystal growth via DLA links to the fractal geometry in Spotlight #15 (Chaos) — both generate structures with non-integer fractal dimensions.

Algorithms & Methods in This Collection