From 42e8ab93a6a52e4d9f1600279ea918a215eeb5d8 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Fri, 7 Aug 2026 11:58:02 +0100
Subject: [PATCH 01/18] Add calculation for outboard midplane near SOL radial
 profile and update output message

---
 process/models/physics/scrape_off_layer.py | 58 ++++++++++++++++++++++
 1 file changed, 58 insertions(+)

diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 4388d79914..45b8ec82e9 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -344,3 +344,61 @@ def calculate_upstream_sol_outboard_parallel_area(
             * len_plasma_sol_power_decay
             * (b_plasma_surface_poloidal_average / b_plasma_outboard_total)
         )
+
+    @staticmethod
+    def calculate_outboard_midplane_near_sol_radial_profile(
+        rmajor: float,
+        rminor: float,
+        len_plasma_sol_power_decay: float,
+        pflux_plasma_outboard_sol_parallel_mw: float,
+        r: float | np.ndarray,
+    ) -> float | np.ndarray:
+        """Calculate the outboard midplane near SOL radial profile (qₗₗ(r)) [MW/m²].
+
+        Parameters
+        ----------
+        rmajor : float
+            Major radius of the plasma (R₀) [m]
+        rminor : float
+            Minor radius of the plasma (a) [m]
+        len_plasma_sol_power_decay : float
+            Power decay length (λ_q) [m]
+        pflux_plasma_outboard_sol_parallel_mw : float
+            Parallel power flux at the outboard midplane (qₗₗ,ᵤ) [MW/m²]
+        r : float|np.ndarray
+            Radial position(s) at which to calculate the SOL profile [m]
+
+        Returns
+        -------
+        float|np.ndarray
+            Outboard midplane SOL radial profile (qₗₗ(r)) [MW/m²]
+
+        Notes
+        -----
+        - The exponential model is highly valid in the "near-SOL" (typically the first
+        few millimeters to a centimeter outside the separatrix). In this region, parallel
+        heat transport is dominated by classical electron heat conduction
+        (Spitzer-Härm conductivity), which is vastly faster than perpendicular diffusion.
+        This competition between fast parallel conduction and slow perpendicular
+        diffusion naturally produces an exponential radial profile.
+
+        - The midplane exponential assumes steady-state H-mode conditions without the
+        massive, transient convective bursts caused by ELMs, which momentarily
+        flatten the entire midplane profile.
+
+        References
+        ----------
+        [1] T. Eich et al., “Scaling of the tokamak near the scrape-off layer H-mode
+        power width and implications for ITER,” Nuclear Fusion, vol. 53, no. 9,
+        p. 093031, Aug. 2013, doi: 10.1088/0029-5515/53/9/093031.
+
+        """
+        if r < (rmajor + rminor):
+            raise ValueError(
+                f"Radial position r={r} must be greater than or equal to the plasma "
+                f"edge (rmajor + rminor)={rmajor + rminor}."
+            )
+
+        return pflux_plasma_outboard_sol_parallel_mw * np.exp(
+            -(r - (rmajor + rminor)) / len_plasma_sol_power_decay
+        )

From 5ddbedc39047825f9082131c1de772c0b188b504 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Fri, 7 Aug 2026 13:43:47 +0100
Subject: [PATCH 02/18] Add function to plot midplane near SOL radial profile
 and update scrape off layer validation

---
 process/core/io/plot/summary.py            | 49 ++++++++++++++++++++++
 process/models/physics/scrape_off_layer.py |  2 +-
 2 files changed, 50 insertions(+), 1 deletion(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index bd34580d43..e30d44f74b 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -85,6 +85,7 @@
 )
 from process.models.physics.profiles import PlasmaProfileShapeType
 from process.models.pulse import PulseTimings
+from process.models.physics.scrape_off_layer import ScrapeOffLayer
 from process.models.superconductors import SuperconductorModel
 from process.models.tfcoil.base import (
     TFCoilShapeModel,
@@ -9349,6 +9350,50 @@ def make_bbox_props(power: float) -> dict[str, Any]:
     axis.get_yaxis().set_ticks([])
 
 
+def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: int):
+    """Function to plot the radial profile of the near SOL at the midplane."""
+    rmajor = mfile.get("rmajor", scan=scan)
+    rminor = mfile.get("rminor", scan=scan)
+    len_plasma_sol_power_decay = mfile.get(
+        "len_plasma_sol_eich13_power_decay", scan=scan
+    )
+    r = np.linspace(
+        (rmajor + rminor), (rmajor + rminor) + (3 * len_plasma_sol_power_decay), 100
+    )
+
+    radial_profile = (
+        ScrapeOffLayer().calculate_outboard_midplane_near_sol_radial_profile(
+            rmajor=rmajor,
+            rminor=rminor,
+            len_plasma_sol_power_decay=mfile.get(
+                "len_plasma_sol_eich13_power_decay", scan=scan
+            ),
+            pflux_plasma_outboard_sol_parallel_mw=mfile.get(
+                "pflux_plasma_outboard_sol_eich13_parallel_mw", scan=scan
+            ),
+            r=r,
+        )
+    )
+
+    axis.axvline(
+        x=rmajor + rminor + len_plasma_sol_power_decay,
+        color="k",
+        linestyle="--",
+        label=r"$\lambda_q$",
+    )
+
+    axis.set_xlim([
+        rmajor + rminor,
+        (rmajor + rminor) + (3 * len_plasma_sol_power_decay),
+    ])
+    axis.plot(r, radial_profile)
+    axis.grid()
+    axis.legend()
+    axis.set_title(r"Upstream Near SOL $q_{\parallel}$ Radial Profile")
+    axis.set_xlabel("Radial Position [m]")
+    axis.set_ylabel(r"$q_{\parallel}$ [MW/m$^2$]")
+
+
 def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed=None):
     """Function to plot a scatter box plot of L-H threshold power comparisons.
 
@@ -16731,6 +16776,10 @@ def _add_page(name: str | None = None):
         pages["plasma_exhaust"].add_subplot(122), m_file, scan, colour_scheme
     )
 
+    plot_midplane_near_sol_radial_profile(
+        _add_page("midplane_near_sol_radial_profile").add_subplot(111), m_file, scan
+    )
+
     plot_debye_length_profile(
         _add_page("microscopic_quantities").add_subplot(232), m_file, scan
     )
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 45b8ec82e9..99aee0bed5 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -393,7 +393,7 @@ def calculate_outboard_midplane_near_sol_radial_profile(
         p. 093031, Aug. 2013, doi: 10.1088/0029-5515/53/9/093031.
 
         """
