diff --git a/documentation/source/physics-models/plasma_composition.md b/documentation/source/physics-models/plasma_composition.md
index 32f6a1b78e..910a66e359 100644
--- a/documentation/source/physics-models/plasma_composition.md
+++ b/documentation/source/physics-models/plasma_composition.md
@@ -36,7 +36,7 @@ using input array `f_nd_impurity_electrons(1,...,14)`. The available species alo
 As stated above, the number density fractions for hydrogen (all isotopes) and
 helium should not be set, as they are calculated by the code. This is to ensure 
 plasma quasi-neutrality taking into account the fuel ratios
-`f_plasma_fuel_deuterium`, `f_plasma_fuel_tritium` and `f_plasma_fuel_helium3`, and the alpha particle fraction `f_nd_alpha_electron` which may 
+`f_plasma_fuel_deuterium`, `f_plasma_fuel_tritium` and `f_plasma_fuel_helium3`, and the alpha particle fraction `f_nd_alpha_thermal_electron` which may 
 be input by the user or selected as an iteration variable.
 
 !!! note "Location of impurities"
@@ -59,8 +59,8 @@ All of the plasma composites are normally given as a fraction of the volume aver
 
 1. **Alpha Ash Portion Calculation**
 
-    - Calculate the number density of alpha particles (`nd_plasma_alphas_vol_avg`) using the electron density (`nd_plasma_electrons_vol_avg`) and the alpha particle to electron ratio (`f_nd_alpha_electron`).
-        - `f_nd_alpha_electron` can be set as an iteration variable (`ixc = 109`) or set directly.
+    - Calculate the number density of thermal alpha particles (`nd_plasma_alphas_thermal_vol_avg`) using the electron density (`nd_plasma_electrons_vol_avg`) and the alpha particle to electron ratio (`f_nd_alpha_thermal_electron`).
+        - `f_nd_alpha_thermal_electron` can be set as an iteration variable (`ixc = 109`) or set directly.
 
     $$
     n_{\alpha} = \mathtt{f\_nd\_alpha\_electron}\times n_{\text{e}}
diff --git a/examples/data/large_tokamak_eval_IN.DAT b/examples/data/large_tokamak_eval_IN.DAT
index 028de03d25..28a7240878 100644
--- a/examples/data/large_tokamak_eval_IN.DAT
+++ b/examples/data/large_tokamak_eval_IN.DAT
@@ -301,7 +301,7 @@ i_plasma_geometry   = 0 * switch for plasma cross-sectional shape calculation;
 kappa    = 1.85 * plasma separatrix elongation (calculated if `i_plasma_geometry = 1-5; 7 or 9-10`)
 q95        = 3.7339078193128556 * Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95
 q0       = 1.0 * safety factor on axis
-f_nd_alpha_electron   = 0.060238763988650204 * thermal alpha density/electron density (`iteration variable 109`)
+f_nd_alpha_thermal_electron   = 0.060238763988650204 * thermal alpha density/electron density (`iteration variable 109`)
 rmajor   = 8.0 * plasma major radius (m) (`iteration variable 3`)
 i_single_null = 1 * switch for single null / double null plasma;
 f_sync_reflect    = 0.6 * synchrotron wall reflectivity factor
diff --git a/examples/data/large_tokamak_varied_min_net_electric_IN.DAT b/examples/data/large_tokamak_varied_min_net_electric_IN.DAT
index 1d40f7e9b7..be0f9d527a 100644
--- a/examples/data/large_tokamak_varied_min_net_electric_IN.DAT
+++ b/examples/data/large_tokamak_varied_min_net_electric_IN.DAT
@@ -230,7 +230,7 @@ boundl(60) = 65000.0
 boundu(60) = 90000.0
 c_tf_turn = 65000.0
 
-* Helium fraction (f_nd_alpha_electron)
+* Helium fraction (f_nd_alpha_thermal_electron)
 ixc = 109
 boundu(109) = 0.1
 
diff --git a/examples/data/large_tokamak_varyrun_IN.DAT b/examples/data/large_tokamak_varyrun_IN.DAT
index f95574a5e1..f65eca9829 100644
--- a/examples/data/large_tokamak_varyrun_IN.DAT
+++ b/examples/data/large_tokamak_varyrun_IN.DAT
@@ -230,7 +230,7 @@ boundl(60) = 65000.0
 boundu(60) = 90000.0
 c_tf_turn = 65000.0
 
-* Helium fraction (f_nd_alpha_electron)
+* Helium fraction (f_nd_alpha_thermal_electron)
 ixc = 109
 boundu(109) = 0.1
 
diff --git a/process/core/init.py b/process/core/init.py
index 73fe6a7b72..1f0a049c4c 100644
--- a/process/core/init.py
+++ b/process/core/init.py
@@ -334,7 +334,7 @@ def check_process(inputs, data):  # noqa: ARG001
     if data.impurity_radiation.f_nd_impurity_electrons[1] != 0.1:  # noqa: RUF069
         raise ProcessValidationError(
             "The thermal alpha/electron density ratio should be controlled using"
-            " f_nd_alpha_electron (itv 109) and not f_nd_impurity_electrons(2)."
+            " f_nd_alpha_thermal_electron (itv 109) and not f_nd_impurity_electrons(2)."
             "f_nd_impurity_electrons(2) should be removed from the input file,"
             " or set to the default value 0.1D0."
         )
