diff --git a/documentation/source/physics-models/plasma_exhaust.md b/documentation/source/physics-models/plasma_exhaust.md
new file mode 100644
index 0000000000..1d62172763
--- /dev/null
+++ b/documentation/source/physics-models/plasma_exhaust.md
@@ -0,0 +1,108 @@
+# Plasma Exhaust | `PlasmaExhaust`
+
+## Separatrix Power | `calculate_separatrix_power()`
+
+The separatrix power is defined as the thermal and conducted power flowing outward from the core plasma across the separatrix into the scrape-off layer.
+
+$$
+\overbrace{P_{\text{sep}}}^{\texttt{p_plasma_separatrix_mw}} = \underbrace{f_{\alpha}P_{\alpha} + P_{\text{c}} + P_{\text{OH}} + P_{\text{HCD}}}_{\text{Plasma heating}} - \underbrace{P_{\text{rad}}}_{\text{Plasma Loss}}
+$$
+
+where $f_{\alpha}$ is the [fraction of alpha power that is coupled to the plasma](../physics-models/fusion_reactions/plasma_reactions.md#coupled-alpha-particle-power), $P_{\alpha}$ is the alpha power, $P_{\text{c}}$ is the charged particle power, $P_{\text{OH}}$ is the ohmic heating power, $P_{\text{HCD}}$ is the plasma heating done by the external heating & current drive systems and $P_{\text{rad}}$ is the total radiation power given off by the plasma.
+
+--------------
+
+## Divertor Protection Metrics
+
+### `calculate_psep_over_r_metric()`
+
+A simple metric for the divertor heat flux challenege is to simply divide by total separatrix power by the plasma major radius. This works as a stand in metric for the total available divertor area in which to spread the heat
+
+
+$$
+\texttt{p_plasma_separatrix_rmajor_mw}=\frac{P_{\text{sep}}}{R_0}
+$$
+
+------------
+
+### EU-DEMO Protection Re-Attachment Metric | `calculate_eu_demo_re_attachment_metric()`
+
+The divertor protection metric is obtained by starting from the target heat flux formula, inserting the Eich heat-flux-width scaling and $q_{95}-B_{\text{p}}$​ relationship, and reducing the resulting dependence to the engineering quantity whose value is constrained to remain below a reference EU-DEMO limit to ensure the divertor can withstand temporary reattachment without[^1][^2].
+
+$$
+\texttt{p_div_bt_q_aspect_rmajor_mw}=\frac{P_{\text{sep}}B_{\text{T}}}{q_{95}AR}
+$$
+
+
+-----------
+
+## Divertor Separatrix Power Heat Splits
+
+### Brunner Model | `calculate_brunner_divertor_power_splits()`
+
+The Brunner model[^brunner_model] was made to estimate the power-sharing between four divertors in a dynamic double null configuration. Consisting of an in-out and up-down sharing scheme. The model considers the separatrix power flow going to each divertor rather than that which reaches the target plate after dissipation. The power split assumes that the fraction of power to the side walls is negligibly small and there is no radiation in the SOL. We thus have:
+
+$$
+P_{\text{sep}} \approx P_{\text{div,tot}} = P_{\text{lower,in}} + P_{\text{lower,out}} + P_{\text{upper,in}} + P_{\text{upper,out}}
+$$
+
+The main determining variable for the power split is $\delta R_{\text{sep}}$ representing the radial distance between the first and second separatrixes at the plasma outboard mid-plane.
+
+A negative value of $\delta R_{\text{sep}}$ represents a dominant lower null configuration, A $\delta R_{\text{sep}}$ of 0 represents a perfect double null configuration and a positive $\delta R_{\text{sep}}$ is a upper null configuration.
+
+It is observed that the majority of the power crosses the plasma boundary at the outer mid-plane. Assuming that	the	resultant heat flux profile	is exponential with	a single e-folding width equal to the scrape off layer power decay length, $\lambda_{\text{q,o}}$.
+
+For the outboard, side the fraction of the outboard power going to the lower target is:
+
+$$
+\frac{P_{\text{lower,out}}}{P_{\text{lower,out}}+P_{\text{upper,out}}} = \frac{1}{1+e^{\frac{\delta R_{\text{sep}}}{\lambda_{\text{q,o}}}}}
+$$
+
+For the upper target:
+
+$$
+\frac{P_{\text{upper,out}}}{P_{\text{lower,out}}+P_{\text{upper,out}}} = \frac{1}{1+e^{-\frac{\delta R_{\text{sep}}}{\lambda_{\text{q,o}}}}}
+$$
+
+For the inboard divertors it was seen that the power sharing follows a similar logistic function trend but with a distinctly different power decay length, $\lambda_{\text{q,i}}$.
+
+For the inboard side the fraction of the inboard power going to the lower target is:
+
+$$
+\frac{P_{\text{lower,in}}}{P_{\text{lower,in}}+P_{\text{upper,in}}} = \frac{1}{1+e^{\frac{\delta R_{\text{sep}}}{\lambda_{\text{q,i}}}}}
+$$
+
+The total inboard split compared to the outboard was found to give a Gaussian like dependence as follows:
+
+$$
+\frac{P_{\text{lower,in}}+P_{\text{lower,out}}}{P_{\text{lower,in}}+P_{\text{upper,in}}+P_{\text{lower,out}}+P_{\text{upper,out}}} = \\ P_{\text{in,0}}+(P_{\text{in,0}}-P_{\text{in,}\infin})\times \left(1-\frac{2}{1+e^{-\left(\frac{\delta R_{\text{sep}}}{\lambda_{\text{q,io}}}\right)^2}}\right)
+$$
+
+where $P_{\text{in,0}}$ is the fraction of power to the inner divertors at $\delta R_{\text{sep}} = 0$, $P_{\text{in,}\infin}$ is the fraction of the power to the inner divertors at $R_{\text{sep}} = \infin$
+
+Fitted values to data from DIII-D[^d3_d_values] have shown the values of $P_{\text{in,0}}$ and $P_{\text{in,}\infin}$ to be 0.16 and 0.41 respectively which is what we use here. We also assume that out Gaussian width in the total power sharing is equal to the outboard power decay length, $\lambda_{\text{q,io}} = \lambda_{\text{q,o}}$.
+
+In `PROCESS` we normally want to control how much of the separatrix power is directed to either the upper or lower divertors, this is done through changing `f_p_div_lower_separatrix`. Fitting to the DIII-D[^d3_d_values] we calculate the value of $\delta R_{\text{sep}}$ via:
+
+$$
+\delta R_{\text{sep}} = -3\times 10^{-3} \tanh^{-1}\left(2(\texttt{f_p_div_lower_separatrix}-0.5)\right)
+$$
+
+Truncated values of the $\delta R_{\text{sep}}$ are set to $1.5\times 10^{-2}\ \text{m}$ and $-1.5\times 10^{-2} \ \text{m}$ respectively based on when the power split gets close to a perfect upper or lower null.
+
+--------------------
+
+
+[^1] M. Siccinio, G. Federici, R. Kembleton, H. Lux, F. Maviglia, and J. Morris,
+"Figure of merit for divertor protection in the preliminary design of the
+EU-DEMO reactor," Nuclear Fusion, vol. 59, no. 10, pp. 106026-106026,
+Jul. 2019, doi: https://doi.org/10.1088/1741-4326/ab3153.
+
+[^2] H. Zohm et al., "A stepladder approach to a tokamak fusion power plant,"
+Nuclear Fusion, vol. 57, no. 8, pp. 086002-086002, May 2017, doi: https://doi.org/10.1088/1741-4326/aa739e.
+
+[^brunner_model] D. Brunner, A. Q. Kuang, B. LaBombard, and J. L. Terry, “The dependence of
+divertor power sharing on magnetic flux balance in near double-null configurations on 
+Alcator C-Mod,” Nuclear Fusion, vol. 58, no. 7, p. 076010, May 2018, doi: https://doi.org/10.1088/1741-4326/aac006.
+
+[^d3_d_values] T. W. Petrie et al., “The effect of divertor magnetic balance on H-mode performance in DIII-D,” Journal of Nuclear Materials, vol. 290-293, pp. 935-939, Mar. 2001, doi: https://doi.org/10.1016/S0022-3115(00)00492-X
\ No newline at end of file
diff --git a/mkdocs.yml b/mkdocs.yml
index 1c80b5d1de..6cfa762788 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -67,6 +67,7 @@ nav:
               - Confinement time: physics-models/plasma_confinement.md
               - L-H transition: physics-models/plasma_h_mode.md
               - Plasma Core Power Balance: physics-models/plasma_power_balance.md
+              - Plasma Exhaust: physics-models/plasma_exhaust.md
               - Plasma Scrape-off Layer: physics-models/plasma_scrape_off_layer.md
               - Detailed Plasma Physics: physics-models/detailed_physics.md
           - Pulsed Plant Operation: physics-models/pulsed-plant.md
diff --git a/process/core/init.py b/process/core/init.py
index 7704fd05e8..a1be430db5 100644
--- a/process/core/init.py
+++ b/process/core/init.py
@@ -718,7 +718,7 @@ def check_process(inputs, data):  # noqa: ARG001
 
         # Check if a single null divertor is used in double null machine
         if i_single_null == DivertorNumberModels.DOUBLE_NULL and (
-            data.physics.f_p_div_lower in {1.0, 0.0}
+            data.physics.f_p_div_lower_separatrix in {1.0, 0.0}
         ):
             logger.error(
                 "Operating with a single null in a double null machine", stacklevel=2
diff --git a/process/core/input.py b/process/core/input.py
index bf30712004..600d953816 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -186,7 +186,7 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "beta_norm_max": InputVariable("physics", float, range=(0.0, 20.0)),
     "beta_poloidal_eps_max": InputVariable("physics", float, range=(0.01, 10.0)),
     "f_p_alpha_plasma_deposited": InputVariable("physics", float, range=(0.0, 1.0)),
-    "f_p_div_lower": InputVariable("physics", float, range=(0.0, 1.0)),
+    "f_p_div_lower_separatrix": InputVariable("physics", float, range=(0.0, 1.0)),
     "f_plasma_fuel_deuterium": InputVariable("physics", float, range=(0.0, 1.0)),
     "ffwal": InputVariable("physics", float, range=(0.0, 10.0)),
     "f_nd_plasma_pedestal_greenwald": InputVariable("physics", float, range=(0.1, 1.5)),
@@ -1181,6 +1181,9 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     ),
     "i_len_sol_outboard_power_decay": InputVariable("physics", int, choices=range(4)),
     "len_sol_outboard_power_decay": InputVariable("physics", float, range=(0.0, 1.0)),
+    "f_len_sol_power_decay_inboard_outboard": InputVariable(
+        "physics", float, range=(0.01, 2.0)
+    ),
 }
 
 