-        if r < (rmajor + rminor):
+        if np.any(r < (rmajor + rminor)):
             raise ValueError(
                 f"Radial position r={r} must be greater than or equal to the plasma "
                 f"edge (rmajor + rminor)={rmajor + rminor}."

From d844c3a982091ca08c92f29f43f6ab73f6630c95 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 17 Aug 2026 13:30:54 +0100
Subject: [PATCH 03/18] Add Eich target heat flux profile calculation to
 ScrapeOffLayer model

---
 process/core/io/plot/summary.py            |  1 +
 process/models/physics/scrape_off_layer.py | 67 ++++++++++++++++++++++
 2 files changed, 68 insertions(+)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index e30d44f74b..f5f34f17cb 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -86,6 +86,7 @@
 from process.models.physics.profiles import PlasmaProfileShapeType
 from process.models.pulse import PulseTimings
 from process.models.physics.scrape_off_layer import ScrapeOffLayer
+from process.models.pulse import PulseTimings
 from process.models.superconductors import SuperconductorModel
 from process.models.tfcoil.base import (
     TFCoilShapeModel,
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 99aee0bed5..408715d737 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -3,6 +3,7 @@
 import logging
 
 import numpy as np
+import scipy
 
 from process.core import constants
 from process.core import process_output as po
@@ -402,3 +403,69 @@ def calculate_outboard_midplane_near_sol_radial_profile(
         return pflux_plasma_outboard_sol_parallel_mw * np.exp(
             -(r - (rmajor + rminor)) / len_plasma_sol_power_decay
         )
+
+    @staticmethod
+    def calculate_eich_target_heat_flux_profile(
+        pflux_plasma_sol_parallel_mw: float,
+        len_plasma_sol_power_decay: float,
+        f_b_div_flux_expansion: float,
+        len_plasma_sol_power_spreading: float,
+        plux_target_background_heat_flux_mw: float,
+        r: float | np.ndarray,
+    ) -> float | np.ndarray:
+        """Calculate the Eich target heat flux profile (qₜ(r)) [MW/m²].
+
+        Parameters
+        ----------
+        pflux_plasma_sol_parallel_mw : float
+            Parallel power flux at the outboard midplane (qₗₗ,ᵤ) [MW/m²]
+        len_plasma_sol_power_decay : float
+            Power decay length (λ_q) [m]
+        f_b_div_flux_expansion : float
+            Divertor flux expansion factor (fₓ) [-]
+        len_plasma_sol_power_spreading : float
+            Power spreading length in the divertor (S) [m]
+        plux_target_background_heat_flux_mw : float
+            Background heat flux at the divertor target [MW/m²]
+        r : float|np.ndarray
+            Radial position(s) at which to calculate the target heat flux profile [m]
+
+        Returns
+        -------
+        float|np.ndarray
+            Eich target heat flux profile (qₜ(r)) [MW/m²]
+
+        Notes
+        -----
+        - The Eich target heat flux profile is derived from the midplane exponential
+        profile, taking into account the magnetic geometry and flux expansion between
+        the midplane and the divertor target. The profile is typically characterized by
+        a combination of an exponential decay and a Gaussian spreading due to cross-field
+        transport in the divertor leg.
+
+        References
+        ----------
+        [1] T. Eich, B. Sieglin, A. Scarabosio, W. Fundamenski, R. J. Goldston, and
+        A. Herrmann, “Inter-ELM Power Decay Length for JET and ASDEX Upgrade: Measurement
+        and Comparison with Heuristic Drift-Based Model,” Physical Review Letters,
+        vol. 107, no. 21, Nov. 2011, doi: https://doi.org/10.1103/PhysRevLett.107.215001
+
+        [2] T. Eich et al., “Scaling of the tokamak near the scrape-off layer H-mode
+        power width and implications for ITER,” Nuclear Fusion, vol. 53, no. 9,
+        p. 093031, Aug. 2013, doi: 10.1088/0029-5515/53/9/093031.
+
+        """
+        return (pflux_plasma_sol_parallel_mw / 2) * np.exp(
+            (
+                (len_plasma_sol_power_spreading)
+                / (2 * len_plasma_sol_power_decay * f_b_div_flux_expansion)
+            )
+            ** 2
+            - (r / (len_plasma_sol_power_spreading * f_b_div_flux_expansion))
+        ) * scipy.special.erfc(
+            (
+                len_plasma_sol_power_spreading
+                / (2 * len_plasma_sol_power_decay * f_b_div_flux_expansion)
+            )
+            - (r / (len_plasma_sol_power_spreading))
+        ) + plux_target_background_heat_flux_mw

From 29b55c68e8ce220ebdafa129cd443ff1de0dbe6e Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 17 Aug 2026 14:31:07 +0100
Subject: [PATCH 04/18] Add function to plot lower outboard Eich target heat
 flux profile and update calculation method

---
 process/core/io/plot/summary.py            | 59 +++++++++++++++++-----
 process/models/physics/scrape_off_layer.py | 18 ++++++-
 2 files changed, 62 insertions(+), 15 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index f5f34f17cb..db5b6c41b3 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9375,24 +9375,53 @@ def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: in
             r=r,
         )
     )
+    
+    axis.plot(r, radial_profile)
+    axis.grid()
+    axis.set_title(r"Midplane Near SOL Radial Profile")
+    axis.set_xlabel("Radial Position [m]")
+    axis.set_ylabel(r"$q_{||}$ [MW/m$^2$]")
 
-    axis.axvline(
-        x=rmajor + rminor + len_plasma_sol_power_decay,
-        color="k",
-        linestyle="--",
-        label=r"$\lambda_q$",
+
+def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, scan: int):
+    """Function to plot the Eich target profile at the lower outboard divertor."""
+    rmajor = mfile.get("rmajor", scan=scan)
+    rminor = mfile.get("rminor", scan=scan)
+    len_plasma_sol_power_decay = mfile.get(
+        "len_plasma_sol_eich13_power_decay", scan=scan
+    )
+    f_b_flux_expansion = 5.0
+    r = np.linspace(
+        (rmajor + rminor)- (1.5 * len_plasma_sol_power_decay) * f_b_flux_expansion,
+        (rmajor + rminor) + (3 * len_plasma_sol_power_decay) * f_b_flux_expansion,
+        200,
     )
 
-    axis.set_xlim([
-        rmajor + rminor,
-        (rmajor + rminor) + (3 * len_plasma_sol_power_decay),
-    ])
-    axis.plot(r, radial_profile)
+    pflux_target_profile = ScrapeOffLayer().calculate_eich_target_heat_flux_profile(
+        rmajor=rmajor,
+        rminor=rminor,
+        pflux_plasma_sol_parallel_mw=mfile.get(
+            "pflux_plasma_outboard_sol_parallel_mw", scan=scan
+        ),
+        len_plasma_sol_power_decay=mfile.get("len_sol_outboard_power_decay", scan=scan),
+        f_b_div_flux_expansion=f_b_flux_expansion,
+        len_plasma_sol_power_spreading=1.5e-3,
+        plux_target_background_heat_flux_mw=0.0,
+        r=r,
+    )
+    peak_idx = np.argmax(pflux_target_profile)
+    peak_r = r[peak_idx]
+    peak_q = pflux_target_profile[peak_idx]
+
+    axis.plot(r, pflux_target_profile)
+    axis.axvline(peak_r, color="black", linestyle="--", linewidth=1)
+    axis.axhline(peak_q, color="black", linestyle="--", linewidth=1)
     axis.grid()
     axis.legend()
-    axis.set_title(r"Upstream Near SOL $q_{\parallel}$ Radial Profile")
+    axis.minorticks_on()
+    axis.set_title(r"Lower Outboard Eich Target Heat Flux Profile")
     axis.set_xlabel("Radial Position [m]")
-    axis.set_ylabel(r"$q_{\parallel}$ [MW/m$^2$]")
+    axis.set_ylabel(r"$q_{||,t}$ [MW/m$^2$]")
 
 
 def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed=None):