diff --git a/process/core/input.py b/process/core/input.py
index 80c35712a6..09138f9072 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -414,7 +414,7 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "f_asym": InputVariable("stellarator", float, range=(0.9, 2.0)),
     "f_fw_peak": InputVariable("fwbs", float, range=(1.0, 100.0)),
     "f_fw_rad_max": InputVariable("constraints", float, range=(0.1, 10)),
-    "f_nd_alpha_electron": InputVariable("physics", float, range=(1e-12, 1.0)),
+    "f_nd_alpha_thermal_electron": InputVariable("physics", float, range=(1e-12, 1.0)),
     "f_nd_beam_electron": InputVariable("physics", float, range=(0.0, 1.0)),
     "f_nd_protium_electrons": InputVariable("physics", float, range=(0.0, 1.0)),
     "f_neut_shield": InputVariable("fwbs", float, range=(0.0, 1.0)),
diff --git a/process/core/io/mfile/comparison.py b/process/core/io/mfile/comparison.py
index 95825a1d15..cd0afce7da 100644
--- a/process/core/io/mfile/comparison.py
+++ b/process/core/io/mfile/comparison.py
@@ -72,7 +72,7 @@
     "p_coolant_pump_elec_total_mw",
     "p_plant_electric_net_mw",
     "pflux_fw_neutron_mw",
-    "f_nd_alpha_electron",
+    "f_nd_alpha_thermal_electron",
     "p_plasma_inner_rad_mw",
     "p_plasma_rad_mw",
     "p_blkt_nuclear_heat_total_mw",
@@ -127,7 +127,7 @@
     "t_energy_confinement",
     "hfact",
     "tauelaw",
-    "f_nd_alpha_electron",
+    "f_nd_alpha_thermal_electron",
     "pflux_fw_neutron_mw",
     "p_plasma_inner_rad_mw",
     "p_plasma_sync_mw",
@@ -221,7 +221,7 @@
     "p_plant_electric_net_mw",
     "pflux_fw_neutron_mw",
     "pflux_fw_rad_mw",
-    "f_nd_alpha_electron",
+    "f_nd_alpha_thermal_electron",
     "p_plasma_inner_rad_mw",
     "p_plasma_rad_mw",
     "f_c_plasma_bootstrap",
diff --git a/process/core/io/obsolete_vars.py b/process/core/io/obsolete_vars.py
index 0e661e845b..bce1b7c86d 100644
--- a/process/core/io/obsolete_vars.py
+++ b/process/core/io/obsolete_vars.py
@@ -155,7 +155,7 @@
     "nimp": "n_impurities",
     "ssync": "f_sync_reflect",
     "rnbeam": "f_nd_beam_electron",
-    "ralpne": "f_nd_alpha_electron",
+    "ralpne": "f_nd_alpha_thermal_electron",
     "protium": "f_nd_protium_electrons",
     "clhsf": "f_z_cryostat",
     "ddwex": "dr_cryostat",
@@ -470,6 +470,7 @@
     "fiooic": "f_j_tf_wp_critical_max",
     "alstroh": "stress_cs_steel_max",
     "i_cs_stress": None,
