@@ -363,6 +363,9 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | `airfoil_id`         | Integer | Index into `ib_airfoil` array for NACA airfoil geometry patches. |
 | `model_id`           | Integer | Index into `stl_models` array for STL/OBJ geometry patches. |
 | `slip`               | Logical | Apply a slip boundary |
+| `thermal_bc`         | Integer | Thermal boundary-condition selector: 0 = zero-normal-gradient temperature, 1 = prescribed wall temperature, 2 = reacting surface energy balance. |
+| `Twall`              | Real    | Prescribed wall temperature used when `thermal_bc = 1`. |
+| `surface_reaction`   | Integer | Heterogeneous surface-reaction flag: 0 = disabled, 1 = enabled. |
 | `moving_ibm`         | Integer | Sets the method used for IB movement. |
 | `vel(i)`             | Real    | Initial velocity of the moving IB in the i-th direction. |
 | `angular_vel(i)`     | Real    | Initial angular velocity of the moving IB in the i-th direction. |
@@ -413,6 +416,12 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `slip` applies a slip boundary to the surface of the patch if true and a no-slip boundary condition to the surface if false.
 
+- `thermal_bc` selects the thermal immersed-boundary condition. A value of 0 applies a zero-normal-gradient temperature condition, 1 prescribes the wall temperature using `Twall`, and 2 solves the reacting-surface energy balance for the surface temperature. The `thermal_bc = 2` option requires `surface_reaction = 1`. A non-zero `thermal_bc` requires `chemistry = T` and cannot be combined with `inj_species > 0`, since the thermal condition is applied by the chemistry ghost-state reconstruction, which an injecting surface bypasses.
+
+- `Twall` specifies the prescribed surface temperature when `thermal_bc = 1` and must be positive in that case.
+
+- `surface_reaction` enables heterogeneous surface chemistry when set to 1. Surface reactions require `chemistry = T` and cannot be combined with `inj_species > 0`.
+
 - For STL/OBJ geometry (geometry 5 or 12), set `model_id` to index into the `stl_models` array and specify `model_filepath`, `model_scale`, `model_translate`, and `model_threshold` on that entry.
 
 - `moving_ibm` sets the method by which movement will be applied to the immersed boundary. Using 0 will result in no movement. Using 1 will result 1-way coupling where the boundary moves at a constant rate and applied forces to the fluid based upon its own motion. In 1-way coupling, the fluid does not apply forces back onto the IB. Using 2 will result in 2-way coupling, where the boundary pushes on the fluid and the fluid pushes back on the boundary via pressure and viscous forces. If external forces are applied, the boundary will also experience those forces.
@@ -1249,13 +1258,17 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 | `chem_params%%adap_substeps`   | Logical | Per-rank adaptive sub-step count driven by local stiffness  |
 | `chem_params%%reaction_substeps_max` | Integer | Sub-step ceiling when `adap_substeps` is enabled       |
 | `cantera_file`                | String  | Cantera-format mechanism file (e.g., .yaml)              |
+| `surface_cantera_file`        | String  | Cantera-format mechanism file for heterogeneous surface chemistry |
+| `surface_phase`               | String  | Cantera interface phase name for heterogeneous surface chemistry |
 
 - `chem_params%%transport_model` specifies the methodology for calculating diffusion coefficients and other transport properties, `1` for mixture-average, `2` for Unity-Lewis
 - `chem_params%%reaction_substeps` controls how the reaction source is integrated. With `0` (default) the net production rates are added to the flow right-hand side and advanced by the flow time stepper (fine for hydrogen). With a value `> 0`, the reaction is instead integrated by operator splitting after each flow update: every cell's constant-density, constant-internal-energy reactor is advanced over the timestep with that many sub-steps of an **α-QSS** (quasi-steady-state) integrator — a matrix-free, Jacobian-free predictor–corrector (Mott/CHEMEQ2) that splits the net rate into creation/destruction parts and applies a Padé α-weighting, so it stays stable on stiff mechanisms where an explicit source diverges. This decouples the (often much faster) chemical timescale from the flow timestep and is required for stiff mechanisms — e.g. hydrocarbons such as GRI-Mech methane, which otherwise diverge on the first step
 - `chem_params%%adap_substeps` (default `F`) makes each rank choose its α-QSS sub-step count per flow step from the largest chemical stiffness among its own cells: the count sits at `reaction_substeps` (the floor) in inert or burned gas and rises toward `reaction_substeps_max` (the ceiling) only across the reaction front. It uses no MPI collectives. When enabled, `reaction_substeps >= 1` and `reaction_substeps_max >= reaction_substeps` are required
 
 - `cantera_file` specifies the chemical mechanism file. If the file is part of the standard Cantera library, only the filename is required. Otherwise, the file must be located in the same directory as your `case.py` file
 