@@ -16778,7 +16807,11 @@ def _add_page(name: str | None = None):
     )
 
     plot_midplane_near_sol_radial_profile(
-        _add_page("midplane_near_sol_radial_profile").add_subplot(111), m_file, scan
+        _add_page("midplane_near_sol_radial_profile").add_subplot(121), m_file, scan
+    )
+
+    plot_div_lower_outboard_eich_target_profile(
+        pages["midplane_near_sol_radial_profile"].add_subplot(122), m_file, scan
     )
 
     plot_debye_length_profile(
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 408715d737..ed06056feb 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -406,6 +406,8 @@ def calculate_outboard_midplane_near_sol_radial_profile(
 
     @staticmethod
     def calculate_eich_target_heat_flux_profile(
+        rmajor: float,
+        rminor: float,
         pflux_plasma_sol_parallel_mw: float,
         len_plasma_sol_power_decay: float,
         f_b_div_flux_expansion: float,
@@ -417,6 +419,10 @@ def calculate_eich_target_heat_flux_profile(
 
         Parameters
         ----------
+        rmajor : float
+            Major radius of the plasma (R₀) [m]
+        rminor : float
+            Minor radius of the plasma (a) [m]
         pflux_plasma_sol_parallel_mw : float
             Parallel power flux at the outboard midplane (qₗₗ,ᵤ) [MW/m²]
         len_plasma_sol_power_decay : float
@@ -461,11 +467,19 @@ def calculate_eich_target_heat_flux_profile(
                 / (2 * len_plasma_sol_power_decay * f_b_div_flux_expansion)
             )
             ** 2
-            - (r / (len_plasma_sol_power_spreading * f_b_div_flux_expansion))
+            - (
+                (r - (rmajor + rminor))
+                * f_b_div_flux_expansion
+                / (len_plasma_sol_power_spreading * f_b_div_flux_expansion)
+            )
         ) * scipy.special.erfc(
             (
                 len_plasma_sol_power_spreading
                 / (2 * len_plasma_sol_power_decay * f_b_div_flux_expansion)
             )
-            - (r / (len_plasma_sol_power_spreading))
+            - (
+                (r - (rmajor + rminor))
+                * f_b_div_flux_expansion
+                / (len_plasma_sol_power_spreading)
+            )
         ) + plux_target_background_heat_flux_mw

From 3075bd534885beadb43805bbf80d0d8a67cabf1e Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 18 Aug 2026 08:53:28 +0100
Subject: [PATCH 05/18] Add function to plot separatrix power flux profiles and
 update main plot structure

---
 process/core/io/plot/summary.py | 54 +++++++++++++++++++++++++++++++--
 1 file changed, 51 insertions(+), 3 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index db5b6c41b3..03c5bce635 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9423,6 +9423,47 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     axis.set_xlabel("Radial Position [m]")
     axis.set_ylabel(r"$q_{||,t}$ [MW/m$^2$]")
 
+def plot_sol_power_flux_profiles(axis: plt.Axes, mfile: MFile, scan: int, colour_scheme):
+    """Plot separatrix power split fractions as a bar chart."""
+    plot_plasma(axis=axis, mfile=mfile, scan=scan, colour_scheme=colour_scheme)
+    rmajor, rminor,kappa= mfile.get_variables(
+        "rmajor",
+        "rminor",
+        "kappa",
+        scan=scan,
+    )
+
+    plasma_scale = max(rminor, abs(kappa * rminor), 1e-6)
+    scale_factor = min(max(plasma_scale / 2.0, 0.7), 1.0)
+    text_fontsize = 9 * scale_factor
+
+
+    
+    outboard_pos = (rmajor + rminor, 0.0)
+
+    axis.text(
+        *outboard_pos,
+        f"$f_{{\\mathrm{{outboard}}}} = {5:.3f}$\n"
+        f"$\\Delta r_{{\\mathrm{{sep}}}} = {6:.3f}$ m",
+        fontsize=text_fontsize,
+        verticalalignment="center",
+        horizontalalignment="center",
+        bbox={
+            "boxstyle": f"round,pad={0.3 * scale_factor:.3f}",
+            "alpha": 1.0,
+            "linewidth": 2 * scale_factor,
+            "edgecolor": "black",
+        },
+        zorder=101,
+    )
+    
+
+    axis.spines["top"].set_visible(False)
+    axis.spines["right"].set_visible(False)
+    axis.spines["bottom"].set_visible(False)
+    axis.spines["left"].set_visible(False)
+    axis.get_xaxis().set_ticks([])
+    axis.get_yaxis().set_ticks([])
 
 def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed=None):
     """Function to plot a scatter box plot of L-H threshold power comparisons.
@@ -16806,14 +16847,21 @@ def _add_page(name: str | None = None):
         pages["plasma_exhaust"].add_subplot(122), m_file, scan, colour_scheme
     )
 
+    
+    
+    plot_sol_power_flux_profiles(_add_page("sol_powerfluxes").add_subplot(121), m_file, scan, colour_scheme)
+
+    
+    fig, (ax1, ax2) = plt.subplots(2, sharex=True)
+    
     plot_midplane_near_sol_radial_profile(
-        _add_page("midplane_near_sol_radial_profile").add_subplot(121), m_file, scan
+        pages["sol_powerfluxes"].add_subplot(324), m_file, scan
     )
 
     plot_div_lower_outboard_eich_target_profile(
-        pages["midplane_near_sol_radial_profile"].add_subplot(122), m_file, scan
+        pages["sol_powerfluxes"].add_subplot(326), m_file, scan
     )
-
+    
     plot_debye_length_profile(
         _add_page("microscopic_quantities").add_subplot(232), m_file, scan
     )

From 8ffc70db7853b02f511fa23ba1d61e5ac764af7c Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 25 Aug 2026 14:43:45 +0100
Subject: [PATCH 06/18] Update midplane near SOL radial profile and lower
 outboard Eich target profile plots with additional data and improved labeling

---
 process/core/io/plot/summary.py | 63 ++++++++++++++++++++++-----------
 1 file changed, 43 insertions(+), 20 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 03c5bce635..3c972ee697 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -84,7 +84,6 @@
     PlasmaShapeModelType,
 )
 from process.models.physics.profiles import PlasmaProfileShapeType
-from process.models.pulse import PulseTimings
 from process.models.physics.scrape_off_layer import ScrapeOffLayer
 from process.models.pulse import PulseTimings
 from process.models.superconductors import SuperconductorModel
@@ -9361,6 +9360,7 @@ def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: in
     r = np.linspace(
         (rmajor + rminor), (rmajor + rminor) + (3 * len_plasma_sol_power_decay), 100
     )
+    len_sol_outboard_power_decay = mfile.get("len_sol_outboard_power_decay", scan=scan)
 
     radial_profile = (
         ScrapeOffLayer().calculate_outboard_midplane_near_sol_radial_profile(
@@ -9375,11 +9375,22 @@ def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: in
             r=r,
         )
     )
-    
-    axis.plot(r, radial_profile)
+
+    x_ref = rmajor + rminor + len_sol_outboard_power_decay
+
+    axis.plot(r, radial_profile, label=r"$q_{||}$ profile")
+    axis.axvline(
+        x_ref,
+        color="black",
+        linestyle="--",
+        linewidth=1,
+        label=r"$\lambda_{q,\mathrm{out}}$",
+    )
     axis.grid()
+    axis.legend()
     axis.set_title(r"Midplane Near SOL Radial Profile")
     axis.set_xlabel("Radial Position [m]")
+    axis.minorticks_on()
     axis.set_ylabel(r"$q_{||}$ [MW/m$^2$]")
 
 
@@ -9392,7 +9403,7 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     )
     f_b_flux_expansion = 5.0
     r = np.linspace(
-        (rmajor + rminor)- (1.5 * len_plasma_sol_power_decay) * f_b_flux_expansion,
+        (rmajor + rminor) - (1.5 * len_plasma_sol_power_decay) * f_b_flux_expansion,
         (rmajor + rminor) + (3 * len_plasma_sol_power_decay) * f_b_flux_expansion,
         200,
     )
@@ -9416,6 +9427,15 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     axis.plot(r, pflux_target_profile)
     axis.axvline(peak_r, color="black", linestyle="--", linewidth=1)
     axis.axhline(peak_q, color="black", linestyle="--", linewidth=1)
+    axis.text(
+        0.02,
+        0.98,
+        rf"$f_x$ = {f_b_flux_expansion:.2f}",
+        transform=axis.transAxes,
+        ha="left",
+        va="top",
+        bbox={"boxstyle": "round", "facecolor": "white", "alpha": 0.8},
+    )
     axis.grid()
     axis.legend()
     axis.minorticks_on()
@@ -9423,28 +9443,34 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     axis.set_xlabel("Radial Position [m]")
     axis.set_ylabel(r"$q_{||,t}$ [MW/m$^2$]")
 
+
 def plot_sol_power_flux_profiles(axis: plt.Axes, mfile: MFile, scan: int, colour_scheme):
     """Plot separatrix power split fractions as a bar chart."""
     plot_plasma(axis=axis, mfile=mfile, scan=scan, colour_scheme=colour_scheme)
-    rmajor, rminor,kappa= mfile.get_variables(
+    rmajor, rminor, kappa = mfile.get_variables(
         "rmajor",
         "rminor",
         "kappa",
         scan=scan,
     )
-
+    len_sol_outboard_power_decay = mfile.get("len_sol_outboard_power_decay", scan=scan)
+    a_plasma_outboard_sol_parallel = mfile.get(
+        "a_plasma_outboard_sol_parallel", scan=scan
+    )
+    pflux_plasma_outboard_sol_parallel_mw = mfile.get(
+        "pflux_plasma_outboard_sol_parallel_mw", scan=scan
+    )
     plasma_scale = max(rminor, abs(kappa * rminor), 1e-6)
     scale_factor = min(max(plasma_scale / 2.0, 0.7), 1.0)
     text_fontsize = 9 * scale_factor
 
-
-    
     outboard_pos = (rmajor + rminor, 0.0)
 
     axis.text(
         *outboard_pos,
-        f"$f_{{\\mathrm{{outboard}}}} = {5:.3f}$\n"
-        f"$\\Delta r_{{\\mathrm{{sep}}}} = {6:.3f}$ m",
+        f"$\\lambda_q = {len_sol_outboard_power_decay * 1e3:.3f}$ mm\n"
+        f"$A_{{||}} = {a_plasma_outboard_sol_parallel:.4f}$ m$^2$\n"
+        f"$q_{{||}} = {pflux_plasma_outboard_sol_parallel_mw:,.2f}$ MW/m$^2$",
         fontsize=text_fontsize,
         verticalalignment="center",
         horizontalalignment="center",
@@ -9456,7 +9482,6 @@ def plot_sol_power_flux_profiles(axis: plt.Axes, mfile: MFile, scan: int, colour
         },
         zorder=101,
     )
-    
 
     axis.spines["top"].set_visible(False)
     axis.spines["right"].set_visible(False)
@@ -9465,6 +9490,7 @@ def plot_sol_power_flux_profiles(axis: plt.Axes, mfile: MFile, scan: int, colour
     axis.get_xaxis().set_ticks([])
     axis.get_yaxis().set_ticks([])
 
+
 def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed=None):
     """Function to plot a scatter box plot of L-H threshold power comparisons.
 
@@ -16847,21 +16873,18 @@ def _add_page(name: str | None = None):
         pages["plasma_exhaust"].add_subplot(122), m_file, scan, colour_scheme
     )
 
-    
-    
-    plot_sol_power_flux_profiles(_add_page("sol_powerfluxes").add_subplot(121), m_file, scan, colour_scheme)
+    plot_sol_power_flux_profiles(
+        _add_page("sol_powerfluxes").add_subplot(121), m_file, scan, colour_scheme
+    )
 
-    
-    fig, (ax1, ax2) = plt.subplots(2, sharex=True)
-    
     plot_midplane_near_sol_radial_profile(
-        pages["sol_powerfluxes"].add_subplot(324), m_file, scan
+        pages["sol_powerfluxes"].add_subplot(336), m_file, scan
     )
 
     plot_div_lower_outboard_eich_target_profile(
-        pages["sol_powerfluxes"].add_subplot(326), m_file, scan
+        pages["sol_powerfluxes"].add_subplot(339), m_file, scan
     )
-    
+
     plot_debye_length_profile(
         _add_page("microscopic_quantities").add_subplot(232), m_file, scan
     )

From e226cd57076567c548a83fe55c92ecdb887f3f14 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 25 Aug 2026 14:49:04 +0100
Subject: [PATCH 07/18] Add outboard lower divertor flux expansion factor to
 PhysicsData class

---
 process/data_structure/physics_variables.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 53d7bf343a..6e4e0c3f91 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1760,6 +1760,9 @@ class PhysicsData:
     - =2 MAST 2014 scaling 1
     - =3 MAST 2014 scaling 2
     """
+    
+    f_b_div_outboard_lower_flux_expansion: float = 5.0
+    """Outboard lower divertor flux expansion factor for the divertor targets (fₓ)"""
 
     dt_power_density_plasma: float = 0.0
     sigmav_dt_average: float = 0.0

From 0fb881e2c2f906a54e7413d317d3b3a034a9b1b4 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 25 Aug 2026 14:56:03 +0100
Subject: [PATCH 08/18] Add outboard lower divertor flux expansion factor to
 PhysicsData and update scrape off layer calculations

---
 process/data_structure/physics_variables.py | 2 +-
 process/models/physics/scrape_off_layer.py  | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 6e4e0c3f91..5276fc669c 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1760,7 +1760,7 @@ class PhysicsData:
     - =2 MAST 2014 scaling 1
     - =3 MAST 2014 scaling 2
     """
-    
+
     f_b_div_outboard_lower_flux_expansion: float = 5.0
     """Outboard lower divertor flux expansion factor for the divertor targets (fₓ)"""
 
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index ed06056feb..51703005ff 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -95,6 +95,7 @@ def run(self):
 
         self.data.physics.pflux_plasma_outboard_sol_parallel_mw = (
             self.data.physics.p_plasma_separatrix_mw
+            * self.data.physics.f_p_div_outboard_separatrix
             / self.data.physics.a_plasma_outboard_sol_parallel
         )
 

From ff2c96256d3a734b8bd03249a99ac382bc465666 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 25 Aug 2026 15:02:36 +0100
Subject: [PATCH 09/18] Add outboard lower divertor power spreading length
 factors to PhysicsData

---
 process/data_structure/physics_variables.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 5276fc669c..10d3241fae 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1764,6 +1764,14 @@ class PhysicsData:
     f_b_div_outboard_lower_flux_expansion: float = 5.0
     """Outboard lower divertor flux expansion factor for the divertor targets (fₓ)"""
 
+    len_div_outboard_lower_scrabosio14_power_spreading: float = 0.0
+    """Scrabosio 2014 H-mode power spreading length/factor in the scrape-off layer scaling
+    (S) [m]"""
+
+    len_div_outboard_lower_power_spreading: float = 0.0
+    """Power spreading length/factor at the outboard lower divertor target
+    (S) [m]"""
+
     dt_power_density_plasma: float = 0.0
     sigmav_dt_average: float = 0.0
     dhe3_power_density: float = 0.0

From f2eaf28652ecb58bb301e68246d26eb42cc58ae1 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 25 Aug 2026 15:04:08 +0100
Subject: [PATCH 10/18] Refactor docstrings for outboard lower divertor power
 spreading length factors in PhysicsData

---
 process/data_structure/physics_variables.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 10d3241fae..1f0a5608b7 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1765,12 +1765,10 @@ class PhysicsData:
     """Outboard lower divertor flux expansion factor for the divertor targets (fₓ)"""
 
     len_div_outboard_lower_scrabosio14_power_spreading: float = 0.0
-    """Scrabosio 2014 H-mode power spreading length/factor in the scrape-off layer scaling
-    (S) [m]"""
+    """Scrabosio 2014 H-mode power spreading length/factor in the scrape-off layer scaling (S) [m]"""
 
     len_div_outboard_lower_power_spreading: float = 0.0
-    """Power spreading length/factor at the outboard lower divertor target
-    (S) [m]"""
+    """Power spreading length/factor at the outboard lower divertor target (S) [m]"""
 
     dt_power_density_plasma: float = 0.0
     sigmav_dt_average: float = 0.0

From 38e417b4a9c25d70cfa5e3404a6d1222b41f3c97 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 25 Aug 2026 15:16:38 +0100
Subject: [PATCH 11/18] Add Scrabosio 2014 power spreading factor calculation
 to ScrapeOffLayer model

---
 .../physics-models/plasma_scrape_off_layer.md | 31 +++++++-
 process/models/physics/scrape_off_layer.py    | 73 +++++++++++++++++++
 2 files changed, 103 insertions(+), 1 deletion(-)

diff --git a/documentation/source/physics-models/plasma_scrape_off_layer.md b/documentation/source/physics-models/plasma_scrape_off_layer.md
index a85fd981f6..5fc7318d8c 100644
--- a/documentation/source/physics-models/plasma_scrape_off_layer.md
+++ b/documentation/source/physics-models/plasma_scrape_off_layer.md
@@ -83,6 +83,33 @@ The $R^2$ value for this fit is 0.55
 
 --------
 
+## Spreading Parameter
+
+The scrape-off layer (SOL) spreading parameter $S$ represents a Gaussian width that quantifies additional perpendicular heat spreading in the divertor leg. It works alongside the upstream heat flux decay length $\lambda_{q}$ to determine total target heat loads on the divertor.
+
+Unlike $\lambda_{q}$, which is governed by robust upstream parallel and perpendicular transport physics at the plasma midplane, $S$ is inherently a "local" divertor parameter. Deriving a single, absolute multi-machine formula for $S$ is incredibly difficult due to several overlapping regional variables:
+
+- Divertor Geometry: The path length from the X-point to the target tile heavily impacts how much the heat spreads radially.
+
+- Plasma Recycling Regimes: Low-recycling, high-recycling, and detached plasma conditions completely alter the cross-field diffusion rates.
+
+- Localized Radiation: Impurity seeding and neutral gas interactions dissipate power unevenly along the divertor leg, altering the effective Gaussian profile width.
+
+-----------
+
+### Scarabosio 2015 | `calculate_scarabosio2015_power_spreading_factor()`
+
+The H-mode SOL spreading factor, $S$ is given in $\text{m}$ by[^scarabosio_2015]:
+
+$$
+S = (0.12(\pm0.07)\times 10^{-3}) P_{\text{sep}}^{0.21(\pm0.11)}R_0^{0.71(\pm0.5)}B_{\text{p}}(a)^{-0.82(\pm0.27)}n_{\text{sep}}^{0.71(\pm0.5)}
+$$
+
+- This was fitted from ASDEX Upgrade and JET outer target data
+- The $R^2$ value of the regression fit was 0.65
+
+------------
+
 [^eich_2013]: T. Eich et al., “Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER,” Nuclear Fusion, vol. 53, no. 9 p. 093031, Aug. 2013, doi: 10.1088/0029-5515/53/9/093031.
 
 [^mast_2014]: A. J. Thornton and A. Kirk, “Scaling of the scrape-off layer width during inter-ELM H modes on MAST as measured by infrared thermography,”
@@ -90,4 +117,6 @@ Plasma Physics and Controlled Fusion, vol. 56, no. 5, p. 055008, Apr. 2014, doi:
 
 [^stangeby_boundary]: P. C. Stangeby, “The Plasma Boundary of Magnetic Fusion Devices,” Jan. 2000, doi: 10.1201/9780367801489.
 
-[^henderson_step]: S. S. Henderson et al., “An overview of the STEP divertor design and the simple models driving the plasma exhaust scenario,” Nuclear Fusion, vol. 65, no. 1, pp. 016033–016033, Nov. 2024, doi: 10.1088/1741-4326/ad93e7.
\ No newline at end of file
+[^henderson_step]: S. S. Henderson et al., “An overview of the STEP divertor design and the simple models driving the plasma exhaust scenario,” Nuclear Fusion, vol. 65, no. 1, pp. 016033–016033, Nov. 2024, doi: 10.1088/1741-4326/ad93e7.
+
+[^scarabosio_2015]: A. Scarabosio et al., “Scaling of the divertor power spreading (S-factor) in open and closed divertor operation in JET and ASDEX Upgrade,” Journal of Nuclear Materials, vol. 463, pp. 49-54, Aug. 2015, doi: 10.1016/j.jnucmat.2014.11.076.
\ No newline at end of file
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 51703005ff..76dfcbe6ff 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -104,6 +104,18 @@ def run(self):
             / self.data.physics.a_plasma_outboard_sol_eich13_parallel
         )
 