+    "f_nd_alpha_electron": "f_nd_alpha_thermal_electron",
 }
 
 OBS_VARS_HELP = {
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index fdffe8a73d..eee3961f52 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -3876,7 +3876,7 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
     scan: int :
 
     """
-    nd_alphas = mfile.get("nd_plasma_alphas_vol_avg", scan=scan)
+    nd_alphas = mfile.get("nd_plasma_alphas_thermal_vol_avg", scan=scan)
     nd_protons = mfile.get("nd_plasma_protons_vol_avg", scan=scan)
     nd_impurities = mfile.get("nd_plasma_impurities_vol_avg", scan=scan)
     nd_ions_total = mfile.get("nd_plasma_ions_total_vol_avg", scan=scan)
@@ -3973,7 +3973,7 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
     ax_main.plot(
         rho,
         density_profiles_plotting[1],
-        label=r"$n_{\alpha}$",
+        label=r"$n_{\alpha,\text{thermal}}$",
         color="#d62728",
         linewidth=1.5,
     )
@@ -4125,8 +4125,8 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
             f"{mfile.get('nd_plasma_fuel_ions_vol_avg', scan=scan):.3e} m$^{{-3}}$"
         ),
         (
-            r"$\langle n_{\text{alpha}} \rangle $: "
-            f"{mfile.get('nd_plasma_alphas_vol_avg', scan=scan):.3e} m$^{{-3}}$"
+            r"$\langle n_{\alpha,\text{thermal}} \rangle $: "
+            f"{mfile.get('nd_plasma_alphas_thermal_vol_avg', scan=scan):.3e} m$^{{-3}}$"
         ),
         (
             r"$\langle n_{\text{impurities}} \rangle $: "
diff --git a/process/core/io/variable_metadata.py b/process/core/io/variable_metadata.py
index 247dc38eb9..91ccf19fc8 100644
--- a/process/core/io/variable_metadata.py
+++ b/process/core/io/variable_metadata.py
@@ -44,7 +44,7 @@ class VariableMetadata:
         description="Effective energy confinement time",
         units="s",
     ),
-    "f_nd_alpha_electron": VariableMetadata(
+    "f_nd_alpha_thermal_electron": VariableMetadata(
         latex=r"$f_\mathrm{\alpha}$",
         description="Alpha particle heating fraction",
         units="",
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index 8e0de61362..490022a270 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -92,7 +92,7 @@ class IterationVariable:
     98: IterationVariable("f_blkt_li6_enrichment", "fwbs", 10.00, 100.0),
     104: IterationVariable("fcwr", "constraints", 0.001, 1.0),
     108: IterationVariable("breeder_f", "fwbs", 0.060, 1.0),
-    109: IterationVariable("f_nd_alpha_electron", "physics", 0.05, 0.15),
+    109: IterationVariable("f_nd_alpha_thermal_electron", "physics", 0.05, 0.15),
     114: IterationVariable("len_fw_channel", "fwbs", 0.001, 1.0e3),
     119: IterationVariable("temp_plasma_separatrix_kev", "physics", 0.0, 1.0e1),
     122: IterationVariable("f_a_cs_turn_steel", "pf_coil", 0.001, 0.950),
diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py
index 837708a1ae..8daf03fb78 100644
--- a/process/data_structure/numerics.py
+++ b/process/data_structure/numerics.py
@@ -471,7 +471,7 @@ class NumericsData:
     * (106) NOT USED
     * (107) NOT USED
     * (108) breeder_f: Volume of Li4SiO4 / (Volume of Be12Ti + Li4SiO4)
-    * (109) f_nd_alpha_electron: thermal alpha density / electron density
+    * (109) f_nd_alpha_thermal_electron: thermal alpha density / electron density
     * (110) NOT USED
     * (111) NOT USED
     * (112) NOT USED
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 0cafda166e..7477961576 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -609,8 +609,8 @@ class PhysicsData:
     )
     """Array of plasma electron density upper limits values (nₑ,max) [/m³]"""
 
-    nd_plasma_alphas_vol_avg: float = 0.0
-    """Plasma volume averaged thermal alpha density (⟨n_α⟩) [/m³]"""  # noqa: RUF001
+    nd_plasma_alphas_thermal_vol_avg: float = 0.0
+    """Plasma volume averaged thermal alpha density (⟨n_αₜₕ⟩) [/m³]"""
 
     nd_beam_ions: float = 0.0
     """Hot beam ion density, variable (⟨n_beam⟩) [/m³]"""
@@ -1258,8 +1258,8 @@ class PhysicsData:
     rad_fraction_total: float = 0.0
     """Radiation fraction total = SoL + LCFS radiation / total power deposited in plasma"""
 
-    f_nd_alpha_electron: float = 0.1
-    """thermal alpha density/electron density (`iteration variable 109`)"""
+    f_nd_alpha_thermal_electron: float = 0.1
+    """Thermal alpha density/electron density (⟨n_αₜₕ⟩/⟨nₑ⟩)"""
 
     f_nd_protium_electrons: float = 0.0
     """Seeded f_nd_protium_electrons density / electron density."""
diff --git a/process/models/physics/current_drive.py b/process/models/physics/current_drive.py
index 9d63db895a..811f6c38bd 100644
--- a/process/models/physics/current_drive.py
+++ b/process/models/physics/current_drive.py
@@ -191,7 +191,7 @@ def iternb(self):
             self.data.current_drive.e_beam_kev / self.data.physics.m_beam_amu,
             self.data.physics.temp_plasma_electron_vol_avg_kev,
             self.data.physics.nd_plasma_electrons_vol_avg,
-            self.data.physics.f_nd_alpha_electron,
+            self.data.physics.f_nd_alpha_thermal_electron,
             self.data.physics.f_nd_plasma_carbon_electron,
             self.data.physics.f_nd_plasma_oxygen_electron,
             self.data.physics.f_nd_plasma_iron_argon_electron,
@@ -283,7 +283,7 @@ def culnbi(self) -> tuple[float, float, float]:
             self.data.current_drive.e_beam_kev / self.data.physics.m_beam_amu,
             self.data.physics.temp_plasma_electron_vol_avg_kev,
             self.data.physics.nd_plasma_electrons_vol_avg,
-            self.data.physics.f_nd_alpha_electron,
+            self.data.physics.f_nd_alpha_thermal_electron,
             self.data.physics.f_nd_plasma_carbon_electron,
             self.data.physics.f_nd_plasma_oxygen_electron,
             self.data.physics.f_nd_plasma_iron_argon_electron,
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 795dc0899a..012dd0b307 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -256,7 +256,7 @@ def run(self):
             self.data.physics.m_ions_total_amu,
             self.data.physics.nd_plasma_ions_total_vol_avg,
             self.data.physics.nd_plasma_fuel_ions_vol_avg,
-            self.data.physics.nd_plasma_alphas_vol_avg,
+            self.data.physics.nd_plasma_alphas_thermal_vol_avg,
             self.data.physics.vol_plasma,
             self.data.physics.nd_plasma_electrons_vol_avg,
         )
@@ -944,7 +944,7 @@ def run(self):
             self.data.physics.fusden_alpha_total,
             self.data.physics.plasma_current,
             sbar,
-            self.data.physics.nd_plasma_alphas_vol_avg,
+            self.data.physics.nd_plasma_alphas_thermal_vol_avg,
             self.data.physics.t_energy_confinement,
             self.data.physics.vol_plasma,
             self.data.physics.burnup_in,
@@ -1155,9 +1155,9 @@ def plasma_composition(self):
           charge.
         """
         # Alpha ash portion
-        self.data.physics.nd_plasma_alphas_vol_avg = (
+        self.data.physics.nd_plasma_alphas_thermal_vol_avg = (
             self.data.physics.nd_plasma_electrons_vol_avg
-            * self.data.physics.f_nd_alpha_electron
+            * self.data.physics.f_nd_alpha_thermal_electron
         )
 
         # ======================================================================
@@ -1172,14 +1172,14 @@ def plasma_composition(self):
             self.data.physics.nd_plasma_protons_vol_avg = max(
                 self.data.physics.f_nd_protium_electrons
                 * self.data.physics.nd_plasma_electrons_vol_avg,
-                self.data.physics.nd_plasma_alphas_vol_avg
+                self.data.physics.nd_plasma_alphas_thermal_vol_avg
                 * (self.data.physics.f_plasma_fuel_helium3 + 1.0e-3),
             )  # rough estimate
         else:
             self.data.physics.nd_plasma_protons_vol_avg = max(
                 self.data.physics.f_nd_protium_electrons
                 * self.data.physics.nd_plasma_electrons_vol_avg,
-                self.data.physics.nd_plasma_alphas_vol_avg
+                self.data.physics.nd_plasma_alphas_thermal_vol_avg
                 * self.data.physics.proton_rate_density
                 / self.data.physics.fusden_alpha_total,
             )
@@ -1220,7 +1220,7 @@ def plasma_composition(self):
         # znfuel is the sum of Zi.ni for the three fuel ions
         znfuel = (
             self.data.physics.nd_plasma_electrons_vol_avg
-            - 2.0 * self.data.physics.nd_plasma_alphas_vol_avg
+            - 2.0 * self.data.physics.nd_plasma_alphas_thermal_vol_avg
             - self.data.physics.nd_plasma_protons_vol_avg
             - self.data.physics.nd_beam_ions
             - znimp
@@ -1259,7 +1259,7 @@ def plasma_composition(self):
             self.data.physics.f_plasma_fuel_helium3
             * self.data.physics.nd_plasma_fuel_ions_vol_avg
             / self.data.physics.nd_plasma_electrons_vol_avg
-            + self.data.physics.f_nd_alpha_electron
+            + self.data.physics.f_nd_alpha_thermal_electron
         )
 
         # ======================================================================
@@ -1278,7 +1278,7 @@ def plasma_composition(self):
         # Total ion density
         self.data.physics.nd_plasma_ions_total_vol_avg = (
             self.data.physics.nd_plasma_fuel_ions_vol_avg
-            + self.data.physics.nd_plasma_alphas_vol_avg
+            + self.data.physics.nd_plasma_alphas_thermal_vol_avg
             + self.data.physics.nd_plasma_protons_vol_avg
             + self.data.physics.nd_beam_ions
             + self.data.physics.nd_plasma_impurities_vol_avg
@@ -1374,7 +1374,10 @@ def plasma_composition(self):
                 self.data.physics.m_fuel_amu
                 * self.data.physics.nd_plasma_fuel_ions_vol_avg
             )
-            + (constants.M_ALPHA_AMU * self.data.physics.nd_plasma_alphas_vol_avg)
+            + (
+                constants.M_ALPHA_AMU
+                * self.data.physics.nd_plasma_alphas_thermal_vol_avg
+            )
             + (self.data.physics.nd_plasma_protons_vol_avg * constants.M_PROTON_AMU)
             + (self.data.physics.m_beam_amu * self.data.physics.nd_beam_ions)
         )
@@ -1411,7 +1414,11 @@ def plasma_composition(self):
                 * self.data.physics.nd_plasma_fuel_ions_vol_avg
                 / constants.M_HELION_AMU
             )
