diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md
index 6032e13737..d5a772bb23 100644
--- a/documentation/source/physics-models/profiles/plasma_profiles.md
+++ b/documentation/source/physics-models/profiles/plasma_profiles.md
@@ -636,6 +636,47 @@ The same function is run from the `i_plasma_pedestal == 0 ` profile case, found
 
 -----
 
+
+
+## Calculate profile volume average | `calculate_vol_avg_of_profile()`
+
+General method to calculate the volume averaged value of any normalised radial profile.
+
+### Derivation 
+
+The integral $I$ of any quantity $G(\rho)$ over the volume of the plasma where $\rho$ is the normalised minor radius is:
+
+$$
+I=\int_{}^{} G(\rho) \ dV
+$$
+
+The volume element is the surface area of the toroidal surface times the thickness $dV=2\pi R_0 2\pi\kappa r  dr  $
+
+$$
+I=\int_{0}^{1} G(\rho) 2\pi R_0 2\pi\kappa r  \ dr 
+$$
+
+Substituting the normalised minor radius coordinate $\rho=r/a$, and bringing the constants out in front,
+
+$$
+I=4 \pi^2 R_0 \kappa a^2 \int_{0}^{1} G(\rho)  \rho  \ d\rho 
+$$
+
+However, if we calculate the volume $V$, the exact same integral is used where $G=1$, 
+
+$$
+V=4 \pi^2 R_0 \kappa a^2 \int_{0}^{1}   \rho  \ d\rho = \frac{1}{2}  4 \pi^2 R_0 \kappa a^2
+$$
+
+Therefore the general volume integral of  $G(\rho)$, where $V$ is the volume, is
+
+$$
+I=2 V \int_{0}^{1} G(\rho)  \rho  \ d\rho 
+$$
+
+If $G$ is a density then $2  \int_{0}^{1} G(\rho)  \rho  \ d\rho$ is the volume averaged density.
+
+
 ## Key Constraints
 
 --------
diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py
index 53342d9abc..d782e19b68 100644
--- a/process/models/physics/profiles.py
+++ b/process/models/physics/profiles.py
@@ -558,3 +558,52 @@ def set_physics_variables(self):
             / self.data.physics.temp_plasma_electron_vol_avg_kev
             * self.data.physics.temp_plasma_electron_on_axis_kev
         )
+
+
+def calculate_vol_avg_of_profile(profile_x: np.ndarray, profile_y: np.ndarray) -> float:
+    """Calculate the volume averaged value (⟨profile_y⟩) of a radially normalised
+    profile.
+
+    Parameters
+    ----------
+    profile_x :
+        The x-values of the profile.
+    profile_y :
+        The y-values of the profile.
+
+    Returns
+    -------
+    float
+        The volume-averaged value (⟨profile_y⟩) of the profile.
+
+    Raises
+    ------
+    ValueError
+        If profile_x is not a 1D array, contains fewer than 2 points,
+        does not span from 0 to 1,or is not strictly increasing.
+
+
+    Notes
+    -----
+    - The 2 factor in the calculation arises using both sides of the profile and the
+    radial normalisation of the profile.
+
+
+    """
+    if profile_x.ndim != 1:
+        raise ValueError("profile_x must be a 1D array.")
+
+    if profile_x.size < 2:
+        raise ValueError("profile_x must contain at least 2 points.")
+
+    if not np.isclose(profile_x[0], 0.0) or not np.isclose(profile_x[-1], 1.0):
+        raise ValueError("profile_x must span from 0 to 1.")
+
+    if np.any(np.diff(profile_x) <= 0):
+        raise ValueError("profile_x must be strictly increasing.")
+
+    return 2.0 * sp.integrate.simpson(
+        profile_y * profile_x,
+        x=profile_x,
+        dx=profile_x[1] - profile_x[0],
+    )