+        self.data.physics.len_div_outboard_lower_scrabosio14_power_spreading = self.calculate_scarabosio2014_power_spreading_factor(  # noqa: E501
+            p_plasma_separatrix_mw=self.data.physics.p_plasma_separatrix_mw,
+            b_plasma_surface_poloidal_average=self.data.physics.b_plasma_surface_poloidal_average,
+            nd_plasma_separatrix_electron_19=self.data.physics.nd_plasma_separatrix_electron
+            / 1e19,
+            rmajor=self.data.physics.rmajor,
+        )
+
+        self.data.physics.len_div_outboard_lower_power_spreading = (
+            self.data.physics.len_div_outboard_lower_scrabosio14_power_spreading
+        )
+
     def output(self) -> None:
         """Output plasma scrape off layer physics information."""
         po.oheadr(self.outfile, "Plasma Scrape Off Layer")
@@ -189,6 +201,22 @@ def output(self) -> None:
             "(pflux_plasma_outboard_sol_eich13_parallel_mw)",
             self.data.physics.pflux_plasma_outboard_sol_eich13_parallel_mw,
         )
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
+        po.osubhd(self.outfile, "Power Spreading Factors (S):")
+
+        po.ovarre(
+            self.outfile,
+            "Outboard lower divertor power spreading factor (S) [m]",
+            "(len_div_outboard_lower_power_spreading)",
+            self.data.physics.len_div_outboard_lower_power_spreading,
+        )
+        po.ovarre(
+            self.outfile,
+            "Scrabosio 2014 H-mode power spreading factor (S) [m]",
+            "(len_div_outboard_lower_scrabosio14_power_spreading)",
+            self.data.physics.len_div_outboard_lower_scrabosio14_power_spreading,
+        )
 
     @staticmethod
     def calculate_eich2013_sol_power_decay_length(
@@ -484,3 +512,48 @@ def calculate_eich_target_heat_flux_profile(
                 / (len_plasma_sol_power_spreading)
             )
         ) + plux_target_background_heat_flux_mw