-            + (4.0 * self.data.physics.nd_plasma_alphas_vol_avg / constants.M_ALPHA_AMU)
+            + (
+                4.0
+                * self.data.physics.nd_plasma_alphas_thermal_vol_avg
+                / constants.M_ALPHA_AMU
+            )
             + (self.data.physics.nd_plasma_protons_vol_avg / constants.M_PROTON_AMU)
             + (
                 (1.0 - self.data.current_drive.f_beam_tritium)
@@ -1448,7 +1455,7 @@ def phyaux(
         fusden_alpha_total: float,
         plasma_current: float,
         sbar: float,
-        nd_plasma_alphas_vol_avg: float,
+        nd_plasma_alphas_thermal_vol_avg: float,
         t_energy_confinement: float,
         vol_plasma: float,
         burnup_in: float,
@@ -1470,7 +1477,7 @@ def phyaux(
             Plasma current (A).
         sbar : float
             Exponent for aspect ratio (normally 1).
-        nd_plasma_alphas_vol_avg : float
+        nd_plasma_alphas_thermal_vol_avg : float
             Alpha ash density (/m3).
         t_energy_confinement : float
             Global energy confinement time (s).
@@ -1508,7 +1515,7 @@ def phyaux(
         t_alpha_confinement = (
             0.0
             if fusden_alpha_total == 0.0  # noqa: RUF069
-            else nd_plasma_alphas_vol_avg / fusden_alpha_total
+            else nd_plasma_alphas_thermal_vol_avg / fusden_alpha_total
         )
 
         # Fractional burnup
@@ -1525,8 +1532,8 @@ def phyaux(
         # Remember that unburnt fuel-ion pairs/m3 = 0.5 * unburnt fuel-ions/m3
         if burnup_in <= 1.0e-9:
             burnup = (
-                nd_plasma_alphas_vol_avg
-                / (nd_plasma_alphas_vol_avg + 0.5 * nd_plasma_fuel_ions_vol_avg)
+                nd_plasma_alphas_thermal_vol_avg
+                / (nd_plasma_alphas_thermal_vol_avg + 0.5 * nd_plasma_fuel_ions_vol_avg)
                 / tauratio
             )
         else:
@@ -2702,16 +2709,16 @@ def output_temperature_density_profile_info(self) -> None:
         )
         po.ovarre(
             self.outfile,
-            "Thermalised alpha volume averaged number density (⟨n_alpha⟩) (/m³)",
-            "(nd_plasma_alphas_vol_avg)",
-            self.data.physics.nd_plasma_alphas_vol_avg,
+            "Thermalised alpha volume averaged number density (⟨n_αₜₕ⟩) (/m³)",
+            "(nd_plasma_alphas_thermal_vol_avg)",
+            self.data.physics.nd_plasma_alphas_thermal_vol_avg,
             "OP ",
         )
         po.ovarre(
             self.outfile,
-            "Thermalised alpha to electron number density ratio (⟨n_alpha⟩/⟨nₑ⟩)",
-            "(f_nd_alpha_electron)",
-            self.data.physics.f_nd_alpha_electron,
+            "Thermalised alpha to electron number density ratio (⟨n_αₜₕ⟩/⟨nₑ⟩)",
+            "(f_nd_alpha_thermal_electron)",
+            self.data.physics.f_nd_alpha_thermal_electron,
         )
         po.ovarre(
             self.outfile,
@@ -3021,7 +3028,7 @@ def calculate_plasma_masses(
         m_ions_total_amu: float,
         nd_plasma_ions_total_vol_avg: float,
         nd_plasma_fuel_ions_vol_avg: float,
-        nd_plasma_alphas_vol_avg: float,
+        nd_plasma_alphas_thermal_vol_avg: float,
         vol_plasma: float,
         nd_plasma_electrons_vol_avg: float,
     ) -> tuple[float, float, float, float, float]:
@@ -3037,7 +3044,7 @@ def calculate_plasma_masses(
             Total ion density (/m3).
         nd_plasma_fuel_ions_vol_avg : float
             Fuel ion density (/m3).
-        nd_plasma_alphas_vol_avg : float
+        nd_plasma_alphas_thermal_vol_avg : float
             Alpha ash density (/m3).
         vol_plasma : float
             Plasma volume (m3).
@@ -3059,7 +3066,9 @@ def calculate_plasma_masses(
             nd_plasma_ions_total_vol_avg * vol_plasma
         )
 
-        m_plasma_alpha = (nd_plasma_alphas_vol_avg * vol_plasma) * constants.ALPHA_MASS
+        m_plasma_alpha = (
+            nd_plasma_alphas_thermal_vol_avg * vol_plasma
+        ) * constants.ALPHA_MASS
 
         m_plasma_electron = constants.ELECTRON_MASS * (
             nd_plasma_electrons_vol_avg * vol_plasma
@@ -5332,7 +5341,7 @@ def run(self):
             nd_plasma_ions=(
                 self.data.physics.nd_plasma_electrons_vol_avg
                 * (
-                    self.data.physics.nd_plasma_alphas_vol_avg
+                    self.data.physics.nd_plasma_alphas_thermal_vol_avg
                     / self.data.physics.nd_plasma_electrons_vol_avg
                 )
             ),
@@ -5345,7 +5354,7 @@ def run(self):
             nd_plasma_ions=(
                 self.plasma_profile.neprofile.profile_y
                 * (
-                    self.data.physics.nd_plasma_alphas_vol_avg
+                    self.data.physics.nd_plasma_alphas_thermal_vol_avg
                     / self.data.physics.nd_plasma_electrons_vol_avg
                 )
             ),
diff --git a/process/models/stellarator/neoclassics.py b/process/models/stellarator/neoclassics.py
index f555675c7a..9eec4689f1 100644
--- a/process/models/stellarator/neoclassics.py
+++ b/process/models/stellarator/neoclassics.py
@@ -180,7 +180,7 @@ def init_profile_values_from_PROCESS(self, rho):
             * (1 - rho**2) ** self.data.physics.alphan
         )
         densa = (
-            self.data.physics.nd_plasma_alphas_vol_avg
+            self.data.physics.nd_plasma_alphas_thermal_vol_avg
             * (1 + self.data.physics.alphan)
             * (1 - rho**2) ** self.data.physics.alphan
         )
@@ -261,7 +261,7 @@ def init_profile_values_from_PROCESS(self, rho):
             * 1.0
             / self.data.physics.rminor
             * rho
-            * self.data.physics.nd_plasma_alphas_vol_avg
+            * self.data.physics.nd_plasma_alphas_thermal_vol_avg
             * (1 + self.data.physics.alphan)
             * (1.0 - rho**2) ** (self.data.physics.alphan - 1.0)
             * self.data.physics.alphan
@@ -550,7 +550,7 @@ def neoclassics_calc_nu_star_fromT(self, iota):
             * self.data.physics.f_plasma_fuel_deuterium,
             self.data.physics.nd_plasma_fuel_ions_vol_avg
             * (1 - self.data.physics.f_plasma_fuel_deuterium),
-            self.data.physics.nd_plasma_alphas_vol_avg,
+            self.data.physics.nd_plasma_alphas_thermal_vol_avg,
         ])
 
         mass = np.array([
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 3482874866..326d41aa3e 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2358,7 +2358,7 @@ def st_phys(self, output):
             self.data.physics.fusden_alpha_total,
             self.data.physics.plasma_current,
             sbar,
-            self.data.physics.nd_plasma_alphas_vol_avg,
+            self.data.physics.nd_plasma_alphas_thermal_vol_avg,
             self.data.physics.t_energy_confinement,
             self.data.physics.vol_plasma,
             self.data.physics.burnup_in,
diff --git a/tests/integration/data/large_tokamak_eval.IN.DAT b/tests/integration/data/large_tokamak_eval.IN.DAT
index aaef0417ce..91cb0d70be 100644
--- a/tests/integration/data/large_tokamak_eval.IN.DAT
+++ b/tests/integration/data/large_tokamak_eval.IN.DAT
@@ -301,7 +301,7 @@ i_plasma_geometry   = 0 * switch for plasma cross-sectional shape calculation;
 kappa    = 1.85 * plasma separatrix elongation (calculated if `i_plasma_geometry = 1-5; 7 or 9-10`)
 q95        = 3.7339078193128556 * Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95
 q0       = 1.0 * safety factor on axis
-f_nd_alpha_electron   = 0.060238763988650204 * thermal alpha density/electron density (`iteration variable 109`)
+f_nd_alpha_thermal_electron   = 0.060238763988650204 * thermal alpha density/electron density (`iteration variable 109`)
 rmajor   = 8.0 * plasma major radius (m) (`iteration variable 3`)
 i_single_null = 1 * switch for single null / double null plasma;
 f_sync_reflect    = 0.6 * synchrotron wall reflectivity factor
diff --git a/tests/regression/input_files/large_tokamak_eval.IN.DAT b/tests/regression/input_files/large_tokamak_eval.IN.DAT
index 8a4b55ca1e..79df0d3a10 100644
--- a/tests/regression/input_files/large_tokamak_eval.IN.DAT
+++ b/tests/regression/input_files/large_tokamak_eval.IN.DAT
@@ -302,7 +302,7 @@ i_plasma_geometry   = 0 * switch for plasma cross-sectional shape calculation;
 kappa    = 1.85 * plasma separatrix elongation (calculated if `i_plasma_geometry = 1-5; 7 or 9-10`)
 q95        = 3.7339078193128556 * Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95
 q0       = 1.0 * safety factor on axis
-f_nd_alpha_electron   = 0.060238763988650204 * thermal alpha density/electron density (`iteration variable 109`)
+f_nd_alpha_thermal_electron   = 0.060238763988650204 * thermal alpha density/electron density (`iteration variable 109`)
 rmajor   = 8.0 * plasma major radius (m) (`iteration variable 3`)
 i_single_null = 1 * switch for single null / double null plasma;
 f_sync_reflect    = 0.6 * synchrotron wall reflectivity factor
diff --git a/tests/regression/input_files/large_tokamak_nof.IN.DAT b/tests/regression/input_files/large_tokamak_nof.IN.DAT
index 2e75a17a8b..fc97d410d4 100644
--- a/tests/regression/input_files/large_tokamak_nof.IN.DAT
+++ b/tests/regression/input_files/large_tokamak_nof.IN.DAT
@@ -241,7 +241,7 @@ boundl(60) = 65000.0
 boundu(60) = 90000.0
 c_tf_turn = 65000.0
 
-* Helium fraction (f_nd_alpha_electron)
+* Helium fraction (f_nd_alpha_thermal_electron)
 ixc = 109
 boundu(109) = 0.1
 
diff --git a/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT b/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT
index 5af73487c4..ff8bf4e594 100644
--- a/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT
+++ b/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT
@@ -286,7 +286,7 @@ f_h_mode_margin = 1.1
 * JUSTIFICATION: Used with icc=15 to enforce H-mode
 
 ixc = 109
-f_nd_alpha_electron = 6.14955494147856815e-02
+f_nd_alpha_thermal_electron = 6.14955494147856815e-02
 * DESCRIPTION:   Thermal alpha density / electron density
 * JUSTIFICATION: Calculate He fraction consistent with tau ratio
 
diff --git a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
index 22d59307f8..7ea0b50584 100644
--- a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
+++ b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
@@ -280,7 +280,7 @@ itartpf  = 1 * switch for Spherical Tokamak PF models;
 kappa    = 2.8 * plasma separatrix elongation (calculated if `i_plasma_geometry = 1-5; 7 or 9-10`)
 q95        = 5.835830999686161 * Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95
 q0       = 2.0 * safety factor on axis
-f_nd_alpha_electron = 0.08870796537675113 * thermal alpha density/electron density (`iteration variable 109`)
+f_nd_alpha_thermal_electron = 0.08870796537675113 * thermal alpha density/electron density (`iteration variable 109`)
 ind_plasma_internal_norm      = 0.3 * plasma normalised internal inductance
 i_ind_plasma_internal_norm = 0 * Normalised plasma intenral induction selection switch (0 = user input) 
 rmajor   = 4.5 * plasma major radius (m) (`iteration variable 3`)
diff --git a/tests/regression/input_files/st_regression.IN.DAT b/tests/regression/input_files/st_regression.IN.DAT
index 0122b7f5e5..165335d9b8 100644
--- a/tests/regression/input_files/st_regression.IN.DAT
+++ b/tests/regression/input_files/st_regression.IN.DAT
@@ -252,7 +252,7 @@ boundu(6) = 5.0E+20
 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 ixc = 109
-f_nd_alpha_electron = 0.2
+f_nd_alpha_thermal_electron = 0.2
 boundl(109) = 0.001
 boundu(109) = 1.0
 * DESCRIPTION:   Thermal alpha density / electron density
diff --git a/tests/regression/input_files/stellarator_helias.IN.DAT b/tests/regression/input_files/stellarator_helias.IN.DAT
index db6a453e34..90133d84fe 100644
--- a/tests/regression/input_files/stellarator_helias.IN.DAT
+++ b/tests/regression/input_files/stellarator_helias.IN.DAT
@@ -100,7 +100,7 @@ boundu(10) = 1.2
 * ixc = 50                * itv_fiooic
 f_j_tf_wp_critical_max = 0.8
 
-ixc = 109 * f_nd_alpha_electron: thermal alpha density / electron density
+ixc = 109 * f_nd_alpha_thermal_electron: thermal alpha density / electron density
 boundl(109) = 0.0001
 boundu(109) = 0.4
 
diff --git a/tests/unit/core/test_mfile.py b/tests/unit/core/test_mfile.py
index ed5961a7ac..019317d0ef 100644
--- a/tests/unit/core/test_mfile.py
+++ b/tests/unit/core/test_mfile.py
@@ -32,7 +32,7 @@ def test_get_mfile_initial_ixc_values(input_file, tmp_path, data_structure_obj):
     assert iteration_variable_values[-1] == pytest.approx(0.5)
 
     # A default not provided in the MFile
-    assert iteration_variable_names[-4] == "f_nd_alpha_electron"
+    assert iteration_variable_names[-4] == "f_nd_alpha_thermal_electron"
     assert iteration_variable_values[-4] == pytest.approx(0.1)
 
 
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index 0c6b09ecf4..897258deff 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -1253,7 +1253,7 @@ class PlasmaCompositionParam(NamedTuple):
 
     f_alpha_ion: Any = None
 
-    f_nd_alpha_electron: Any = None
+    f_nd_alpha_thermal_electron: Any = None
 
     f_nd_beam_electron: Any = None
 
@@ -1281,7 +1281,7 @@ class PlasmaCompositionParam(NamedTuple):
 
     f_plasma_fuel_helium3: Any = None
 
-    nd_plasma_alphas_vol_avg: Any = None
+    nd_plasma_alphas_thermal_vol_avg: Any = None
 
     nd_plasma_electrons_vol_avg: Any = None
 
@@ -1386,7 +1386,7 @@ class PlasmaCompositionParam(NamedTuple):
             f_nd_plasma_carbon_electron=0,
             f_nd_plasma_oxygen_electron=0,
             f_alpha_ion=0,
-            f_nd_alpha_electron=0.10000000000000001,
+            f_nd_alpha_thermal_electron=0.10000000000000001,
             f_nd_beam_electron=0,
             n_charge_plasma_effective_vol_avg=0,
             nd_plasma_impurities_vol_avg=0,
@@ -1400,7 +1400,7 @@ class PlasmaCompositionParam(NamedTuple):
             alphan=1,
             nd_beam_ions=0,
             f_plasma_fuel_helium3=0,
-            nd_plasma_alphas_vol_avg=0,
+            nd_plasma_alphas_thermal_vol_avg=0,
             nd_plasma_electrons_vol_avg=7.5e19,
             nd_plasma_protons_vol_avg=0,
             first_call=1,
@@ -1492,7 +1492,7 @@ class PlasmaCompositionParam(NamedTuple):
             f_nd_plasma_carbon_electron=0,
             f_nd_plasma_oxygen_electron=0,
             f_alpha_ion=0.3154069116809366,
-            f_nd_alpha_electron=0.10000000000000001,
+            f_nd_alpha_thermal_electron=0.10000000000000001,
             f_nd_beam_electron=0,
             n_charge_plasma_effective_vol_avg=2.0909945616489103,
             nd_plasma_impurities_vol_avg=28875000000000004,
@@ -1506,7 +1506,7 @@ class PlasmaCompositionParam(NamedTuple):
             alphan=1,
             nd_beam_ions=0,
             f_plasma_fuel_helium3=0,
-            nd_plasma_alphas_vol_avg=7.5e18,
+            nd_plasma_alphas_thermal_vol_avg=7.5e18,
             nd_plasma_electrons_vol_avg=7.5e19,
             nd_plasma_protons_vol_avg=7500000000000000,
             first_call=0,
@@ -1588,7 +1588,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
         "f_nd_plasma_carbon_electron",
         "f_nd_plasma_oxygen_electron",
         "f_alpha_ion",
-        "f_nd_alpha_electron",
+        "f_nd_alpha_thermal_electron",
         "f_nd_beam_electron",
         "n_charge_plasma_effective_vol_avg",
         "nd_plasma_impurities_vol_avg",
@@ -1602,7 +1602,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
         "alphan",
         "nd_beam_ions",
         "f_plasma_fuel_helium3",
-        "nd_plasma_alphas_vol_avg",
+        "nd_plasma_alphas_thermal_vol_avg",
         "nd_plasma_electrons_vol_avg",
         "nd_plasma_protons_vol_avg",
         "first_call",
@@ -1658,7 +1658,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
         plasmacompositionparam.expected_m_beam_amu
     )
 
-    assert physics.data.physics.nd_plasma_alphas_vol_avg == pytest.approx(
+    assert physics.data.physics.nd_plasma_alphas_thermal_vol_avg == pytest.approx(
         plasmacompositionparam.expected_nd_alphas
     )
 
@@ -1837,7 +1837,7 @@ class PhyauxParam(NamedTuple):
 
     nd_plasma_fuel_ions_vol_avg: Any = None
 
-    nd_plasma_alphas_vol_avg: Any = None
+    nd_plasma_alphas_thermal_vol_avg: Any = None
 
     fusden_total: Any = None
 
@@ -1876,7 +1876,7 @@ class PhyauxParam(NamedTuple):
             burnup_in=0,
             aspect=3,
             nd_plasma_fuel_ions_vol_avg=5.8589175702454272e19,
-            nd_plasma_alphas_vol_avg=7.5e18,
+            nd_plasma_alphas_thermal_vol_avg=7.5e18,
             fusden_total=1.9852091609123786e17,
             fusden_alpha_total=1.973996644759543e17,
             plasma_current=18398455.678867526,
@@ -1895,7 +1895,7 @@ class PhyauxParam(NamedTuple):
             burnup_in=0,
             aspect=3,
             nd_plasma_fuel_ions_vol_avg=5.8576156204039725e19,
-            nd_plasma_alphas_vol_avg=7.5e18,
+            nd_plasma_alphas_thermal_vol_avg=7.5e18,
             fusden_total=1.9843269653375773e17,
             fusden_alpha_total=1.9731194318497056e17,
             plasma_current=18398455.678867526,
@@ -1940,7 +1940,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics):
     ) = physics.phyaux(
         aspect=phyauxparam.aspect,
         nd_plasma_fuel_ions_vol_avg=phyauxparam.nd_plasma_fuel_ions_vol_avg,
-        nd_plasma_alphas_vol_avg=phyauxparam.nd_plasma_alphas_vol_avg,
+        nd_plasma_alphas_thermal_vol_avg=phyauxparam.nd_plasma_alphas_thermal_vol_avg,
         fusden_total=phyauxparam.fusden_total,
         fusden_alpha_total=phyauxparam.fusden_alpha_total,
         plasma_current=phyauxparam.plasma_current,
@@ -3183,7 +3183,7 @@ def test_calculate_plasma_masses(physics):
     m_ions_total_amu = 3.0
     nd_plasma_ions_total_vol_avg = 1.0e20
     nd_plasma_fuel_ions_vol_avg = 0.8e20
-    nd_plasma_alphas_vol_avg = 0.1e20
+    nd_plasma_alphas_thermal_vol_avg = 0.1e20
     vol_plasma = 100.0
     nd_plasma_electrons_vol_avg = 1.0e20
 
@@ -3198,7 +3198,7 @@ def test_calculate_plasma_masses(physics):
         m_ions_total_amu=m_ions_total_amu,
         nd_plasma_ions_total_vol_avg=nd_plasma_ions_total_vol_avg,
         nd_plasma_fuel_ions_vol_avg=nd_plasma_fuel_ions_vol_avg,
-        nd_plasma_alphas_vol_avg=nd_plasma_alphas_vol_avg,
+        nd_plasma_alphas_thermal_vol_avg=nd_plasma_alphas_thermal_vol_avg,
         vol_plasma=vol_plasma,
         nd_plasma_electrons_vol_avg=nd_plasma_electrons_vol_avg,
     )
@@ -3489,7 +3489,7 @@ def test_detailed_physics_run_computes_profiles(monkeypatch, physics, process_mo
     )
     monkeypatch.setattr(
         physics.data.physics,
-        "nd_plasma_alphas_vol_avg",
+        "nd_plasma_alphas_thermal_vol_avg",
         np.mean(plasma.neprofile.profile_y) * 0.1,
     )
 
diff --git a/tests/unit/models/stellarator/test_neoclassics.py b/tests/unit/models/stellarator/test_neoclassics.py
index e789d0dc39..bdc5cbaa1a 100644
--- a/tests/unit/models/stellarator/test_neoclassics.py
+++ b/tests/unit/models/stellarator/test_neoclassics.py
@@ -22,7 +22,7 @@ class InitNeoclassicsParam(NamedTuple):
     temp_plasma_ion_on_axis_kev: Any = None
     nd_plasma_ions_on_axis: Any = None
     f_plasma_fuel_deuterium: Any = None
-    nd_plasma_alphas_vol_avg: Any = None
+    nd_plasma_alphas_thermal_vol_avg: Any = None
     rminor: Any = None
     rmajor: Any = None
     b_plasma_toroidal_on_axis: Any = None
@@ -71,7 +71,7 @@ class InitNeoclassicsParam(NamedTuple):
             temp_plasma_ion_on_axis_kev=12.579710000000002,
             nd_plasma_ions_on_axis=2.3930858160000005e20,
             f_plasma_fuel_deuterium=0.5,
-            nd_plasma_alphas_vol_avg=2.9820384000000004e19,
+            nd_plasma_alphas_thermal_vol_avg=2.9820384000000004e19,
             rminor=1.7993820274145451,
             rmajor=22.16,
             b_plasma_toroidal_on_axis=5.2400000000000002,
@@ -997,8 +997,8 @@ def test_init_neoclassics(initneoclassicsparam, monkeypatch, neoclassics):
     )
     monkeypatch.setattr(
         neoclassics.data.physics,
-        "nd_plasma_alphas_vol_avg",
-        initneoclassicsparam.nd_plasma_alphas_vol_avg,
+        "nd_plasma_alphas_thermal_vol_avg",
+        initneoclassicsparam.nd_plasma_alphas_thermal_vol_avg,
     )
     monkeypatch.setattr(neoclassics.data.physics, "rminor", initneoclassicsparam.rminor)
     monkeypatch.setattr(neoclassics.data.physics, "rmajor", initneoclassicsparam.rmajor)
diff --git a/tracking/tracking_data.py b/tracking/tracking_data.py
index ba6d05c0bb..68608fbd5c 100644
--- a/tracking/tracking_data.py
+++ b/tracking/tracking_data.py
@@ -90,7 +90,7 @@
     "Physics.triang95",
     "Physics.p_plasma_inner_rad_mw",
     "Physics.temp_plasma_separatrix_kev",
-    "Physics.f_nd_alpha_electron",
+    "Physics.f_nd_alpha_thermal_electron",
     "Physics.pflux_fw_neutron_mw",
     "Physics.aspect",
     "Physics.rminor",