From a7e51a5a5c777ef8849a8736119b4a2be07344bf 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/19] 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 b7ea9ebf40..e51bf4c828 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -481,3 +481,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 62e2c8800211a36439d2dde16ad786efcb9c5af8 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/19] 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 a611ac63c2..644ee8d87c 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.superconductors import SuperconductorModel
 from process.models.tfcoil.base import (
     TFCoilShapeModel,
@@ -9357,6 +9358,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.
 
@@ -16739,6 +16784,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 e51bf4c828..c24bfa2224 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -530,7 +530,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 5947667f79064dd45407dc465584f231b6e0dcff 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/19] 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 644ee8d87c..95892e1db2 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -87,6 +87,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 c24bfa2224..9ee6bc3fd7 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
@@ -539,3 +540,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 96d1bec38d197ebb32c781c62ad4109646138d95 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/19] 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 95892e1db2..2fd3a3a9f8 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9383,24 +9383,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):
@@ -16786,7 +16815,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 9ee6bc3fd7..17be0c98c7 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -543,6 +543,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,
@@ -554,6 +556,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
@@ -598,11 +604,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 299c74282cc080aeb0f82b861b48591165fe7f13 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/19] 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 2fd3a3a9f8..69c19ba38e 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9431,6 +9431,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.
@@ -16814,14 +16855,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 4ec477a4189f59e4f002a7a2183084e2d02eb2e1 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/19] 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 69c19ba38e..a4e1b35343 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -85,7 +85,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
@@ -9369,6 +9368,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(
@@ -9383,11 +9383,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$]")
 
 
@@ -9400,7 +9411,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,
     )
@@ -9424,6 +9435,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()
@@ -9431,28 +9451,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",
@@ -9464,7 +9490,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)
@@ -9473,6 +9498,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.
 
@@ -16855,21 +16881,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 320e5dc968b9bc49cf511a6b9af1b876f001f7ca 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/19] 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 b09ce1a7cf..2cecd2d7c1 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1768,6 +1768,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 7b41314d3ff5d9548ae303507197ae304d6cada6 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/19] 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 2cecd2d7c1..02a1b5ccfd 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1768,7 +1768,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 17be0c98c7..d80227d3b3 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -130,6 +130,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 5e1706433c0febaa16a628a6b944a5fa4ecb725e 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/19] 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 02a1b5ccfd..977f85115f 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1772,6 +1772,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 3787086b892f37c6dfe3bd479a3fcf380da609f8 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/19] 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 977f85115f..dd2ca2f276 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1773,12 +1773,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 7c1fde407dd614838d43ab7b87bbeaba5e13d758 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/19] 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 d08badc3af..d5f0eb00c7 100644
--- a/documentation/source/physics-models/plasma_scrape_off_layer.md
+++ b/documentation/source/physics-models/plasma_scrape_off_layer.md
@@ -115,6 +115,33 @@ The scaling is done for type-I ELMy H-mode plasmas
 
 ------------------
 
+## 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,”
@@ -124,4 +151,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.
 
-[^eich_2011]: T. Eich, B. Sieglin, A. Scarabosio, W. Fundamenski, Robert James 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
+[^eich_2011]: T. Eich, B. Sieglin, A. Scarabosio, W. Fundamenski, Robert James 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. 
+
+[^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 d80227d3b3..04c09a6dd9 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -139,6 +139,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")
@@ -236,6 +248,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(
@@ -621,3 +649,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 947b712428bea713f20df33efab1e32407b76818 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/19] 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 a4e1b35343..47af5addb5 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9362,23 +9362,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,
         )
@@ -9409,6 +9404,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,
@@ -9424,8 +9422,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)
@@ -9438,7 +9436,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 04c09a6dd9..53f40cee23 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -578,10 +578,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
         ----------
@@ -597,7 +597,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]
@@ -605,11 +605,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
@@ -648,7 +648,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 1bd935ac2c90747b00ce732b097fd4de8803e797 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/19] 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 d5f0eb00c7..736e26e155 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,
 
 ------------------
 
@@ -153,4 +177,6 @@ Plasma Physics and Controlled Fusion, vol. 56, no. 5, p. 055008, Apr. 2014, doi:
 
 [^eich_2011]: T. Eich, B. Sieglin, A. Scarabosio, W. Fundamenski, Robert James 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. 
 
-[^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 b62a6110fdc9f73f00c0c0f259a83b6d70d4814d 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/19] 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 53f40cee23..3d5ad3efa0 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -540,6 +540,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
@@ -609,11 +614,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 e9015d55c8..2b90e5a677 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
@@ -189,3 +190,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 fd18d7ef82153a016557d6eb2ea8e4d9768d1748 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/19] 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 47af5addb5..5fe4987ed5 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9358,13 +9358,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 = (
@@ -9380,7 +9382,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,
@@ -9391,9 +9401,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$]")
 
 
@@ -9434,18 +9448,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$]")
 
@@ -16883,13 +16897,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 4ae4621429c211fd59ba94713fef83f63ea17db0 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/19] 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 5fe4987ed5..e7d66c98ad 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9421,7 +9421,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 3d5ad3efa0..3d20dfcd8e 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -264,6 +264,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 1a734c58fe0c1d5bb0ae542993cde0dbaa8132a4 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/19] :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 736e26e155..311d3dff90 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 e7d66c98ad..910849a31c 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -9423,7 +9423,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,
     )
@@ -9448,7 +9449,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,
@@ -9457,10 +9458,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 3d20dfcd8e..d9f67caec9 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -650,19 +650,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 2b90e5a677..088c678f1a 100644
--- a/tests/unit/models/physics/test_scrape_off_layer.py
+++ b/tests/unit/models/physics/test_scrape_off_layer.py
@@ -240,18 +240,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",
     [
@@ -289,4 +277,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 b98c327bc80f666cd7d46ecbac2891520855e11a 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/19] 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 dd2ca2f276..769cb15490 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1772,8 +1772,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 d9f67caec9..7db3f1cc08 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -139,7 +139,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
@@ -148,7 +148,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:
@@ -260,9 +260,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, "----------------------------")
@@ -661,13 +661,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
         ----------
@@ -683,7 +683,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
         -----

From a6ba7fb0214e686793d1b707be579175af53ba0a Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Fri, 11 Sep 2026 17:21:26 +0100
Subject: [PATCH 19/19] Update H-mode SOL spreading factor equation to correct
 exponent for n_sep

---
 documentation/source/physics-models/plasma_scrape_off_layer.md | 2 +-
 1 file changed, 1 insertion(+), 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 311d3dff90..4247d7359b 100644
--- a/documentation/source/physics-models/plasma_scrape_off_layer.md
+++ b/documentation/source/physics-models/plasma_scrape_off_layer.md
@@ -170,7 +170,7 @@ Unlike $\lambda_{q}$, which is governed by robust upstream parallel and perpendi
 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)}
+S = \left(0.12(\pm0.07)\times 10^{-3}\right) 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.02(\pm0.23)}
 $$
 
 - This was fitted from ASDEX Upgrade and JET outer target data