+
+    @staticmethod
+    def calculate_scarabosio2014_power_spreading_factor(
+        p_plasma_separatrix_mw: float,
+        b_plasma_surface_poloidal_average: float,
+        nd_plasma_separatrix_electron_19: float,
+        rmajor: float,
+    ) -> float:
+        """Calculate the Scrabosio 2014 H-mode power spreading factor (S).
+
+        Parameters
+        ----------
+        p_plasma_separatrix_mw : float
+            Power crossing the separatrix (Pₛₑₚ) [MW]
+        b_plasma_surface_poloidal_average : float
+            Poloidal magnetic field at the plasma surface (Bₚₒₗ(a))  [T]
+        nd_plasma_separatrix_electron_19 : float
+            Electron density at the separatrix (nₑ,ₛₑₚ) [10¹⁹ m⁻³]
+        rmajor : float
+            Major radius of the plasma (R₀) [m]
+
+        Returns
+        -------
+        float
+            Scrabosio 2014 H-mode power spreading factor (S) [m]
+
+        Notes
+        -----
+        - The R² for the fit is 0.65
+
+        References
+        ----------
+        [1] A. Scarabosio et al., “Scaling of the divertor power spreading (S-factor) in
+        open and closed divertor operation in JET and ASDEX Upgrade,”
+        Journal of Nuclear Materials, vol. 463, pp. 49-54, Aug. 2015,
+        doi: 10.1016/j.jnucmat.2014.11.076.
+
+        """
+        return (
+            0.12e-3
+            * p_plasma_separatrix_mw**0.21
+            * b_plasma_surface_poloidal_average**-0.82
+            * nd_plasma_separatrix_electron_19**-0.02
+            * rmajor**0.71
+        )

From 048eaeb64c66e3b81b3150631b643deb5241eb41 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Tue, 25 Aug 2026 15:53:00 +0100
Subject: [PATCH 12/18] Refactor power decay length variables in summary and
 scrape off layer models

---
 process/core/io/plot/summary.py            | 22 ++++++++++------------
 process/models/physics/scrape_off_layer.py | 12 ++++++------
 2 files changed, 16 insertions(+), 18 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 3c972ee697..d7b2c8284a 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9354,23 +9354,18 @@ def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: in
     """Function to plot the radial profile of the near SOL at the midplane."""
     rmajor = mfile.get("rmajor", scan=scan)
     rminor = mfile.get("rminor", scan=scan)
-    len_plasma_sol_power_decay = mfile.get(
-        "len_plasma_sol_eich13_power_decay", scan=scan
-    )
+    len_sol_outboard_power_decay = mfile.get("len_sol_outboard_power_decay", scan=scan)
     r = np.linspace(
-        (rmajor + rminor), (rmajor + rminor) + (3 * len_plasma_sol_power_decay), 100
+        (rmajor + rminor), (rmajor + rminor) + (3 * len_sol_outboard_power_decay), 100
     )
-    len_sol_outboard_power_decay = mfile.get("len_sol_outboard_power_decay", scan=scan)
 
     radial_profile = (
         ScrapeOffLayer().calculate_outboard_midplane_near_sol_radial_profile(
             rmajor=rmajor,
             rminor=rminor,
-            len_plasma_sol_power_decay=mfile.get(
-                "len_plasma_sol_eich13_power_decay", scan=scan
-            ),
+            len_plasma_sol_power_decay=len_sol_outboard_power_decay,
             pflux_plasma_outboard_sol_parallel_mw=mfile.get(
-                "pflux_plasma_outboard_sol_eich13_parallel_mw", scan=scan
+                "pflux_plasma_outboard_sol_parallel_mw", scan=scan
             ),
             r=r,
         )
@@ -9401,6 +9396,9 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     len_plasma_sol_power_decay = mfile.get(
         "len_plasma_sol_eich13_power_decay", scan=scan
     )
+    len_div_outboard_lower_power_spreading = mfile.get(
+        "len_div_outboard_lower_power_spreading", scan=scan
+    )
     f_b_flux_expansion = 5.0
     r = np.linspace(
         (rmajor + rminor) - (1.5 * len_plasma_sol_power_decay) * f_b_flux_expansion,
@@ -9416,8 +9414,8 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
         ),
         len_plasma_sol_power_decay=mfile.get("len_sol_outboard_power_decay", scan=scan),
         f_b_div_flux_expansion=f_b_flux_expansion,
-        len_plasma_sol_power_spreading=1.5e-3,
-        plux_target_background_heat_flux_mw=0.0,
+        len_plasma_sol_power_spreading=len_div_outboard_lower_power_spreading,
+        pflux_target_background_heat_flux_mw=0.0,
         r=r,
     )
     peak_idx = np.argmax(pflux_target_profile)
