diff --git a/TODO.md b/TODO.md
index f921824..7a4f5a0 100644
--- a/TODO.md
+++ b/TODO.md
@@ -8,7 +8,7 @@
 - [x] Flesh out condenser behavior: vacuum pump limits, cooling water temperature coupling, and dynamic back-pressure with fouling.
 - [ ] Dashboard polish: compact turbine/generator rows, color critical warnings (SCRAM/heat-sink), and reduce repeated log noise.
 - [ ] Dashboard multi-page view (F1/F2): retain numeric view on F1; future F2 schematic should mirror real PWR layout with ASCII art, flow/relief status, and minimal animations; add help/status hints and size checks; keep perf sane.
-- [ ] Core thermal realism: add a simple radial fuel model (pellet/rim/clad temps) with burnup-driven conductivity drop and gap conductance; upgrade CHF/DNB correlation (e.g., Groeneveld/Chen) parameterized by pressure, mass flux, and quality, then calibrate margins.
+- [x] Core thermal realism: add a simple radial fuel model (pellet/rim/clad temps) with burnup-driven conductivity drop and gap conductance; upgrade CHF/DNB correlation (e.g., Groeneveld/Chen) parameterized by pressure, mass flux, and quality, then calibrate margins.
 - [ ] Transient protection ladder: add dP/dt and dT/dt trips for SG overfill/depressurization and rod-run-in alarms; implement graded warn/arm/trip stages surfaced on the dashboard.
 - [ ] Chemistry & fouling: track dissolved oxygen/boron and corrosion/fouling that degrade HX efficiency and condenser vacuum; let feedwater temperature/chemistry affect steam purity/back-pressure.
 - [ ] Balance-of-plant dynamics: steam-drum level controller with shrink/swell, feedwater valve model, turbine throttle governor/overspeed trip, and improved load-follow tied to grid demand ramps.
diff --git a/src/reactor_sim/state.py b/src/reactor_sim/state.py
index 6e498b5..0856734 100644
--- a/src/reactor_sim/state.py
+++ b/src/reactor_sim/state.py
@@ -21,6 +21,8 @@ class CoreState:
     power_output_mw: float  # MW thermal
     burnup: float  # fraction of fuel consumed
     clad_temperature: float | None = None  # Kelvin
+    pellet_center_temperature: float | None = None  # Kelvin, peak centerline surrogate
+    gap_conductance: float = 1.0  # surrogate factor 0-1
     dnb_margin: float | None = None  # ratio to critical heat flux surrogate
     subcooling_margin: float | None = None  # K until boiling
     xenon_inventory: float = 0.0
@@ -32,6 +34,8 @@ class CoreState:
     def __post_init__(self) -> None:
         if self.clad_temperature is None:
             self.clad_temperature = self.fuel_temperature
+        if self.pellet_center_temperature is None:
+            self.pellet_center_temperature = self.fuel_temperature
         if self.dnb_margin is None:
             self.dnb_margin = 1.0
         if self.subcooling_margin is None:
@@ -131,6 +135,8 @@ class PlantState:
         core_blob.pop("dnb_margin", None)
         core_blob.pop("subcooling_margin", None)
         core_blob.setdefault("clad_temperature", core_blob.get("fuel_temperature", 295.0))
+        core_blob.setdefault("pellet_center_temperature", core_blob.get("fuel_temperature", 295.0))
+        core_blob.setdefault("gap_conductance", 1.0)
         turbines_blob = data.get("turbines")
         if turbines_blob is None:
             # Compatibility with previous single-turbine snapshots.
diff --git a/src/reactor_sim/thermal.py b/src/reactor_sim/thermal.py
index 3fd5a91..56663be 100644
--- a/src/reactor_sim/thermal.py
+++ b/src/reactor_sim/thermal.py
@@ -119,13 +119,18 @@ class ThermalSolver:
         heating = (0.002 * max(0.0, power_mw) + 0.01 * max(0.0, residual_power_mw)) * dt
         cooling = 0.025 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt
         core.fuel_temperature += heating - cooling
-        # Simple clad/fuel split: clad lags fuel and is cooled by coolant.
+        # Simple radial split: fuel -> clad conduction with burnup-driven conductivity drop and gap penalty.
+        gap_penalty = 1.0 + 2.0 * core.burnup
+        conductivity = 0.03 / gap_penalty
+        conduction = conductivity * max(0.0, core.fuel_temperature - (core.clad_temperature or primary.temperature_out)) * dt
         clad = core.clad_temperature or primary.temperature_out
-        conduction = 0.02 * max(0.0, core.fuel_temperature - clad) * dt
-        clad_cooling = 0.05 * max(0.0, clad - primary.temperature_out) * dt
+        clad_cooling = 0.06 * max(0.0, clad - primary.temperature_out) * dt
         clad = max(primary.temperature_out, clad + conduction - clad_cooling)
         core.fuel_temperature = max(primary.temperature_out, core.fuel_temperature - conduction)
-        # Keep fuel temperature bounded and never below the coolant outlet temperature.
+        # Peak pellet centerline surrogate based on power and burnup conductivity loss.
+        peak_delta = 20.0 * (core.power_output_mw / max(1.0, constants.NORMAL_CORE_POWER_MW)) * (1.0 + core.burnup)
+        core.pellet_center_temperature = min(constants.MAX_CORE_TEMPERATURE, core.fuel_temperature + peak_delta)
+        # Keep temperatures bounded and never below coolant outlet.
         core.fuel_temperature = min(core.fuel_temperature, constants.MAX_CORE_TEMPERATURE)
         core.clad_temperature = min(clad, constants.MAX_CORE_TEMPERATURE)
         core.subcooling_margin = max(0.0, saturation_temperature(primary.pressure) - primary.temperature_out)
@@ -176,13 +181,15 @@ class ThermalSolver:
         self._resolve_secondary_state(secondary)
 
     def _critical_heat_flux(self, primary: CoolantLoopState) -> float:
-        """Rough CHF surrogate using mass flux and pressure."""
-        # Use a coarse mass-flux and pressure scaling to emulate higher CHF with more flow/pressure.
-        mass_flux = max(1.0, primary.mass_flow_rate / 50.0)  # kg/m2-s surrogate
-        flux_factor = max(0.5, min(3.0, (mass_flux / 200.0) ** 0.5))
-        pressure_factor = 0.5 + 0.5 * min(1.5, primary.pressure / max(0.1, constants.MAX_PRESSURE))
-        base_chf = 1.0e7  # W/m2 surrogate
-        return base_chf * flux_factor * pressure_factor
+        """CHF surrogate with pressure, mass flux, and subcooling influence (Chen-like)."""
+        area = self._core_surface_area()
+        mass_flux = max(1.0, primary.mass_flow_rate / max(1.0, area))  # kg/s per m2 surrogate
+        pressure_factor = 0.6 + 0.4 * min(1.0, primary.pressure / max(0.1, constants.MAX_PRESSURE))
+        subcool = max(0.0, saturation_temperature(primary.pressure) - primary.temperature_out)
+        subcool_factor = max(0.2, min(1.0, subcool / 30.0))
+        flux_factor = max(0.4, min(3.0, (mass_flux / 500.0) ** 0.5))
+        base_chf = 1.2e7  # W/m2 surrogate baseline
+        return base_chf * flux_factor * pressure_factor * subcool_factor
 
     def _core_surface_area(self) -> float:
         # Simple surrogate: area scaling with volume^(2/3)