diff --git a/process/core/io/obsolete_vars.py b/process/core/io/obsolete_vars.py
index 7227a4b773..a9734f8f13 100644
--- a/process/core/io/obsolete_vars.py
+++ b/process/core/io/obsolete_vars.py
@@ -301,7 +301,7 @@
     "irfcd": "i_hcd_calculations",
     "wave_mode": "i_ecrh_wave_mode",
     "pheatfix": "p_hcd_secondary_extra_heat_mw",
-    "ftar": "f_p_div_lower",
+    "ftar": "f_p_div_lower_separatrix",
     "beta_div": "deg_div_field_plate",
     "divclfr": "f_vol_div_coolant",
     "divdens": "den_div_structure",
@@ -478,6 +478,7 @@
     "minmax": "i_figure_merit",
     "neqns": "n_equality_constraints",
     "nineqns": "n_inequality_constraints",
+    "f_p_div_lower": "f_p_div_lower_separatrix",
 }
 
 OBS_VARS_HELP = {
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 5c53ac2933..887eea4f09 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -23,7 +23,10 @@
 from process.data_structure.impurity_radiation_variables import N_IMPURITIES
 from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode
 from process.data_structure.pfcoil_variables import NFIXMX
-from process.data_structure.physics_variables import ConfinementTimeModel
+from process.data_structure.physics_variables import (
+    ConfinementTimeModel,
+    DivertorNumberModels,
+)
 from process.data_structure.superconducting_tf_coil_variables import TFWPIntegerTurnType
 from process.models.build import Build
 from process.models.engineering.materials import (
@@ -64,6 +67,7 @@
     ElectronCyclotron,
 )
 from process.models.physics.density_limit import DensityLimitModel
+from process.models.physics.exhaust import calculate_brunner_divertor_power_splits
 from process.models.physics.impurity_radiation import read_impurity_file
 from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel
 from process.models.physics.physics import (
@@ -9087,6 +9091,263 @@ def plot_sol_power_decay_length_comparison(axis: plt.Axes, mfile: MFile, scan: i
     axis.set_facecolor("#f0f0f0")
 
 
+def plot_brunner_divertor_power_split_comparison_stackplot(
+    axis: plt.Axes, mfile: MFile, scan: int
+):
+    """Plot Brunner divertor power split fractions as a stack plot over dr_sep."""
+    # Use the case decay length when available; fall back to 1 mm if absent.
+
+    len_plasma_sol_outboard_pd = mfile.get("len_sol_outboard_power_decay", scan=scan)
+    len_plasma_sol_inboard_pd = mfile.get("len_sol_inboard_power_decay", scan=scan)
+    colors = plt.cm.plasma(np.linspace(0.15, 0.85, 4))
+
+    dr_sep_values = np.linspace(
+        -5 * len_plasma_sol_outboard_pd,
+        5 * len_plasma_sol_outboard_pd,
+        200,
+    )
+    f_p_inboard_lower = np.zeros_like(dr_sep_values)
+    f_p_inboard_upper = np.zeros_like(dr_sep_values)
+    f_p_outboard_lower = np.zeros_like(dr_sep_values)
+    f_p_outboard_upper = np.zeros_like(dr_sep_values)
+
+    for idx, dr_sep in enumerate(dr_sep_values):
+        div_power_splits = calculate_brunner_divertor_power_splits(
+            dr_outboard_midplane_sep=dr_sep,
+            len_plasma_sol_outboard_power_decay=len_plasma_sol_outboard_pd,
+            len_plasma_sol_inboard_power_decay=len_plasma_sol_inboard_pd,
+        )
+        f_p_inboard_lower[idx] = div_power_splits.f_p_div_inboard_lower_separatrix
+        f_p_inboard_upper[idx] = div_power_splits.f_p_div_inboard_upper_separatrix
+        f_p_outboard_lower[idx] = div_power_splits.f_p_div_outboard_lower_separatrix
+        f_p_outboard_upper[idx] = div_power_splits.f_p_div_outboard_upper_separatrix
+
+    axis.stackplot(
+        dr_sep_values,
+        f_p_inboard_lower,
+        f_p_inboard_upper,
+        f_p_outboard_lower,
+        f_p_outboard_upper,
+        labels=[
+            "$f_{P,\\mathrm{in,lower}}$",
+            "$f_{P,\\mathrm{in,upper}}$",
+            "$f_{P,\\mathrm{out,lower}}$",
+            "$f_{P,\\mathrm{out,upper}}$",
+        ],
+        colors=colors,
+        alpha=0.9,
+    )
+
+    axis.axvline(
+        mfile.get("dr_plasma_outboard_midplane_separatrix_separation", scan=scan),
+        color="k",
+        linestyle="--",
+        linewidth=1.0,
+        alpha=0.5,
+        label="$\u0394 r_{\\mathrm{sep}}$",
+    )
+    axis.set_ylim([0.0, 1.0])
+    axis.set_xlim([
+        -5 * len_plasma_sol_outboard_pd,
+        5 * len_plasma_sol_outboard_pd,
+    ])
+    axis.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.35)
+    axis.set_title("Brunner Divertor Power Split Fractions")
+    axis.set_xlabel("$\\Delta r_{\\mathrm{sep}}$ [m]")
+    axis.set_ylabel("Power split fraction, $f_P$")
+    axis.legend(loc="upper left", fontsize=8)
+
+
+def plot_separatrix_power_split(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, dr_sep = mfile.get_variables(
+        "rmajor",
+        "rminor",
+        "kappa",
+        "dr_plasma_outboard_midplane_separatrix_separation",
+        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
+
+    is_double_null = (
+        DivertorNumberModels(mfile.get("i_single_null", scan=scan))
+        == DivertorNumberModels.DOUBLE_NULL
+    )
+    p_sep = mfile.get("p_plasma_separatrix_mw", scan=scan)
+    f_outboard = mfile.get("f_p_div_outboard_separatrix", scan=scan)
+    f_inboard = mfile.get("f_p_div_inboard_separatrix", scan=scan)
+    p_outboard = p_sep * f_outboard
+    p_inboard = p_sep * f_inboard
+    p_lower_inboard = mfile.get("p_div_lower_inboard_separatrix_mw", scan=scan)
+    p_lower_outboard = mfile.get("p_div_lower_outboard_separatrix_mw", scan=scan)
+
+    power_values = [
+        p_sep,
+        p_outboard,
+        p_inboard,
+        p_lower_inboard,
+        p_lower_outboard,
+    ]
+
+    p_upper_inboard = None
+    p_upper_outboard = None
+    if is_double_null:
+        p_upper_inboard = mfile.get("p_div_upper_inboard_separatrix_mw", scan=scan)
+        p_upper_outboard = mfile.get("p_div_upper_outboard_separatrix_mw", scan=scan)
+        power_values.extend([p_upper_inboard, p_upper_outboard])
+
+    power_min = min(power_values)
+    power_max = max(power_values)
+    colour_map = mpl.colormaps["coolwarm"]
+
+    def make_bbox_props(power: float) -> dict[str, Any]:
+        norm_power = (
+            1.0
+            if np.isclose(power_max, power_min)
+            else (power - power_min) / (power_max - power_min)
+        )
+        return {
+            "boxstyle": f"round,pad={0.3 * scale_factor:.3f}",
+            "facecolor": colour_map(norm_power),
+            "alpha": 1.0,
+            "linewidth": 2 * scale_factor,
+            "edgecolor": "black",
+        }
+
+    centre_pos = (rmajor, 0.0)
+    outboard_pos = (rmajor + rminor, 0.0)
+    inboard_pos = (rmajor - rminor, 0.0)
+    lower_inboard_pos = (rmajor - rminor, -kappa * rminor)
+    lower_outboard_pos = (rmajor + rminor, -kappa * rminor)
+    upper_inboard_pos = (rmajor - rminor, kappa * rminor)
+    upper_outboard_pos = (rmajor + rminor, kappa * rminor)
+
+    axis.text(
+        *centre_pos,
+        f"$P_{{\\mathrm{{sep}}}} = {p_sep:.3f}$ MW",
+        fontsize=text_fontsize,
+        verticalalignment="center",
+        horizontalalignment="center",
+        bbox=make_bbox_props(p_sep),
+        zorder=101,
+    )
+    axis.text(
+        *outboard_pos,
+        f"$f_{{\\mathrm{{outboard}}}} = {f_outboard:.3f}$\n"
+        f"$\\Delta r_{{\\mathrm{{sep}}}} = {dr_sep:.3f}$ m",
+        fontsize=text_fontsize,
+        verticalalignment="center",
+        horizontalalignment="center",
+        bbox=make_bbox_props(p_outboard),
+        zorder=101,
+    )
+    axis.text(
+        *inboard_pos,
+        f"$f_{{\\mathrm{{inboard}}}} = {f_inboard:.3f}$",
+        fontsize=text_fontsize,
+        verticalalignment="center",
+        horizontalalignment="center",
+        bbox=make_bbox_props(p_inboard),
+        zorder=101,
+    )
+    axis.text(
+        *lower_inboard_pos,
+        f"$f_{{\\mathrm{{lower\\ inboard}}}} = {mfile.get('f_p_div_lower_inboard_separatrix', scan=scan):.3f}$\n"
+        f"$P_{{\\mathrm{{lower\\ inboard}}}} = {p_lower_inboard:.3f}$ MW",
+        fontsize=text_fontsize,
+        verticalalignment="center",
+        horizontalalignment="center",
+        bbox=make_bbox_props(p_lower_inboard),
+        zorder=101,
+    )
+    axis.text(
+        *lower_outboard_pos,
+        f"$f_{{\\mathrm{{lower\\ outboard}}}} = {mfile.get('f_p_div_lower_outboard_separatrix', scan=scan):.3f}$\n"
+        f"$P_{{\\mathrm{{lower\\ outboard}}}} = {p_lower_outboard:.3f}$ MW",
+        fontsize=text_fontsize,
+        verticalalignment="center",
+        horizontalalignment="center",
+        bbox=make_bbox_props(p_lower_outboard),
+        zorder=101,
+    )
+    if is_double_null:
+        axis.text(
+            *upper_inboard_pos,
+            f"$f_{{\\mathrm{{upper\\ inboard}}}} = {mfile.get('f_p_div_upper_inboard_separatrix', scan=scan):.3f}$\n"
+            f"$P_{{\\mathrm{{upper\\ inboard}}}} = {p_upper_inboard:.3f}$ MW",
+            fontsize=text_fontsize,
+            verticalalignment="center",
+            horizontalalignment="center",
+            bbox=make_bbox_props(p_upper_inboard),
+            zorder=101,
+        )
+        axis.text(
+            *upper_outboard_pos,
+            f"$f_{{\\mathrm{{upper\\ outboard}}}} = {mfile.get('f_p_div_upper_outboard_separatrix', scan=scan):.3f}$\n"
+            f"$P_{{\\mathrm{{upper\\ outboard}}}} = {p_upper_outboard:.3f}$ MW",
+            fontsize=text_fontsize,
+            verticalalignment="center",
+            horizontalalignment="center",
+            bbox=make_bbox_props(p_upper_outboard),
+            zorder=101,
+        )
+
+    arrow_props = {
+        "arrowstyle": "->",
+        "color": "red",
+        "linewidth": 3 * scale_factor,
+        "shrinkA": 14 * scale_factor,
+        "shrinkB": 14 * scale_factor,
+        "mutation_scale": 12 * scale_factor,
+    }
+    axis.annotate(
+        "", xy=outboard_pos, xytext=centre_pos, arrowprops=arrow_props, zorder=1
+    )
+    axis.annotate(
+        "", xy=inboard_pos, xytext=centre_pos, arrowprops=arrow_props, zorder=1
+    )
+    axis.annotate(
+        "",
+        xy=lower_outboard_pos,
+        xytext=outboard_pos,
+        arrowprops={**arrow_props, "connectionstyle": "angle3,angleA=0,angleB=-90"},
+        zorder=102,
+    )
+    axis.annotate(
+        "",
+        xy=lower_inboard_pos,
+        xytext=inboard_pos,
+        arrowprops={**arrow_props, "connectionstyle": "angle3,angleA=180,angleB=-90"},
+        zorder=102,
+    )
+    if is_double_null:
+        axis.annotate(
+            "",
+            xy=upper_outboard_pos,
+            xytext=outboard_pos,
+            arrowprops={**arrow_props, "connectionstyle": "angle3,angleA=0,angleB=90"},
+            zorder=102,
+        )
+        axis.annotate(
+            "",
+            xy=upper_inboard_pos,
+            xytext=inboard_pos,
+            arrowprops={**arrow_props, "connectionstyle": "angle3,angleA=180,angleB=90"},
+            zorder=102,
+        )
+
+    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.
 
@@ -16445,6 +16706,14 @@ def _add_page(name: str | None = None):
         _add_page("plasma_compare_3").add_subplot(221), m_file, scan
     )
 
+    plot_brunner_divertor_power_split_comparison_stackplot(
+        _add_page("plasma_exhaust").add_subplot(121), m_file, scan
+    )
+
+    plot_separatrix_power_split(
+        pages["plasma_exhaust"].add_subplot(122), m_file, scan, colour_scheme
+    )
+
     plot_debye_length_profile(
         _add_page("microscopic_quantities").add_subplot(232), m_file, scan
     )
diff --git a/process/core/scan.py b/process/core/scan.py
index d62a84a769..f7fa304e8e 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -191,7 +191,7 @@ def get_val(self, mfile, scan):
     n_tf_wp_pancakes = (48, Area.T)
     n_tf_wp_layers = (49, Area.T)
     f_nd_impurity_electrons__13 = (50, Area.IR)
-    f_p_div_lower = (51, Area.P)
+    f_p_div_lower_separatrix = (51, Area.P)
     rad_fraction_sol = (52, Area.P)
     boundu__157 = (53, Area.NUM)
     b_crit_upper_nbti = (54, Area.T)
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 910ce574c1..11539af97a 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -430,30 +430,49 @@ class PhysicsData:
     t_energy_confinement_beta: float = 0.0
     """[s]"""
 
-    ptarmw: float = 0.0
+    f_len_sol_power_decay_inboard_outboard: float = 0.7
+    """Ratio of inboard to outboard scrape off layer power decay length (λᵢₙ_q/λₒᵤₜ_q)"""
+
+    len_sol_inboard_power_decay: float = 0.0
+    """Inboard scrape off layer power decay length (λᵢₙ_q) [m]"""
 
     len_sol_outboard_power_decay: float = 0.0
-    """Outboard scrape off layer power decay length (λ_q) [m]"""
+    """Outboard scrape off layer power decay length (λₒᵤₜ_q) [m]"""
+
+    dr_plasma_outboard_midplane_separatrix_separation: float = 0.0
+    """Radial distance between the first and second plasma separatrixes at the outer midplane (δR_sep) [m]"""
 
-    drsep: float = 0.0
+    f_p_div_inboard_separatrix: float = 0.0
+    """Fraction of the total plasma separatrix power going to the inboard divertor
+    target(s)"""
 
-    fio: float = 0.0
+    f_p_div_outboard_separatrix: float = 0.0
+    """Fraction of the total plasma separatrix power going to the outboard divertor
+    target(s)"""
 
-    fli: float = 0.0
+    f_p_div_lower_inboard_separatrix: float = 0.0
+    """Fraction of plasma separatrix power to the lower inboard divertor target"""
 
-    flo: float = 0.0
+    f_p_div_lower_outboard_separatrix: float = 0.0
+    """Fraction of plasma separatrix power to the lower outboard divertor target"""
 
-    fui: float = 0.0
+    f_p_div_upper_inboard_separatrix: float = 0.0
+    """Fraction of plasma separatrix power to the upper inboard divertor target"""
 
-    fuo: float = 0.0
+    f_p_div_upper_outboard_separatrix: float = 0.0
+    """Fraction of plasma separatrix power to the upper outboard divertor target"""
 
-    plimw: float = 0.0
+    p_div_lower_inboard_separatrix_mw: float = 0.0
+    """Plasma separatrix power directed to the lower inboard divertor target [MW]"""
 
-    plomw: float = 0.0
+    p_div_lower_outboard_separatrix_mw: float = 0.0
+    """Plasma separatrix power directed to the lower outboard divertor target [MW]"""
 
-    puimw: float = 0.0
+    p_div_upper_inboard_separatrix_mw: float = 0.0
+    """Plasma separatrix power directed to the upper inboard divertor target [MW]"""
 
-    puomw: float = 0.0
+    p_div_upper_outboard_separatrix_mw: float = 0.0
+    """Plasma separatrix power directed to the upper outboard divertor target [MW]"""
 
     rho_star: float = 0.0
 
@@ -737,7 +756,7 @@ class PhysicsData:
     f_plasma_fuel_deuterium: float = 0.5
     """Plasma deuterium fuel fraction"""
 
-    f_p_div_lower: float = 1.0
+    f_p_div_lower_separatrix: float = 1.0
     """fraction of power to the lower divertor in double null configuration
     (`i_single_null = 0` only) (default assumes SN)
     """
diff --git a/process/models/build.py b/process/models/build.py
index d0f86fa049..de7167d7c2 100644
--- a/process/models/build.py
+++ b/process/models/build.py
@@ -1471,6 +1471,13 @@ def divgeom(self, output: bool):
                     self.outfile,
                     "ERROR: null value not supported, check i_single_null value.",
                 )
+            po.ovarre(
+                self.outfile,
+                "Divertor poloidal angle subtended by plasma (degrees)",
+                "(deg_div_poloidal_plasma)",
+                self.data.divertor.deg_div_poloidal_plasma,
+                "OP ",
+            )
         return divht
 
     @staticmethod
diff --git a/process/models/divertor.py b/process/models/divertor.py
index a9563e4383..2723bd8c20 100644
--- a/process/models/divertor.py
+++ b/process/models/divertor.py
@@ -98,7 +98,7 @@ def run(self, output: bool = False):
                 self.data.physics.nd_plasma_separatrix_electron,
                 self.data.divertor.deg_div_field_plate,
                 self.data.physics.rad_fraction_sol,
-                self.data.physics.f_p_div_lower,
+                self.data.physics.f_p_div_lower_separatrix,
                 output=output,
             )
             return
@@ -281,7 +281,7 @@ def divwade(
         nd_plasma_separatrix_electron: float,
         deg_div_field_plate: float,
         rad_fraction_sol: float,
-        f_p_div_lower: float,
+        f_p_div_lower_separatrix: float,
         output: bool,
     ) -> float:
         """Divertor heat load model (Wade 2020)
@@ -316,7 +316,7 @@ def divwade(
             field line angle wrt divertor target plate (degrees)
         rad_fraction_sol : float
             SOL radiation fraction
-        f_p_div_lower : float
+        f_p_div_lower_separatrix : float
             fraction of power to the lower divertor in double null configuration
 
         Returns
@@ -375,8 +375,8 @@ def divwade(
         # For double null, calculate heat loads to upper and lower divertors
         # and use the highest
         if self.data.divertor.n_divertors == 2:
-            hldiv_lower = f_p_div_lower * hldiv_base
-            hldiv_upper = (1.0 - f_p_div_lower) * hldiv_base
+            hldiv_lower = f_p_div_lower_separatrix * hldiv_base
+            hldiv_upper = (1.0 - f_p_div_lower_separatrix) * hldiv_base
             self.data.divertor.pflux_div_heat_load_mw = max(hldiv_lower, hldiv_upper)
         else:
             self.data.divertor.pflux_div_heat_load_mw = hldiv_base
diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py
index c61f9c66aa..95a12fee60 100644
--- a/process/models/physics/exhaust.py
+++ b/process/models/physics/exhaust.py
@@ -1,14 +1,41 @@
 """Module for plasma exhaust calculations and analysis."""
 
 import logging
+from dataclasses import dataclass
+
+import numpy as np
 
 from process.core import constants
 from process.core import process_output as po
 from process.core.model import Model
+from process.data_structure.physics_variables import DivertorNumberModels
 
 logger = logging.getLogger(__name__)
 
 
+@dataclass
+class DivertorSeparatrixPowerSplits:
+    """Dataclass to hold the power splits to the divertor targets."""
+
+    f_p_div_inboard_separatrix: float = 0.0
+    """Fraction of total separatrix power to inboard divertor targets"""
+
+    f_p_div_outboard_separatrix: float = 0.0
+    """Fraction of total separatrix power to outboard divertor targets"""
+
+    f_p_div_inboard_lower_separatrix: float = 0.0
+    """Fraction of total separatrix power to inboard lower divertor target"""
+
+    f_p_div_inboard_upper_separatrix: float = 0.0
+    """Fraction of total separatrix power to inboard upper divertor target"""
+
+    f_p_div_outboard_lower_separatrix: float = 0.0
+    """Fraction of total separatrix power to outboard lower divertor target"""
+
+    f_p_div_outboard_upper_separatrix: float = 0.0
+    """Fraction of total separatrix power to outboard upper divertor target"""
+
+
 class PlasmaExhaust(Model):
     """Class to hold plasma exhaust calculations for plasma processing."""
 
@@ -51,7 +78,7 @@ def output(self):
             # Double null divertor configuration
             po.ovarre(
                 self.outfile,
-                "Plasma separatrix power over major radius (Pₛₑₚ / R₀) (MW/m) "
+                "Plasma separatrix power over major radius (Pₛₑₚ / R₀) [MW/m] "
                 "(On peak divertor)",
                 "(p_plasma_separatrix_rmajor_mw)",
                 self.data.physics.p_plasma_separatrix_rmajor_mw,
@@ -60,7 +87,7 @@ def output(self):
             po.ovarre(
                 self.outfile,
                 "EU-DEMO divertor protection re-attachment metric (PₛₑₚBₜ / q₉₅AR₀) "
-                "(MWT/m) (On peak divertor)",
+                "[MWT/m] (On peak divertor)",
                 "(p_div_bt_q_aspect_rmajor_mw)",
                 self.data.physics.p_div_bt_q_aspect_rmajor_mw,
                 "OP ",
@@ -69,7 +96,7 @@ def output(self):
             # Single null divertor configuration
             po.ovarre(
                 self.outfile,
-                "Plasma separatrix power over major radius (Pₛₑₚ / R₀) (MW/m)",
+                "Plasma separatrix power over major radius (Pₛₑₚ / R₀) [MW/m]",
                 "(p_plasma_separatrix_rmajor_mw)",
                 self.data.physics.p_plasma_separatrix_rmajor_mw,
                 "OP ",
@@ -77,13 +104,14 @@ def output(self):
             po.ovarre(
                 self.outfile,
                 "EU-DEMO divertor protection re-attachment metric (PₛₑₚBₜ / q₉₅AR₀) "
-                "(MWT/m)",
+                "[MWT/m]",
                 "(p_div_bt_q_aspect_rmajor_mw)",
                 self.data.physics.p_div_bt_q_aspect_rmajor_mw,
                 "OP ",
             )
-
         po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
+        self.output_brunner_divertor_power_splits()
 
     @staticmethod
     def calculate_separatrix_power(
@@ -218,3 +246,200 @@ def calculate_radiation_fraction(
             return 0.0
 
         return p_plasma_rad_mw / p_plasma_heating_mw
+
+    def output_brunner_divertor_power_splits(self):
+        """Output the Brunner divertor power splits to the output file."""
+        if self.data.stellarator.istell == 0:
+            po.osubhd(self.outfile, "Brunner Divertor Power Splits:")
+
+            for op in [
+                (
+                    (
+                        "Requested fraction of power to the lower divertor in double "
+                        "null configuration"
+                    ),
+                    "(f_p_div_lower_separatrix)",
+                    self.data.physics.f_p_div_lower_separatrix,
+                ),
+                (
+                    (
+                        "Required distance between the first and second plasma "
+                        "separatrixes at the outer midplane (δR_sep) [m]"
+                    ),
+                    "(dr_plasma_outboard_midplane_separatrix_separation)",
+                    self.data.physics.dr_plasma_outboard_midplane_separatrix_separation,
+                ),
+                None,
+                (
+                    "Outboard side heat flux decay length (m)",
+                    "(len_sol_outboard_power_decay)",
+                    self.data.physics.len_sol_outboard_power_decay,
+                ),
+                None,
+                (
+                    "Fraction of separatrix power on the inner target(s)",
+                    "(f_p_div_inboard_separatrix)",
+                    self.data.physics.f_p_div_inboard_separatrix,
+                ),
+                (
+                    "Fraction of separatrix power on the outer target(s)",
+                    "(f_p_div_outboard_separatrix)",
+                    self.data.physics.f_p_div_outboard_separatrix,
+                ),
+                None,
+                (
+                    "Fraction of separatrix power on the inner lower target",
+                    "(f_p_div_lower_inboard_separatrix)",
+                    self.data.physics.f_p_div_lower_inboard_separatrix,
+                ),
+                (
+                    "Separatrix power on the inner lower target",
+                    "(p_div_lower_inboard_separatrix_mw)",
+                    self.data.physics.p_div_lower_inboard_separatrix_mw,
+                ),
+                None,
+                (
+                    "Fraction of separatrix power on the outer lower target",
+                    "(f_p_div_lower_outboard_separatrix)",
+                    self.data.physics.f_p_div_lower_outboard_separatrix,
+                ),
+                (
+                    "Separatrix power on the outer lower target",
+                    "(p_div_lower_outboard_separatrix_mw)",
+                    self.data.physics.p_div_lower_outboard_separatrix_mw,
+                ),
+            ]:
+                if op is None:
+                    po.oblnkl(self.outfile)
+                else:
+                    desc, var, val = op
+                    po.ovarre(self.outfile, desc, var, val, "OP ")
+
+            if (
+                DivertorNumberModels(self.data.physics.i_single_null)
+                == DivertorNumberModels.DOUBLE_NULL
+            ):
+                po.oblnkl(self.outfile)
+                po.ocmmnt(self.outfile, "----------------------------")
+                po.oblnkl(self.outfile)
+
+                for op in [
+                    (
+                        "Fraction of separatrix power on the inner upper target",
+                        "(f_p_div_upper_inboard_separatrix)",
+                        self.data.physics.f_p_div_upper_inboard_separatrix,
+                    ),
+                    (
+                        "Separatrix power on the inner upper target",
+                        "(p_div_upper_inboard_separatrix_mw)",
+                        self.data.physics.p_div_upper_inboard_separatrix_mw,
+                    ),
+                    None,
+                    (
+                        "Fraction of separatrix power on the outer upper target",
+                        "(f_p_div_upper_outboard_separatrix)",
+                        self.data.physics.f_p_div_upper_outboard_separatrix,
+                    ),
+                    (
+                        "Separatrix power on the outer upper target",
+                        "(p_div_upper_outboard_separatrix_mw)",
+                        self.data.physics.p_div_upper_outboard_separatrix_mw,
+                    ),
+                ]:
+                    if op is None:
+                        po.oblnkl(self.outfile)
+                    else:
+                        desc, var, val = op
+                        po.ovarre(self.outfile, desc, var, val, "OP ")
+
+
+def calculate_brunner_divertor_power_splits(
+    dr_outboard_midplane_sep: float,
+    len_plasma_sol_outboard_power_decay: float,
+    len_plasma_sol_inboard_power_decay: float,
+) -> DivertorSeparatrixPowerSplits:
+    """
+    Calculate the power splits to the divertor targets using Brunner's method.
+
+    Parameters
+    ----------
+    dr_outboard_midplane_sep : float
+        Radial separation of the plasma outboard midplane separatrix (δR_sep) [m].
+    len_plasma_sol_outboard_power_decay : float
+        Power decay length in the scrape-off layer (λ_q) [m].
+    len_plasma_sol_inboard_power_decay : float
+        Power decay length in the scrape-off layer for the inboard side (λᵢₙ_q) [m].
+
+    Returns
+    -------
+    DivertorSeparatrixPowerSplits
+        Dataclass containing the fractions of total separatrix power to each
+        divertor target.
+
+    Notes
+    -----
+    - The fitted value for `f_p_inner_sep_0` and `f_p_inner_sep_infinity` are taken
+    from the fit given by Petrie et al. [2]
+
+    References
+    ----------
+    [1] D. Brunner, A. Q. Kuang, B. LaBombard, and J. L. Terry, “The dependence of
+    divertor power sharing on magnetic flux balance in near double-null
+    configurations on Alcator C-Mod,” Nuclear Fusion, vol. 58, no. 7, p. 076010,
+    May 2018, doi: https://doi.org/10.1088/1741-4326/aac006.
+
+    [2] T. W. Petrie et al., “The effect of divertor magnetic balance on H-mode
+    performance in DIII-D,” Journal of Nuclear Materials, vol. 290-293, pp. 935-939,
+    Mar. 2001, doi: https://doi.org/10.1016/S0022-3115(00)00492-X
+    """
+    # Fraction of the power to the inner divertors at δR_sep = 0 and δR_sep → ∞
+    f_p_inner_sep_0 = 0.16e0
+    f_p_inner_sep_infinity = 0.41e0
+
+    # Fractions of total outboard power going to each target
+    # Outboard lower divertor
+    f_p_outboard_lower = 1 / (
+        1 + np.exp(dr_outboard_midplane_sep / len_plasma_sol_outboard_power_decay)
+    )
+
+    # Outboard upper divertor
+    f_p_outboard_upper = 1 / (
+        1 + np.exp(-dr_outboard_midplane_sep / len_plasma_sol_outboard_power_decay)
+    )
+
+    # Fractions of the total inboard power going to each target
+    f_p_inboard_lower = 1 / (
+        1 + np.exp(dr_outboard_midplane_sep / (len_plasma_sol_inboard_power_decay))
+    )
+
+    f_p_total_inboard = f_p_inner_sep_0 + (f_p_inner_sep_0 - f_p_inner_sep_infinity) * (
+        1.0e0
+        - (
+            2.0e0
+            / (
+                1.0e0
+                + np.exp(
+                    -(
+                        (dr_outboard_midplane_sep / len_plasma_sol_outboard_power_decay)
+                        ** 2
+                    )
+                )
+            )
+        )
+    )
+
+    f_p_total_outboard = 1.0e0 - f_p_total_inboard
+
+    f_p_inboard_upper = f_p_total_inboard * (1.0e0 - f_p_inboard_lower)
+    f_p_outboard_upper = f_p_total_outboard * (1.0e0 - f_p_outboard_lower)
+    f_p_inboard_lower = f_p_total_inboard * f_p_inboard_lower
+    f_p_outboard_lower = f_p_total_outboard * f_p_outboard_lower
+
+    return DivertorSeparatrixPowerSplits(
+        f_p_div_inboard_separatrix=f_p_total_inboard,
+        f_p_div_outboard_separatrix=f_p_total_outboard,
+        f_p_div_inboard_lower_separatrix=f_p_inboard_lower,
+        f_p_div_inboard_upper_separatrix=f_p_inboard_upper,
+        f_p_div_outboard_lower_separatrix=f_p_outboard_lower,
+        f_p_div_outboard_upper_separatrix=f_p_outboard_upper,
+    )
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 95b6d126a8..a198f75358 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -24,6 +24,7 @@
 )
 from process.models.physics import impurity_radiation
 from process.models.physics.bootstrap_current import BootstrapCurrentFractionModel
+from process.models.physics.exhaust import calculate_brunner_divertor_power_splits
 from process.models.physics.profiles import (
     DensityProfilePedestalType,
     PlasmaProfileShapeType,
@@ -845,15 +846,16 @@ def run(self):
         )
 
         # if double null configuration share the power
-        # over the upper and lower divertor, where self.data.physics.f_p_div_lower gives
-        # the factor of power conducted to the lower divertor
+        # over the upper and lower divertor, where
+        # self.data.physics.f_p_div_lower_separatrix gives the factor of power conducted
+        # to the lower divertor
         if self.data.divertor.n_divertors == 2:
             self.data.physics.p_div_lower_separatrix_mw = (
-                self.data.physics.f_p_div_lower
+                self.data.physics.f_p_div_lower_separatrix
                 * self.data.physics.p_plasma_separatrix_mw
             )
             self.data.physics.p_div_upper_separatrix_mw = (
-                1.0e0 - self.data.physics.f_p_div_lower
+                1.0e0 - self.data.physics.f_p_div_lower_separatrix
             ) * self.data.physics.p_plasma_separatrix_mw
             self.data.physics.p_div_separatrix_max_mw = max(
                 self.data.physics.p_div_lower_separatrix_mw,
@@ -993,79 +995,67 @@ def run(self):
 
         # ============================================================
 
-        # Calculate the target imbalances
-        # find the total power into the targets
-        self.data.physics.ptarmw = self.data.physics.p_plasma_separatrix_mw * (
-            1.0e0 - self.data.physics.rad_fraction_sol
-        )
-        # use self.data.physics.f_p_div_lower to find deltarsep
-        # Parameters taken from double null machine
-        # D. Brunner et al
-
-        # Issue #1559 Infinities in physics_module.drsep when running single null in a
-        # double null machine
+        # Issue #1559 Infinities in
+        # physics_module.dr_plasma_outboard_midplane_separatrix_separation when running
+        # single null in a double null machine
         # C W Ashe
-        if self.data.physics.f_p_div_lower < 4.5e-5:
-            self.data.physics.drsep = 1.5e-2
-        elif self.data.physics.f_p_div_lower > (1.0e0 - 4.5e-5):
-            self.data.physics.drsep = -1.5e-2
+        if self.data.physics.f_p_div_lower_separatrix < 4.5e-5:
+            self.data.physics.dr_plasma_outboard_midplane_separatrix_separation = 1.5e-2
+        elif self.data.physics.f_p_div_lower_separatrix > (1.0e0 - 4.5e-5):
+            self.data.physics.dr_plasma_outboard_midplane_separatrix_separation = -1.5e-2
         else:
-            self.data.physics.drsep = (
+            self.data.physics.dr_plasma_outboard_midplane_separatrix_separation = (
                 -2.0e0
                 * 1.5e-3
-                * math.atanh(2.0e0 * (self.data.physics.f_p_div_lower - 0.5e0))
-            )
-        # Model Taken from D3-D paper for conventional divertor
-        # Journal of Nuclear Materials
-        # Volumes 290-293, March 2001, Pages 935-939
-        # Find the innner and outer lower target imbalance
-
-        self.data.physics.fio = 0.16e0 + (0.16e0 - 0.41e0) * (
-            1.0e0
-            - (
-                2.0e0
-                / (
-                    1.0e0
-                    + np.exp(
-                        -(
-                            (
-                                self.data.physics.drsep
-                                / self.data.physics.len_sol_outboard_power_decay
-                            )
-                            ** 2
-                        )
-                    )
+                * math.atanh(
+                    2.0e0 * (self.data.physics.f_p_div_lower_separatrix - 0.5e0)
                 )
             )
+
+        div_power_plits = calculate_brunner_divertor_power_splits(
+            dr_outboard_midplane_sep=self.data.physics.dr_plasma_outboard_midplane_separatrix_separation,
+            len_plasma_sol_outboard_power_decay=self.data.physics.len_sol_outboard_power_decay,
+            len_plasma_sol_inboard_power_decay=self.data.physics.len_sol_inboard_power_decay,
+        )
+
+        self.data.physics.f_p_div_inboard_separatrix = (
+            div_power_plits.f_p_div_inboard_separatrix
+        )
+        self.data.physics.f_p_div_outboard_separatrix = (
+            div_power_plits.f_p_div_outboard_separatrix
+        )
+        self.data.physics.f_p_div_lower_inboard_separatrix = (
+            div_power_plits.f_p_div_inboard_lower_separatrix
+        )
+        self.data.physics.f_p_div_lower_outboard_separatrix = (
+            div_power_plits.f_p_div_outboard_lower_separatrix
+        )
+        self.data.physics.f_p_div_upper_inboard_separatrix = (
+            div_power_plits.f_p_div_inboard_upper_separatrix
+        )
+        self.data.physics.f_p_div_upper_outboard_separatrix = (
+            div_power_plits.f_p_div_outboard_upper_separatrix
+        )
+
+        self.data.physics.p_div_lower_inboard_separatrix_mw = (
+            self.data.physics.f_p_div_lower_inboard_separatrix
+            * self.data.physics.p_plasma_separatrix_mw
+        )
+
+        self.data.physics.p_div_lower_outboard_separatrix_mw = (
+            self.data.physics.f_p_div_lower_outboard_separatrix
+            * self.data.physics.p_plasma_separatrix_mw
+        )
+
+        self.data.physics.p_div_upper_inboard_separatrix_mw = (
+            self.data.physics.f_p_div_upper_inboard_separatrix
+            * self.data.physics.p_plasma_separatrix_mw
+        )
+
+        self.data.physics.p_div_upper_outboard_separatrix_mw = (
+            self.data.physics.f_p_div_upper_outboard_separatrix
+            * self.data.physics.p_plasma_separatrix_mw
         )
-        if self.data.divertor.n_divertors == 2:
-            # Double Null configuration
-            # Find all the power fractions accross the targets
-            # Taken from D3-D conventional divertor design
-            self.data.physics.fli = (
-                self.data.physics.f_p_div_lower * self.data.physics.fio
-            )
-            self.data.physics.flo = self.data.physics.f_p_div_lower * (
-                1.0e0 - self.data.physics.fio
-            )
-            self.data.physics.fui = (
-                1.0e0 - self.data.physics.f_p_div_lower
-            ) * self.data.physics.fio
-            self.data.physics.fuo = (1.0e0 - self.data.physics.f_p_div_lower) * (
-                1.0e0 - self.data.physics.fio
-            )
-            # power into each target
-            self.data.physics.plimw = self.data.physics.fli * self.data.physics.ptarmw
-            self.data.physics.plomw = self.data.physics.flo * self.data.physics.ptarmw
-            self.data.physics.puimw = self.data.physics.fui * self.data.physics.ptarmw
-            self.data.physics.puomw = self.data.physics.fuo * self.data.physics.ptarmw
-        else:
-            # Single null configuration
-            self.data.physics.fli = self.data.physics.fio
-            self.data.physics.flo = 1.0e0 - self.data.physics.fio
-            # power into each target
-            self.data.physics.plimw = self.data.physics.fli * self.data.physics.ptarmw
-            self.data.physics.plomw = self.data.physics.flo * self.data.physics.ptarmw
 
         # Calculate some derived quantities that may not have been defined earlier
         self.data.physics.p_plasma_heating_total_mw = (
@@ -2193,112 +2183,6 @@ def outplas(self):
             "OP ",
         )
 
-        if self.data.stellarator.istell == 0:
-            po.oblnkl(self.outfile)
-            po.ovarre(
-                self.outfile,
-                "Power incident on the divertor targets (MW)",
-                "(ptarmw)",
-                self.data.physics.ptarmw,
-                "OP ",
-            )
-            po.ovarre(
-                self.outfile,
-                "Divertor poloidal angle subtended by plasma (degrees)",
-                "(deg_div_poloidal_plasma)",
-                self.data.divertor.deg_div_poloidal_plasma,
-                "OP ",
-            )
-            po.ovarre(
-                self.outfile,
-                "Fraction of power to the lower divertor",
-                "(f_p_div_lower)",
-                self.data.physics.f_p_div_lower,
-                "IP ",
-            )
-            po.ovarre(
-                self.outfile,
-                "Outboard side heat flux decay length (m)",
-                "(len_sol_outboard_power_decay)",
-                self.data.physics.len_sol_outboard_power_decay,
-                "OP ",
-            )
-            if self.data.divertor.n_divertors == 2:
-                po.ovarre(
-                    self.outfile,
-                    "Midplane separation of the two magnetic closed flux surfaces (m)",
-                    "(drsep)",
-                    self.data.physics.drsep,
-                    "OP ",
-                )
-
-            po.ovarre(
-                self.outfile,
-                "Fraction of power on the inner targets",
-                "(fio)",
-                self.data.physics.fio,
-                "OP ",
-            )
-            po.ovarre(
-                self.outfile,
-                "Fraction of power incident on the lower inner target",
-                "(fLI)",
-                self.data.physics.fli,
-                "OP ",
-            )
-            po.ovarre(
-                self.outfile,
-                "Fraction of power incident on the lower outer target",
-                "(fLO)",
-                self.data.physics.flo,
-                "OP ",
-            )
-            if self.data.divertor.n_divertors == 2:
-                po.ovarre(
-                    self.outfile,
-                    "Fraction of power incident on the upper inner target",
-                    "(fUI)",
-                    self.data.physics.fui,
-                    "OP ",
-                )
-                po.ovarre(
-                    self.outfile,
-                    "Fraction of power incident on the upper outer target",
-                    "(fUO)",
-                    self.data.physics.fuo,
-                    "OP ",
-                )
-
-            po.ovarre(
-                self.outfile,
-                "Power incident on the lower inner target (MW)",
-                "(pLImw)",
-                self.data.physics.plimw,
-                "OP ",
-            )
-            po.ovarre(
-                self.outfile,
-                "Power incident on the lower outer target (MW)",
-                "(pLOmw)",
-                self.data.physics.plomw,
-                "OP ",
-            )
-            if self.data.divertor.n_divertors == 2:
-                po.ovarre(
-                    self.outfile,
-                    "Power incident on the upper innner target (MW)",
-                    "(pUImw)",
-                    self.data.physics.puimw,
-                    "OP ",
-                )
-                po.ovarre(
-                    self.outfile,
-                    "Power incident on the upper outer target (MW)",
-                    "(pUOmw)",
-                    self.data.physics.puomw,
-                    "OP ",
-                )
-
         po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
diff --git a/process/models/physics/scrape_off_layer.py b/process/models/physics/scrape_off_layer.py
index 37fc3d06a6..4388d79914 100644
--- a/process/models/physics/scrape_off_layer.py
+++ b/process/models/physics/scrape_off_layer.py
@@ -71,6 +71,11 @@ def run(self):
                 self.data.physics.len_plasma_sol_mast14_power_decay_2
             )
 
+        self.data.physics.len_sol_inboard_power_decay = (
+            self.data.physics.f_len_sol_power_decay_inboard_outboard
+            * self.data.physics.len_sol_outboard_power_decay
+        )
+
         self.data.physics.a_plasma_outboard_sol_parallel = self.calculate_upstream_sol_outboard_parallel_area(  # noqa: E501
             rmajor=self.data.physics.rmajor,
             rminor=self.data.physics.rminor,
@@ -105,7 +110,7 @@ def output(self) -> None:
 
         po.ovarre(
             self.outfile,
-            "Outboard SOL power decay length (λ_q) [m]",
+            "Outboard SOL power decay length (λₒᵤₜ_q) [m]",
             "(len_sol_outboard_power_decay)",
             self.data.physics.len_sol_outboard_power_decay,
         )
@@ -118,6 +123,19 @@ def output(self) -> None:
             + " ",
         )
         po.oblnkl(self.outfile)
+        po.ovarre(
+            self.outfile,
+            "Inboard to outboard SOL power decay length ratio (λᵢₙ_q/λₒᵤₜ_q)",
+            "(f_len_sol_power_decay_inboard_outboard)",
+            self.data.physics.f_len_sol_power_decay_inboard_outboard,
+        )
+        po.ovarre(
+            self.outfile,
+            "Inboard SOL power decay length (λᵢₙ_q) [m]",
+            "(len_sol_inboard_power_decay)",
+            self.data.physics.len_sol_inboard_power_decay,
+        )
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "Eich 2013 SOL power decay length (λ_q) [m]",
diff --git a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
index 93f3541265..b3e92c2df9 100644
--- a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
+++ b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT
@@ -263,7 +263,7 @@ b_plasma_toroidal_on_axis       = 3.0 * toroidal field on axis (T) (`iteration v
 nd_plasma_electrons_vol_avg     = 9.69888313737236e+19 * electron density (/m3) (`iteration variable 6`)
 beta_norm_max = 5.0 * Troyon-like coefficient for beta scaling
 i_beta_norm_max = 0 * Normalised beta max selection (0 = user input)
-f_p_div_lower     = 0.5 * fraction of power to the lower divertor in double null configuration
+f_p_div_lower_separatrix     = 0.5 * fraction of power to the lower divertor in double null configuration
 f_nd_plasma_pedestal_greenwald   = 0.1 * fraction of Greenwald density to set as pedestal-top density; If `<0`; pedestal-top
 f_nd_plasma_separatrix_greenwald   = 0.1 * fraction of Greenwald density to set as separatrix density; If `<0`; separatrix
 hfact    = 1.2 * H factor on energy confinement times; radiation corrected (`iteration variable 10`);
diff --git a/tests/regression/input_files/st_regression.IN.DAT b/tests/regression/input_files/st_regression.IN.DAT
index a34de1f2e9..2a346db5ac 100644
--- a/tests/regression/input_files/st_regression.IN.DAT
+++ b/tests/regression/input_files/st_regression.IN.DAT
@@ -631,7 +631,7 @@ f_nd_impurity_electrons(14) = 5e-05    *Tungsten
 *-------------------------------- *Divertor* ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------*
 *‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾*
 
-f_p_div_lower = 0.5 
+f_p_div_lower_separatrix = 0.5 
 * DESCRIPTION:   Fraction of power to the lower divertor
 * JUSTIFICATION: All power to lower divertor as Single-Null
 
diff --git a/tests/unit/models/physics/test_plasma_exhaust.py b/tests/unit/models/physics/test_plasma_exhaust.py
new file mode 100644
index 0000000000..1bce228473
--- /dev/null
+++ b/tests/unit/models/physics/test_plasma_exhaust.py
@@ -0,0 +1,46 @@
+import pytest
+
+from process.models.physics.exhaust import (
+    calculate_brunner_divertor_power_splits,
+)
+
+
+def test_calculate_brunner_divertor_power_splits_at_zero_separatrix_separation():
+    result = calculate_brunner_divertor_power_splits(
+        dr_outboard_midplane_sep=0.0,
+        len_plasma_sol_outboard_power_decay=0.01,
+        len_plasma_sol_inboard_power_decay=0.01,
+    )
+
+    assert result.f_p_div_inboard_separatrix == pytest.approx(0.16)
+    assert result.f_p_div_outboard_separatrix == pytest.approx(0.84)
+    assert result.f_p_div_inboard_lower_separatrix == pytest.approx(0.08)
+    assert result.f_p_div_inboard_upper_separatrix == pytest.approx(0.08)
+    assert result.f_p_div_outboard_lower_separatrix == pytest.approx(0.42)
+    assert result.f_p_div_outboard_upper_separatrix == pytest.approx(0.42)
+
+
+def test_calculate_brunner_divertor_power_splits_asymmetric_case():
+    result = calculate_brunner_divertor_power_splits(
+        dr_outboard_midplane_sep=0.01,
+        len_plasma_sol_outboard_power_decay=0.005,
+        len_plasma_sol_inboard_power_decay=0.0075,
+    )
+
+    assert result.f_p_div_inboard_separatrix == pytest.approx(0.4010068950189542)
+    assert result.f_p_div_outboard_separatrix == pytest.approx(0.5989931049810458)
+    assert result.f_p_div_inboard_lower_separatrix == pytest.approx(0.08365345781749392)
+    assert result.f_p_div_inboard_upper_separatrix == pytest.approx(0.31735343720146025)
+    assert result.f_p_div_outboard_lower_separatrix == pytest.approx(0.07140172838484167)
+    assert result.f_p_div_outboard_upper_separatrix == pytest.approx(0.527591376596204)
+
+    assert (
+        result.f_p_div_inboard_lower_separatrix + result.f_p_div_inboard_upper_separatrix
+    ) == pytest.approx(result.f_p_div_inboard_separatrix)
+    assert (
+        result.f_p_div_outboard_lower_separatrix
+        + result.f_p_div_outboard_upper_separatrix
+    ) == pytest.approx(result.f_p_div_outboard_separatrix)
+    assert (
+        result.f_p_div_inboard_separatrix + result.f_p_div_outboard_separatrix
+    ) == pytest.approx(1.0)
diff --git a/tests/unit/models/test_divertor.py b/tests/unit/models/test_divertor.py
index 803df7cfc7..1ca1bf2ac3 100644
--- a/tests/unit/models/test_divertor.py
+++ b/tests/unit/models/test_divertor.py
@@ -86,7 +86,7 @@ def test_divwade(monkeypatch, divertor):
         nd_plasma_separatrix_electron = 1.0e19
         deg_div_field_plate = 5.0
         rad_fraction_sol = 8.0e-1
-        f_p_div_lower = 1.0
+        f_p_div_lower_separatrix = 1.0
 
         expected_pflux_div_heat_load_mw = 0.58898578
 
@@ -101,7 +101,7 @@ def test_divwade(monkeypatch, divertor):
             nd_plasma_separatrix_electron,
             deg_div_field_plate,
             rad_fraction_sol,
-            f_p_div_lower,
+            f_p_div_lower_separatrix,
             False,
         )