diff --git a/process/core/input.py b/process/core/input.py
index 31b6c6f6cd..bd490eead0 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -1186,6 +1186,7 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "f_len_sol_power_decay_inboard_outboard": InputVariable(
         "physics", float, range=(0.01, 2.0)
     ),
+    "molflow_plasma_fuelling": InputVariable("physics", float, range=(1.0e20, 1.0e23)),
 }
 
 
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 212f2f1e94..3a27cce530 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -2976,7 +2976,7 @@ def plot_main_plasma_information(
         f"   - Average mass of all plasma ions: {mfile.get('m_ions_total_amu', scan=scan):.3f} amu\n"
         f"Fuel mass: {mfile.get('m_plasma_fuel_ions', scan=scan) * 1000:.4f} g\n"
         f"   - Average mass of all fuel ions: {mfile.get('m_fuel_amu', scan=scan):.3f} amu\n\n"
-        f"Fueling rate: {mfile.get('molflow_plasma_fuelling_required', scan=scan):.3e} nucleus-pairs/s\n"
+        f"Fueling rate: {mfile.get('molflow_plasma_fuelling', scan=scan):.3e} nucleus-pairs/s\n"
         f"Fuel burn-up rate: {mfile.get('rndfuel', scan=scan):.3e} reactions/s \n"
         f"Burn-up fraction: {mfile.get('burnup', scan=scan):.4f} \n"
     )
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index c5d786d0f8..08792a461e 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -1981,6 +1981,20 @@ def constraint_equation_92(constraint_registration, data):
     )
 
 
+@ConstraintManager.register_constraint(93, "", "=")
+def constraint_equation_93(constraint_registration, data):
+    """Fuel-ion equilibrium constraint.
+
+    Requires the plasma fuelling rate to equal the fuelling rate required
+    for fuel-ion equilibrium.
+    """
+    return eq(
+        data.physics.molflow_plasma_fuelling,
+        data.physics.molflow_plasma_fuelling_equilibrium,
+        constraint_registration,
+    )
+
+
 def constraint_eqns(m: int, ieqn: int, data: DataStructure):
     """Evaluates the constraints given the current state of PROCESS.
 
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index fd468c5537..5a2dc0ee68 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -239,6 +239,12 @@ class IterationVariable:
     175: IterationVariable("kappa", "physics", 0.00, 10.00),
     176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30),
     177: IterationVariable("f_a_tf_turn_cable_space_extra_void", "tfcoil", 0.01, 1.0),
+    178: IterationVariable(
+        "molflow_plasma_fuelling",
+        "physics",
+        1.0e20,
+        1.0e23,
+    ),
 }
 
 
diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py
index 85785576a4..cc2c175b03 100644
--- a/process/data_structure/numerics.py
+++ b/process/data_structure/numerics.py
@@ -292,6 +292,7 @@ class NumericsData:
             "CS achievable stress load cycles lower limit           ",
             "ECRH ignitability                ",  # Stellarator constraint
             "Fuel composition consistency     ",
+            "Fuel ion equilibrium             ",
         ]
     )
     """Labels describing constraint equations (corresponding itvs)<UL>
@@ -394,6 +395,7 @@ class NumericsData:
     * (90) Lower Limit on number of stress load cycles for CS
     * (91) Checking if the design point is ECRH ignitable
     * (92) D/T/He3 ratio in fuel sums to 1
+    * (93) Fuel ion equilibrium
     """
 
     ixc: list[int] = field(
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index dad94478af..58597f4c6c 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1456,9 +1456,12 @@ class PhysicsData:
     """Plasma safety factor at 95% flux surface (q₉₅) (`iteration variable 18`)
     """
 
-    molflow_plasma_fuelling_required: float = 0.0
+    molflow_plasma_fuelling: float = 0.0
     """plasma fuelling rate (nucleus-pairs/s)"""
 
+    molflow_plasma_fuelling_equilibrium: float = 0.0
+    """Plasma fuelling rate required for fuel-ion equilibrium (nucleus-pairs/s)."""
+
     tauratio: float = 1.0
     """tauratio /1.0/ : ratio of He and pellet particle confinement times"""
 
diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py
index 208f98cbe0..df73ff27f0 100644
--- a/process/models/costs/costs.py
+++ b/process/models/costs/costs.py
@@ -2346,10 +2346,10 @@ def acc2272(self):
         This routine evaluates the Account 2272 - Fuel processing
         """
         if self.data.ife.ife != 1:
-            #  Previous calculation, using molflow_plasma_fuelling_required in Amps:
+            #  Previous calculation, using molflow_plasma_fuelling in Amps:
             #  1.3 should have been
             # self.data.physics.m_fuel_amu*umass/electron_charge*1000*s/day = 2.2
-            # wtgpd = burnup * molflow_plasma_fuelling_required * 1.3e0
+            # wtgpd = burnup * molflow_plasma_fuelling * 1.3e0
 
             #  New calculation: 2 nuclei * reactions/sec * kg/nucleus * g/kg * sec/day
             self.data.physics.wtgpd = (
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 43975e3f0e..17fd226706 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -946,7 +946,7 @@ def run(self):
             self.data.physics.burnup,
             self.data.physics.figmer,
             self.data.physics.fusrat,
-            self.data.physics.molflow_plasma_fuelling_required,
+            self.data.physics.molflow_plasma_fuelling_equilibrium,
             self.data.physics.rndfuel,
             self.data.physics.t_alpha_confinement,
             self.data.physics.f_t_alpha_energy_confinement,
@@ -1484,58 +1484,55 @@ def phyaux(
         vol_plasma: float,
         burnup_in: float,
         tauratio: float,
-    ) -> tuple[float, float, float, float, float, float, float, float]:
-        """Auxiliary physics quantities
+    ) -> tuple[float, float, float, float, float, float, float]:
+        """Calculate auxiliary plasma physics quantities.
 
         Parameters
         ----------
         aspect : float
             Plasma aspect ratio.
         nd_plasma_fuel_ions_vol_avg : float
-            Fuel ion density (/m3).
+            Volume-averaged fuel ion density (/m3).
         fusden_total_vol_avg : float
-            Fusion reaction rate from plasma and beams (/m3/s).
+            Volume-averaged fusion reaction rate density from plasma and beams (/m3/s).
         fusden_alpha_total_vol_avg : float
-            Alpha particle production rate (/m3/s).
+            Volume-averaged alpha particle production rate density (/m3/s).
         plasma_current : float
             Plasma current (A).
         sbar : float
             Exponent for aspect ratio (normally 1).
         nd_plasma_alphas_thermal_vol_avg : float
-            Alpha ash density (/m3).
+            Volume-averaged thermal alpha ash density (/m3).
         t_energy_confinement : float
             Global energy confinement time (s).
         vol_plasma : float
             Plasma volume (m3).
-        burnup_in: float
-            fractional plasma burnup user input
-        tauratio: float
-            ratio of He and pellet particle confinement times
+        burnup_in : float
+            Fractional plasma burnup user input.
+        tauratio : float
+            Ratio of helium ash to fuel particle confinement times.
 
         Returns
         -------
         tuple
             A tuple containing:
-            - burnup (float): Fractional plasma burnup.
-            - figmer (float): Physics figure of merit.
-            - fusrat (float): Number of fusion reactions per second.
-            - molflow_plasma_fuelling_required (float): Fuelling rate for D-T
-              (nucleus-pairs/sec).
-            - rndfuel (float): Fuel burnup rate (reactions/s).
-            - t_alpha_confinement (float): Alpha particle confinement time (s).
-            - f_t_alpha_energy_confinement (float): Fraction of alpha energy confinement.
-            This subroutine calculates extra physics related items needed by other
-            parts of the code.
 
+            - burnup: fractional plasma burnup.
+            - figmer: physics figure of merit.
+            - fusrat: total number of fusion reactions per second.
+            - molflow_plasma_fuelling_equilibrium: plasma fuelling rate required
+            for fuel-ion equilibrium (nucleus-pairs/s).
+            - rndfuel: fuel burnup rate (reactions/s).
+            - t_alpha_confinement: alpha particle confinement time (s).
+            - f_t_alpha_energy_confinement: fraction of alpha energy confinement.
         """
-        figmer = 1e-6 * plasma_current * aspect**sbar
+        figmer = 1.0e-6 * plasma_current * aspect**sbar
 
         # Fusion reactions per second
         fusrat = fusden_total_vol_avg * vol_plasma
 
         # Alpha particle confinement time (s)
         # Number of alphas / alpha production rate
-        # only likely if DD is only active fusion reaction
         t_alpha_confinement = (
             0.0
             if fusden_alpha_total_vol_avg == 0.0  # noqa: RUF069
@@ -1545,15 +1542,8 @@ def phyaux(
         # Fractional burnup
         # (Consider detailed model in: G. L. Jackson, V. S. Chan, R. D. Stambaugh,
         # Fusion Science and Technology, vol.64, no.1, July 2013, pp.8-12)
-        # The ratio of ash to fuel particle confinement times is given by
-        # tauratio
-        # Possible logic...
-        # burnup = fuel ion-pairs burned/m3 / initial fuel ion-pairs/m3;
-        # fuel ion-pairs burned/m3 = alpha particles/m3 (for both D-T and
-        # D-He3 reactions)
-        # initial fuel ion-pairs/m3 = burnt fuel ion-pairs/m3 + unburnt fuel-ion
-        # pairs/m3
-        # Remember that unburnt fuel-ion pairs/m3 = 0.5 * unburnt fuel-ions/m3
+        #
+        # burnup = fuel ion-pairs burned / initial fuel ion-pairs
         if burnup_in <= 1.0e-9:
             burnup = (
                 nd_plasma_alphas_thermal_vol_avg
@@ -1563,11 +1553,11 @@ def phyaux(
         else:
             burnup = burnup_in
 
-        # Fuel burnup rate (reactions/second) (previously Amps)
+        # Fuel burnup rate (reactions/s)
         rndfuel = fusrat
 
-        # Required fuelling rate (fuel ion pairs/second) (previously Amps)
-        molflow_plasma_fuelling_required = rndfuel / burnup
+        # Plasma fuelling rate required to satisfy fuel-ion equilibrium
+        molflow_plasma_fuelling_equilibrium = rndfuel / burnup
 
         f_t_alpha_energy_confinement = t_alpha_confinement / t_energy_confinement
 
@@ -1575,7 +1565,7 @@ def phyaux(
             burnup,
             figmer,
             fusrat,
-            molflow_plasma_fuelling_required,
+            molflow_plasma_fuelling_equilibrium,
             rndfuel,
             t_alpha_confinement,
             f_t_alpha_energy_confinement,
@@ -2435,8 +2425,15 @@ def outplas(self):
         po.ovarre(
             self.outfile,
             "Fuelling rate (nucleus-pairs/s)",
-            "(molflow_plasma_fuelling_required)",
-            self.data.physics.molflow_plasma_fuelling_required,
+            "(molflow_plasma_fuelling)",
+            self.data.physics.molflow_plasma_fuelling,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Fuelling rate required for equilibrium (nucleus-pairs/s)",
+            "(molflow_plasma_fuelling_equilibrium)",
+            self.data.physics.molflow_plasma_fuelling_equilibrium,
             "OP ",
         )
         po.ovarre(
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index d9b58d781c..37c98f051d 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2384,7 +2384,7 @@ def st_phys(self, output):
             self.data.physics.burnup,
             self.data.physics.figmer,
             _fusrat,
-            self.data.physics.molflow_plasma_fuelling_required,
+            self.data.physics.molflow_plasma_fuelling_equilibrium,
             self.data.physics.rndfuel,
             self.data.physics.t_alpha_confinement,
             self.data.physics.f_t_alpha_energy_confinement,
diff --git a/process/models/vacuum.py b/process/models/vacuum.py
index 4c55cd17cd..88b407f384 100644
--- a/process/models/vacuum.py
+++ b/process/models/vacuum.py
@@ -50,7 +50,7 @@ def run(self, output: bool = False):
         # MDK Check this!!
         gasld = (
             2.0e0
-            * self.data.physics.molflow_plasma_fuelling_required
+            * self.data.physics.molflow_plasma_fuelling
             * self.data.physics.m_fuel_amu
             * constants.UMASS
         )
@@ -115,7 +115,7 @@ def vacuum_simple(self, output) -> float:
         # One ITER torus cryopump has a throughput of 50 Pa m3/s = 1.2155e+22 molecules/s
         # Issue #304
         n_iter_vacuum_pumps = (
-            self.data.physics.molflow_plasma_fuelling_required
+            self.data.physics.molflow_plasma_fuelling
             / self.data.vacuum.molflow_vac_pumps
         )
 
@@ -164,8 +164,8 @@ def _vacuum_simple_output(self, n_iter_vacuum_pumps, npumpdown, npump):
         process_output.ovarre(
             self.outfile,
             "Plasma fuelling rate (nucleus-pairs/s)",
-            "(molflow_plasma_fuelling_required)",
-            self.data.physics.molflow_plasma_fuelling_required,
+            "(molflow_plasma_fuelling)",
+            self.data.physics.molflow_plasma_fuelling,
             "OP ",
         )
 
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index 516404b354..b8d0309942 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -1873,7 +1873,7 @@ class PhyauxParam(NamedTuple):
 
     expected_fusrat: Any = None
 
-    expected_molflow_plasma_fuelling_required: Any = None
+    expected_molflow_plasma_fuelling_equilibrium: Any = None
 
     expected_rndfuel: Any = None
 
@@ -1898,7 +1898,7 @@ class PhyauxParam(NamedTuple):
             expected_burnup=0.20383508579699033,
             expected_figmer=55.195367036602576,
             expected_fusrat=3.7484146722826997e20,
-            expected_molflow_plasma_fuelling_required=1.838944781084418e21,
+            expected_molflow_plasma_fuelling_equilibrium=1.838944781084418e21,
             expected_rndfuel=3.7484146722826997e20,
             expected_t_alpha_confinement=37.993985551650177,
         ),
@@ -1917,7 +1917,7 @@ class PhyauxParam(NamedTuple):
             expected_burnup=0.20387039462081086,
             expected_figmer=55.195367036602576,
             expected_fusrat=3.7467489360461772e20,
-            expected_molflow_plasma_fuelling_required=1.8378092331723546e21,
+            expected_molflow_plasma_fuelling_equilibrium=1.8378092331723546e21,
             expected_rndfuel=3.7467489360461772e20,
             expected_t_alpha_confinement=38.010876984618747,
         ),
@@ -1945,7 +1945,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics):
         burnup,
         figmer,
         fusrat,
-        molflow_plasma_fuelling_required,
+        molflow_plasma_fuelling_equilibrium,
         rndfuel,
         t_alpha_confinement,
         _,
@@ -1969,8 +1969,8 @@ def test_phyaux(phyauxparam, monkeypatch, physics):
 
     assert fusrat == pytest.approx(phyauxparam.expected_fusrat)
 
-    assert molflow_plasma_fuelling_required == pytest.approx(
-        phyauxparam.expected_molflow_plasma_fuelling_required
+    assert molflow_plasma_fuelling_equilibrium == pytest.approx(
+        phyauxparam.expected_molflow_plasma_fuelling_equilibrium
     )
 
     assert rndfuel == pytest.approx(phyauxparam.expected_rndfuel)
diff --git a/tests/unit/models/test_vacuum.py b/tests/unit/models/test_vacuum.py
index f1913200a2..60d59cd492 100644
--- a/tests/unit/models/test_vacuum.py
+++ b/tests/unit/models/test_vacuum.py
@@ -40,7 +40,7 @@ def test_simple_model(monkeypatch, vacuum):
         """
         monkeypatch.setattr(
             vacuum.data.physics,
-            "molflow_plasma_fuelling_required",
+            "molflow_plasma_fuelling",
             7.5745668997694112e22,
         )
         monkeypatch.setattr(vacuum.data.physics, "a_plasma_surface", 1500.3146527709359)