+- `surface_cantera_file` and `surface_phase` specify the Cantera mechanism file and interface phase used for heterogeneous surface chemistry. These parameters must be specified together when a surface mechanism is used.
+
 MFC generates and compiles the mechanism's Fortran routines itself. Supported mechanism features and the Cantera-only mixing-layer initialization are described in @ref thermochemistry "Thermochemistry implementation".
 
 ### 18. Chemistry-Specific Boundary Conditions
 
@@ -104,6 +104,15 @@ Each of these fails quietly rather than loudly.
   carrying both labels at once and also matches nothing. An empty selection then exits
   **143**, which reads like an external kill rather than an empty filter. Pass UUIDs when
   you want the union of several groups.
+- **`Chemistry` is the one label not read off the trace.** Any case with ``chemistry='T'``
+  answers to it, because it selects a *build* and not just a test: Frontier AMD's GPU lane
+  compiles its chemistry binaries in a separate SLURM job invoked with `-o Chemistry`
+  (`.github/workflows/common/build.sh`) and then tests with `--no-build`, so a chemistry
+  case the filter misses is never compiled there and fails with a missing binary. Examples
+  are auto-registered as `<dim> -> Example -> <dirname>` and so can never carry the label
+  by hand. If you add a chemistry case, you get this for free; do not re-add the label to
+  a trace to compensate, since the UUID is a hash of the trace and renaming orphans the
+  golden directory.
 - **Sibling `define_case_d` calls at the same stack level are never combined.** Two switches
   that only matter together therefore get no effective coverage unless one is pushed onto
   the stack and the other defined beneath it — `avg_state=1`, for instance, is only read
 
@@ -19,6 +19,7 @@ renders the Fortran template. Neither is a runtime solver dependency.
 | `thermochem/expressions.py` | NASA7, reaction-rate, equilibrium and transport expressions |
 | `thermochem/fortran.py` | Supported-feature checks, expression formatting and module generation |
 | `thermochem/module.fpp.mako` | Fortran interface and numerical routines |
+| `thermochem/surface.py`, `surface.fpp.mako` | Heterogeneous surface-chemistry checks and module |
 | `thermochem/fingerprint.py` | Mechanism and generator content identities for build reuse |
 | `run/input.py` | Mechanism resolution and generation for each target |
 
@@ -58,6 +59,22 @@ can therefore differ in roundoff from the former double-literal expressions.
 Troe falloff guards the logarithm at zero reduced pressure, including when a
 compiler evaluates both arguments of a Fortran `merge` expression.
 
+# Surface chemistry
+
+Reacting immersed boundaries (`surface_cantera_file`, `surface_phase`) also use
+`generate_surface_fortran`, which writes `m_surface_thermochem.fpp` for the simulation
+target. It shares the gas generator's rate-coefficient and NASA7 expressions and its
+`wp`/`$:GPU_ROUTINE` conventions, and calls `m_thermochem` for concentrations and gas
+enthalpies. The module provides the gas-species net production rates and the heat
+released by the surface reactions; without a surface mechanism both return zero.
+
+The supported model is a bulk solid reacting with gas species: irreversible
+(interface-)Arrhenius reactions whose rates depend on gas concentrations, with
+explicit orders allowed. Adjacent bulk phases, such as graphite, have unit activity
+and enter the reaction enthalpy through their NASA7 thermodynamics at the reference
+pressure. Sticking, Blowers-Masel and coverage-dependent rates, reversible surface
+reactions, and surface-site species are rejected before emission.
+
 # Mixing-layer initial conditions
 
 The temporal 2D, spatial 2D and temporal 3D reacting mixing-layer examples share
@@ -92,7 +109,8 @@ After bootstrapping the toolchain, run:
 
 ```sh
 PYTHONPATH=toolchain build/venv/bin/pytest -q \
-    toolchain/mfc/test_thermochem.py toolchain/mfc/test_flamelet.py
+    toolchain/mfc/test_thermochem.py toolchain/mfc/test_flamelet.py \
+    toolchain/mfc/test_surface_chemistry_codegen.py
 ./mfc.sh test --no-mpi -j 8 --only Chemistry
 ```
 
@@ -106,5 +124,9 @@ with OpenACC and OpenMP directives. They run the generated source through Fypp w
 MFC's macros, as the build does. The directive tests
 execute on the host; they do not validate GPU offload on accelerator hardware.
 
+The surface tests compile the example carbon mechanism's module and compare gas
+production rates and reaction heat with Cantera's interface kinetics, check the
+gasified carbon mass balance, and cover each rejected rate law.
+
 Initialization tests check stream limits, normalization, elemental composition,
 enthalpy, density and a burning hot profile for both mixing-layer mechanisms.