@@ -9430,7 +9428,7 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     axis.text(
         0.02,
         0.98,
-        rf"$f_x$ = {f_b_flux_expansion:.2f}",
+        f"$f_x$ = {f_b_flux_expansion:.2f}\n$S$ = {len_div_outboard_lower_power_spreading * 1e3:.3f} mm",
         transform=axis.transAxes,
         ha="left",
         va="top",
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 76dfcbe6ff..cecd77244f 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -441,10 +441,10 @@ def calculate_eich_target_heat_flux_profile(
         len_plasma_sol_power_decay: float,
         f_b_div_flux_expansion: float,
         len_plasma_sol_power_spreading: float,
-        plux_target_background_heat_flux_mw: float,
+        pflux_target_background_heat_flux_mw: float,
         r: float | np.ndarray,
     ) -> float | np.ndarray:
-        """Calculate the Eich target heat flux profile (qₜ(r)) [MW/m²].
+        """Calculate the Eich parallel target heat flux profile (qₗₗ,ₜ(r)) [MW/m²].
 
         Parameters
         ----------
@@ -460,7 +460,7 @@ def calculate_eich_target_heat_flux_profile(
             Divertor flux expansion factor (fₓ) [-]
         len_plasma_sol_power_spreading : float
             Power spreading length in the divertor (S) [m]
-        plux_target_background_heat_flux_mw : float
+        pflux_target_background_heat_flux_mw : float
             Background heat flux at the divertor target [MW/m²]
         r : float|np.ndarray
             Radial position(s) at which to calculate the target heat flux profile [m]
@@ -468,11 +468,11 @@ def calculate_eich_target_heat_flux_profile(
         Returns
         -------
         float|np.ndarray
-            Eich target heat flux profile (qₜ(r)) [MW/m²]
+            Eich parallel target heat flux profile (qₗₗ,ₜ(r)) [MW/m²]
 
         Notes
         -----
-        - The Eich target heat flux profile is derived from the midplane exponential
+        - The Eich parallel target heat flux profile is derived from the midplane exponential
         profile, taking into account the magnetic geometry and flux expansion between
         the midplane and the divertor target. The profile is typically characterized by
         a combination of an exponential decay and a Gaussian spreading due to cross-field
@@ -511,7 +511,7 @@ def calculate_eich_target_heat_flux_profile(
                 * f_b_div_flux_expansion
                 / (len_plasma_sol_power_spreading)
             )
-        ) + plux_target_background_heat_flux_mw
+        ) + pflux_target_background_heat_flux_mw
 
     @staticmethod
     def calculate_scarabosio2014_power_spreading_factor(

From 9e4515f44da20f8988eb3b4a7fd8dce26f1fad7b Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 26 Aug 2026 09:24:15 +0100
Subject: [PATCH 13/18] Add upstream radial decay and Eich heat flux profile
 sections to plasma scrape-off layer documentation

---
 .../physics-models/plasma_scrape_off_layer.md | 28 ++++++++++++++++++-
 1 file changed, 27 insertions(+), 1 deletion(-)

diff --git a/documentation/source/physics-models/plasma_scrape_off_layer.md b/documentation/source/physics-models/plasma_scrape_off_layer.md
index 5fc7318d8c..f6acf14979 100644
--- a/documentation/source/physics-models/plasma_scrape_off_layer.md
+++ b/documentation/source/physics-models/plasma_scrape_off_layer.md
@@ -25,7 +25,31 @@ $$
 A_{\parallel,u} = 2\pi\lambda_{\text{q,u}}R_{\text{u}}\frac{B_{\text{p,u}}}{B_{\text{Tot,u}}}
 $$
 
+---------------
 
+## Upstream radial decay | `calculate_outboard_midplane_near_sol_radial_profile()`
+
+The radial decay length $\lambda_{\text{q}}$ of the scrape-off layer (SOL) at the outer midplane of a tokamak is defined as the e-folding distance over which plasma heat and particle fluxes decay exponentially outside the last closed flux surface. Therefore the decay of the total heat flux outside the separatrix towards the vessel walls can be modelled as [^eich_2011] [^eich_2013]:
+
+$$
+q_{\text{u}}(r) = q_{\parallel,\text{u}}e^{\frac{-r}{\lambda_{\text{q}}}}
+$$
+
+where $r = R - R_{\text{sep}}$, $R_{\text{sep}}$ being the major radius of the separatrix, $\lambda_{\text{q}}$ the [power decay length](#power-decay-lengths) and $q_{\parallel}$ the [upstream energy flux density](#upstream-radial-decay--calculate_outboard_midplane_near_sol_radial_profile)
+
+----------------
+
+## Eich parallel flux at target | `calculate_eich_target_heat_flux_profile()`
+
+The Eich formula (often called the standard SOL heat flux profile) is the primary mathematical model used to describe the distribution of heat target loads on tokamak divertor plates. It convolutionally connects the physics of the plasma edge at the outer midplane with the geometric projection of the heat hitting the divertor surface [^eich_2011] [^eich_2013].
+
+Heat transport into the private flux region is modeled by convolving the power profile $q_{\text{u}}(r)$ with a Gaussian function of width $S$ known as the [spreading parameter](#spreading-parameter). 
+
+$$
+q_{\parallel,t} = \frac{q_0}{2}\times \exp\left(\left(\frac{S}{2\lambda_{\text{q}}}\right)- \frac{\overline{s}}{\lambda_q f_x}\right) \times \operatorname{erfc}\left(\frac{S}{2\lambda_{\text{q}}}- \frac{\overline{s}}{S f_{x}}\right) + q_{\text{BG}}
+$$
+
+where $\overline{s} = s- s_0 = (R_{\text{sep}} - R) \times f_x $. $\operatorname{erfc}$ is the complementary error function, $q_{\text{BG}}$ is the background heat flux, $\lambda_{\text{q}}$ is the [power decay length](#power-decay-lengths), $f_x$ is the effective flux expansion in the region,
 
 ------------------
 
@@ -119,4 +143,6 @@ Plasma Physics and Controlled Fusion, vol. 56, no. 5, p. 055008, Apr. 2014, doi:
 
 [^henderson_step]: S. S. Henderson et al., “An overview of the STEP divertor design and the simple models driving the plasma exhaust scenario,” Nuclear Fusion, vol. 65, no. 1, pp. 016033–016033, Nov. 2024, doi: 10.1088/1741-4326/ad93e7.
 
-[^scarabosio_2015]: A. Scarabosio et al., “Scaling of the divertor power spreading (S-factor) in open and closed divertor operation in JET and ASDEX Upgrade,” Journal of Nuclear Materials, vol. 463, pp. 49-54, Aug. 2015, doi: 10.1016/j.jnucmat.2014.11.076.
\ No newline at end of file
+[^scarabosio_2015]: A. Scarabosio et al., “Scaling of the divertor power spreading (S-factor) in open and closed divertor operation in JET and ASDEX Upgrade,” Journal of Nuclear Materials, vol. 463, pp. 49-54, Aug. 2015, doi: 10.1016/j.jnucmat.2014.11.076.
+
+[^eich_2011]: T. Eich, B. Sieglin, A. Scarabosio, W. Fundamenski, R. J. Goldston, and A. Herrmann, “Inter-ELM Power Decay Length for JET and ASDEX Upgrade: Measurement and Comparison with Heuristic Drift-Based Model,” Physical Review Letters, vol. 107, no. 21, Nov. 2011, doi: https://doi.org/10.1103/PhysRevLett.107.215001
\ No newline at end of file

From c311de822b9d736474e7ca74b9b1c0b3a622b686 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 26 Aug 2026 09:33:39 +0100
Subject: [PATCH 14/18] Add tests for outboard midplane near SOL radial profile
 and Eich target heat flux profile

---
 process/models/physics/scrape_off_layer.py    |  15 ++-
 .../models/physics/test_scrape_off_layer.py   | 101 ++++++++++++++++++
 2 files changed, 111 insertions(+), 5 deletions(-)

diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index cecd77244f..219e9cc71c 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -403,6 +403,11 @@ def calculate_outboard_midplane_near_sol_radial_profile(
         float|np.ndarray
             Outboard midplane SOL radial profile (qₗₗ(r)) [MW/m²]
 
+        Raises
+        ------
+        ValueError
+            If any radial position r is inside the plasma edge (r < rmajor + rminor)
+
         Notes
         -----
         - The exponential model is highly valid in the "near-SOL" (typically the first
@@ -472,11 +477,11 @@ def calculate_eich_target_heat_flux_profile(
 
         Notes
         -----
-        - The Eich parallel target heat flux profile is derived from the midplane exponential
-        profile, taking into account the magnetic geometry and flux expansion between
-        the midplane and the divertor target. The profile is typically characterized by
-        a combination of an exponential decay and a Gaussian spreading due to cross-field
-        transport in the divertor leg.
+        - The Eich parallel target heat flux profile is derived from the midplane
+        exponential profile, taking into account the magnetic geometry and flux expansion
+        between the midplane and the divertor target. The profile is typically
+        characterized by a combination of an exponential decay and a Gaussian spreading
+        due to cross-field transport in the divertor leg.
 
         References
         ----------
diff --git a/tests/unit/models/physics/test_scrape_off_layer.py b/tests/unit/models/physics/test_scrape_off_layer.py
index bbcd78a553..476f74a3d2 100644
--- a/tests/unit/models/physics/test_scrape_off_layer.py
+++ b/tests/unit/models/physics/test_scrape_off_layer.py
@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 
 from process.models.physics.scrape_off_layer import ScrapeOffLayer
@@ -112,3 +113,103 @@ def test_calculate_upstream_sol_outboard_parallel_area_exact():
     )
     assert isinstance(result, float)
     assert pytest.approx(result) == 0.006283185307179587
+
+
+@pytest.mark.parametrize(
+    "r",
+    [
+        8.001,
+        8.01,
+        8.1,
+    ],
+)
+def test_calculate_outboard_midplane_near_sol_radial_profile(r):
+    """Test outboard midplane near SOL radial profile with various parameters."""
+    result = ScrapeOffLayer.calculate_outboard_midplane_near_sol_radial_profile(
+        rmajor=6.0,
+        rminor=2.0,
+        len_plasma_sol_power_decay=0.001,
+        pflux_plasma_outboard_sol_parallel_mw=10.0,
+        r=r,
+    )
+    assert isinstance(result, float)
+    assert result > 0
+
+
+def test_calculate_outboard_midplane_near_sol_radial_profile_exact():
+    """Test outboard midplane near SOL radial profile with exact value check."""
+    result = ScrapeOffLayer.calculate_outboard_midplane_near_sol_radial_profile(
+        rmajor=6.0,
+        rminor=2.0,
+        len_plasma_sol_power_decay=0.001,
+        pflux_plasma_outboard_sol_parallel_mw=10.0,
+        r=8.001,
+    )
+    assert isinstance(result, float)
+    assert pytest.approx(result) == 3.678794411714423
+
+
+def test_calculate_outboard_midplane_near_sol_radial_profile_array():
+    """Test outboard midplane near SOL radial profile with array input."""
+    r = np.array([8.001, 8.002, 8.003])
+    result = ScrapeOffLayer.calculate_outboard_midplane_near_sol_radial_profile(
+        rmajor=6.0,
+        rminor=2.0,
+        len_plasma_sol_power_decay=0.001,
+        pflux_plasma_outboard_sol_parallel_mw=10.0,
+        r=r,
+    )
+    assert isinstance(result, np.ndarray)
+    assert np.all(result > 0)
+
+
+def test_calculate_outboard_midplane_near_sol_radial_profile_invalid_r():
+    """Test outboard midplane near SOL radial profile raises for r inside plasma edge."""
+    with pytest.raises(ValueError, match=r"inside plasma edge|outside plasma"):
+        ScrapeOffLayer.calculate_outboard_midplane_near_sol_radial_profile(
+            rmajor=6.0,
+            rminor=2.0,
+            len_plasma_sol_power_decay=0.001,
+            pflux_plasma_outboard_sol_parallel_mw=10.0,
+            r=7.0,
+        )
+
+
+@pytest.mark.parametrize(
+    "r",
+    [
+        8.001,
+        8.01,
+        8.1,
+    ],
+)
+def test_calculate_eich_target_heat_flux_profile(r):
+    """Test Eich target heat flux profile with various parameters."""
+    result = ScrapeOffLayer.calculate_eich_target_heat_flux_profile(
+        rmajor=6.0,
+        rminor=2.0,
+        pflux_plasma_sol_parallel_mw=10.0,
+        len_plasma_sol_power_decay=0.001,
+        f_b_div_flux_expansion=2.0,
+        len_plasma_sol_power_spreading=0.001,
+        pflux_target_background_heat_flux_mw=0.01,
+        r=r,
+    )
+    assert isinstance(result, float)
+    assert result > 0
+
+
+def test_calculate_eich_target_heat_flux_profile_exact():
+    """Test Eich target heat flux profile with exact value check."""
+    result = ScrapeOffLayer.calculate_eich_target_heat_flux_profile(
+        rmajor=6.0,
+        rminor=2.0,
+        pflux_plasma_sol_parallel_mw=10.0,
+        len_plasma_sol_power_decay=0.001,
+        f_b_div_flux_expansion=2.0,
+        len_plasma_sol_power_spreading=0.001,
+        pflux_target_background_heat_flux_mw=0.01,
+        r=8.001,
+    )
+    assert isinstance(result, float)
+    assert pytest.approx(result) == 3.8999590240461988

From 76a8d9141f01034b5736237f53464695918c7af5 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 26 Aug 2026 15:35:01 +0100
Subject: [PATCH 15/18] Enhance midplane near SOL radial profile plot with
 colour scheme and plasma boundary visualization; update Eich target heat flux
 profile title and adjust label positions for clarity.

---
 process/core/io/plot/summary.py | 36 ++++++++++++++++++++++++---------
 1 file changed, 26 insertions(+), 10 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index d7b2c8284a..3cff87729a 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9350,13 +9350,15 @@ def make_bbox_props(power: float) -> dict[str, Any]:
     axis.get_yaxis().set_ticks([])
 
 
-def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: int):
+def plot_midplane_near_sol_radial_profile(
+    axis: plt.Axes, mfile: MFile, scan: int, colour_scheme: int
+):
     """Function to plot the radial profile of the near SOL at the midplane."""
     rmajor = mfile.get("rmajor", scan=scan)
     rminor = mfile.get("rminor", scan=scan)
     len_sol_outboard_power_decay = mfile.get("len_sol_outboard_power_decay", scan=scan)
     r = np.linspace(
-        (rmajor + rminor), (rmajor + rminor) + (3 * len_sol_outboard_power_decay), 100
+        (rmajor + rminor), (rmajor + rminor) + (7 * len_sol_outboard_power_decay), 100
     )
 
     radial_profile = (
@@ -9372,7 +9374,15 @@ def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: in
     )
 
     x_ref = rmajor + rminor + len_sol_outboard_power_decay
+    plasma_boundary = rmajor + rminor
 
+    axis.axvspan(
+        0,
+        plasma_boundary,
+        color=PLASMA_COLOUR[colour_scheme - 1],
+        alpha=0.35,
+        label="Plasma",
+    )
     axis.plot(r, radial_profile, label=r"$q_{||}$ profile")
     axis.axvline(
         x_ref,
@@ -9383,9 +9393,13 @@ def plot_midplane_near_sol_radial_profile(axis: plt.Axes, mfile: MFile, scan: in
     )
     axis.grid()
     axis.legend()
+    axis.set_xlim(
+        (rmajor + rminor - (3 * len_sol_outboard_power_decay)),
+        (rmajor + rminor) + (7 * len_sol_outboard_power_decay),
+    )
     axis.set_title(r"Midplane Near SOL Radial Profile")
-    axis.set_xlabel("Radial Position [m]")
     axis.minorticks_on()
+    axis.tick_params(axis="x", labelbottom=False)
     axis.set_ylabel(r"$q_{||}$ [MW/m$^2$]")
 
 
@@ -9426,18 +9440,18 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     axis.axvline(peak_r, color="black", linestyle="--", linewidth=1)
     axis.axhline(peak_q, color="black", linestyle="--", linewidth=1)
     axis.text(
-        0.02,
-        0.98,
+        0.6,
+        0.9,
         f"$f_x$ = {f_b_flux_expansion:.2f}\n$S$ = {len_div_outboard_lower_power_spreading * 1e3:.3f} mm",
         transform=axis.transAxes,
         ha="left",
         va="top",
-        bbox={"boxstyle": "round", "facecolor": "white", "alpha": 0.8},
+        bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0},
     )
     axis.grid()
     axis.legend()
     axis.minorticks_on()
-    axis.set_title(r"Lower Outboard Eich Target Heat Flux Profile")
+    axis.set_title(r"Lower Outboard Eich Target Parallel Heat Flux Profile")
     axis.set_xlabel("Radial Position [m]")
     axis.set_ylabel(r"$q_{||,t}$ [MW/m$^2$]")
 
@@ -16875,13 +16889,15 @@ def _add_page(name: str | None = None):
         _add_page("sol_powerfluxes").add_subplot(121), m_file, scan, colour_scheme
     )
 
+    ax_midplane_near_sol = pages["sol_powerfluxes"].add_subplot(336)
     plot_midplane_near_sol_radial_profile(
-        pages["sol_powerfluxes"].add_subplot(336), m_file, scan
+        ax_midplane_near_sol, m_file, scan, colour_scheme
     )
 
-    plot_div_lower_outboard_eich_target_profile(
-        pages["sol_powerfluxes"].add_subplot(339), m_file, scan
+    ax_div_lower_outboard = pages["sol_powerfluxes"].add_subplot(
+        339, sharex=ax_midplane_near_sol
     )
+    plot_div_lower_outboard_eich_target_profile(ax_div_lower_outboard, m_file, scan)
 
     plot_debye_length_profile(
         _add_page("microscopic_quantities").add_subplot(232), m_file, scan

From ae3176256ee7ac3dc8241b7b7478cdbea89c314f Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 26 Aug 2026 15:43:57 +0100
Subject: [PATCH 16/18] Update flux expansion factor retrieval in
 plot_div_lower_outboard_eich_target_profile and enhance output logging in
 ScrapeOffLayer model

---
 process/core/io/plot/summary.py            |  2 +-
 process/models/physics/scrape_off_layer.py | 10 ++++++++++
 2 files changed, 11 insertions(+), 1 deletion(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 3cff87729a..6a517b086d 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9413,7 +9413,7 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     len_div_outboard_lower_power_spreading = mfile.get(
         "len_div_outboard_lower_power_spreading", scan=scan
     )
-    f_b_flux_expansion = 5.0
+    f_b_flux_expansion = mfile.get("f_b_div_outboard_lower_flux_expansion", scan=scan)
     r = np.linspace(
         (rmajor + rminor) - (1.5 * len_plasma_sol_power_decay) * f_b_flux_expansion,
         (rmajor + rminor) + (3 * len_plasma_sol_power_decay) * f_b_flux_expansion,
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 219e9cc71c..30bdc0f691 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -217,6 +217,16 @@ def output(self) -> None:
             "(len_div_outboard_lower_scrabosio14_power_spreading)",
             self.data.physics.len_div_outboard_lower_scrabosio14_power_spreading,
         )
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
+        po.oblnkl(self.outfile)
+        po.ovarre(
+            self.outfile,
+            "Outboard lower divertor flux expansion factor for the divertor targets "
+            "(fₓ)",
+            "(f_b_div_outboard_lower_flux_expansion)",
+            self.data.physics.f_b_div_outboard_lower_flux_expansion,
+        )
 
     @staticmethod
     def calculate_eich2013_sol_power_decay_length(

From 4bf39e844936481bc028ef50dcee4721df696a30 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 26 Aug 2026 16:29:30 +0100
Subject: [PATCH 17/18] :bug: Fix term placement in Eich profile

---
 .../physics-models/plasma_scrape_off_layer.md    | 16 ++++++++++++++--
 process/core/io/plot/summary.py                  |  7 ++++---
 process/models/physics/scrape_off_layer.py       |  9 ++-------
 .../unit/models/physics/test_scrape_off_layer.py | 14 +-------------
 4 files changed, 21 insertions(+), 25 deletions(-)

diff --git a/documentation/source/physics-models/plasma_scrape_off_layer.md b/documentation/source/physics-models/plasma_scrape_off_layer.md
index f6acf14979..b10f2780b6 100644
--- a/documentation/source/physics-models/plasma_scrape_off_layer.md
+++ b/documentation/source/physics-models/plasma_scrape_off_layer.md
@@ -46,10 +46,22 @@ The Eich formula (often called the standard SOL heat flux profile) is the primar
 Heat transport into the private flux region is modeled by convolving the power profile $q_{\text{u}}(r)$ with a Gaussian function of width $S$ known as the [spreading parameter](#spreading-parameter). 
 
 $$
-q_{\parallel,t} = \frac{q_0}{2}\times \exp\left(\left(\frac{S}{2\lambda_{\text{q}}}\right)- \frac{\overline{s}}{\lambda_q f_x}\right) \times \operatorname{erfc}\left(\frac{S}{2\lambda_{\text{q}}}- \frac{\overline{s}}{S f_{x}}\right) + q_{\text{BG}}
+q_{\parallel,t}(s) = \frac{q_0}{2}\times \exp\left[\left(\frac{S}{2\lambda_{\text{q}}f_x}\right)^2- \frac{s-s_0}{\lambda_q f_x}\right] \times \operatorname{erfc}\left(\frac{S}{2\lambda_{\text{q}}f_x}- \frac{s-s_0}{S}\right) + q_{\text{BG}}
 $$
 
-where $\overline{s} = s- s_0 = (R_{\text{sep}} - R) \times f_x $. $\operatorname{erfc}$ is the complementary error function, $q_{\text{BG}}$ is the background heat flux, $\lambda_{\text{q}}$ is the [power decay length](#power-decay-lengths), $f_x$ is the effective flux expansion in the region,
+where $s$ is the coordinate along the divertor target, $s_0$ is the strike-point location on the target, $\operatorname{erfc}$ is the complementary error function, $q_{\text{BG}}$ is the background heat flux, $\lambda_{\text{q}}$ is the [power decay length](#power-decay-lengths), $f_x$ is the effective flux expansion in the region,
+
+A compact equivalent form is:
+
+$$
+q_{\parallel,t}(\overline{s}) = \frac{q_0}{2}\times \exp\left[\left(\frac{S}{2\lambda_{\text{q}}f_x}\right)^2- \frac{\overline{s}}{\lambda_q f_x}\right] \times \operatorname{erfc}\left(\frac{S}{2\lambda_{\text{q}}f_x}- \frac{\overline{s}}{S}\right) + q_{\text{BG}}
+$$
+
+The connection to upstream midplane coordinates is usually:
+
+$$
+\overline{s} = f_x(R-R_{\text{sep}})
+$$
 
 ------------------
 
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 6a517b086d..7e25108f51 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9415,7 +9415,8 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     )
     f_b_flux_expansion = mfile.get("f_b_div_outboard_lower_flux_expansion", scan=scan)
     r = np.linspace(
-        (rmajor + rminor) - (1.5 * len_plasma_sol_power_decay) * f_b_flux_expansion,
+        (rmajor + rminor)
+        - ((f_b_flux_expansion / 2) * len_plasma_sol_power_decay) * f_b_flux_expansion,
         (rmajor + rminor) + (3 * len_plasma_sol_power_decay) * f_b_flux_expansion,
         200,
     )
@@ -9440,7 +9441,7 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
     axis.axvline(peak_r, color="black", linestyle="--", linewidth=1)
     axis.axhline(peak_q, color="black", linestyle="--", linewidth=1)
     axis.text(
-        0.6,
+        0.8,
         0.9,
         f"$f_x$ = {f_b_flux_expansion:.2f}\n$S$ = {len_div_outboard_lower_power_spreading * 1e3:.3f} mm",
         transform=axis.transAxes,
@@ -9449,10 +9450,10 @@ def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, sc
         bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0},
     )
     axis.grid()
-    axis.legend()
     axis.minorticks_on()
     axis.set_title(r"Lower Outboard Eich Target Parallel Heat Flux Profile")
     axis.set_xlabel("Radial Position [m]")
+    axis.set_xlim(r[0], r[-1])
     axis.set_ylabel(r"$q_{||,t}$ [MW/m$^2$]")
 
 
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 30bdc0f691..8cd5ef980f 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -513,19 +513,14 @@ def calculate_eich_target_heat_flux_profile(
             ** 2
             - (
                 (r - (rmajor + rminor))
-                * f_b_div_flux_expansion
-                / (len_plasma_sol_power_spreading * f_b_div_flux_expansion)
+                / (len_plasma_sol_power_decay * f_b_div_flux_expansion)
             )
         ) * scipy.special.erfc(
             (
                 len_plasma_sol_power_spreading
                 / (2 * len_plasma_sol_power_decay * f_b_div_flux_expansion)
             )
-            - (
-                (r - (rmajor + rminor))
-                * f_b_div_flux_expansion
-                / (len_plasma_sol_power_spreading)
-            )
+            - ((r - (rmajor + rminor)) / (len_plasma_sol_power_spreading))
         ) + pflux_target_background_heat_flux_mw
 
     @staticmethod
diff --git a/tests/unit/models/physics/test_scrape_off_layer.py b/tests/unit/models/physics/test_scrape_off_layer.py
index 476f74a3d2..7f7c843a18 100644
--- a/tests/unit/models/physics/test_scrape_off_layer.py
+++ b/tests/unit/models/physics/test_scrape_off_layer.py
@@ -163,18 +163,6 @@ def test_calculate_outboard_midplane_near_sol_radial_profile_array():
     assert np.all(result > 0)
 
 
-def test_calculate_outboard_midplane_near_sol_radial_profile_invalid_r():
-    """Test outboard midplane near SOL radial profile raises for r inside plasma edge."""
-    with pytest.raises(ValueError, match=r"inside plasma edge|outside plasma"):
-        ScrapeOffLayer.calculate_outboard_midplane_near_sol_radial_profile(
-            rmajor=6.0,
-            rminor=2.0,
-            len_plasma_sol_power_decay=0.001,
-            pflux_plasma_outboard_sol_parallel_mw=10.0,
-            r=7.0,
-        )
-
-
 @pytest.mark.parametrize(
     "r",
     [
@@ -212,4 +200,4 @@ def test_calculate_eich_target_heat_flux_profile_exact():
         r=8.001,
     )
     assert isinstance(result, float)
-    assert pytest.approx(result) == 3.8999590240461988
+    assert pytest.approx(result) == 5.534025566786268

From a47f2a2f8f0d9d0569fe62199066fa0781850778 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 27 Aug 2026 10:51:41 +0100
Subject: [PATCH 18/18] Update power spreading factor references from Scrabosio
 2014 to Scarabosio 2015 in physics variables and scrape-off layer model

---
 process/data_structure/physics_variables.py |  4 ++--
 process/models/physics/scrape_off_layer.py  | 16 ++++++++--------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 1f0a5608b7..f41d4f0241 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1764,8 +1764,8 @@ class PhysicsData:
     f_b_div_outboard_lower_flux_expansion: float = 5.0
     """Outboard lower divertor flux expansion factor for the divertor targets (fₓ)"""
 
-    len_div_outboard_lower_scrabosio14_power_spreading: float = 0.0
-    """Scrabosio 2014 H-mode power spreading length/factor in the scrape-off layer scaling (S) [m]"""
+    len_div_outboard_lower_scarabosio15_power_spreading: float = 0.0
+    """Scarabosio 2015 H-mode power spreading length/factor in the scrape-off layer scaling (S) [m]"""
 
     len_div_outboard_lower_power_spreading: float = 0.0
     """Power spreading length/factor at the outboard lower divertor target (S) [m]"""
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 8cd5ef980f..554d20c307 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -104,7 +104,7 @@ def run(self):
             / self.data.physics.a_plasma_outboard_sol_eich13_parallel
         )
 
-        self.data.physics.len_div_outboard_lower_scrabosio14_power_spreading = self.calculate_scarabosio2014_power_spreading_factor(  # noqa: E501
+        self.data.physics.len_div_outboard_lower_scarabosio15_power_spreading = self.calculate_scarabosio2015_power_spreading_factor(  # noqa: E501
             p_plasma_separatrix_mw=self.data.physics.p_plasma_separatrix_mw,
             b_plasma_surface_poloidal_average=self.data.physics.b_plasma_surface_poloidal_average,
             nd_plasma_separatrix_electron_19=self.data.physics.nd_plasma_separatrix_electron
@@ -113,7 +113,7 @@ def run(self):
         )
 
         self.data.physics.len_div_outboard_lower_power_spreading = (
-            self.data.physics.len_div_outboard_lower_scrabosio14_power_spreading
+            self.data.physics.len_div_outboard_lower_scarabosio15_power_spreading
         )
 
     def output(self) -> None:
@@ -213,9 +213,9 @@ def output(self) -> None:
         )
         po.ovarre(
             self.outfile,
-            "Scrabosio 2014 H-mode power spreading factor (S) [m]",
-            "(len_div_outboard_lower_scrabosio14_power_spreading)",
-            self.data.physics.len_div_outboard_lower_scrabosio14_power_spreading,
+            "Scarabosio 2015 H-mode power spreading factor (S) [m]",
+            "(len_div_outboard_lower_scarabosio15_power_spreading)",
+            self.data.physics.len_div_outboard_lower_scarabosio15_power_spreading,
         )
         po.oblnkl(self.outfile)
         po.ocmmnt(self.outfile, "----------------------------")
@@ -524,13 +524,13 @@ def calculate_eich_target_heat_flux_profile(
         ) + pflux_target_background_heat_flux_mw
 
     @staticmethod
-    def calculate_scarabosio2014_power_spreading_factor(
+    def calculate_scarabosio2015_power_spreading_factor(
         p_plasma_separatrix_mw: float,
         b_plasma_surface_poloidal_average: float,
         nd_plasma_separatrix_electron_19: float,
         rmajor: float,
     ) -> float:
-        """Calculate the Scrabosio 2014 H-mode power spreading factor (S).
+        """Calculate the Scarabosio 2015 H-mode power spreading factor (S).
 
         Parameters
         ----------
@@ -546,7 +546,7 @@ def calculate_scarabosio2014_power_spreading_factor(
         Returns
         -------
         float
-            Scrabosio 2014 H-mode power spreading factor (S) [m]
+            Scarabosio 2015 H-mode power spreading factor (S) [m]
 
         Notes
         -----