From cce07535f448be581bfe48365485bbf3d7fff198 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Sun, 16 Aug 2026 19:40:41 +0530
Subject: [PATCH 1/3] Add Zhong 2008 geoid postprocessing

Add pure Appendix A no-self-gravity and self-gravity response operators with explicit internal-load scaling, density-contrast naming, dimensional constant documentation, and harmonic-degree validation.

Provide a rotated-Stokes adapter that delegates normal-traction recovery to the existing Stokes.boundary_normal_traction API, projects only the requested P_l^0 response, and avoids duplicate CBF, constrained, or dynamic-topography implementations.

Keep the feature in uw.postprocessing rather than adding a Zhong-specific facade to the generic Stokes solver. Add focused formula, serial end-to-end, and two/four-rank MPI validation.
---
 docs/developer/CHANGELOG.md                   |   4 +
 ...ry-stress-and-projection-postprocessing.md |  47 +++
 src/underworld3/__init__.py                   |   1 +
 src/underworld3/postprocessing/__init__.py    |  26 ++
 src/underworld3/postprocessing/geoid.py       | 387 ++++++++++++++++++
 .../test_1071_zhong2008_geoid_parallel.py     |  76 ++++
 tests/test_1070_postprocessing_geoid.py       | 171 ++++++++
 7 files changed, 712 insertions(+)
 create mode 100644 src/underworld3/postprocessing/__init__.py
 create mode 100644 src/underworld3/postprocessing/geoid.py
 create mode 100644 tests/parallel/test_1071_zhong2008_geoid_parallel.py
 create mode 100644 tests/test_1070_postprocessing_geoid.py

diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md
index 0b6f16071..799e9b2a3 100644
--- a/docs/developer/CHANGELOG.md
+++ b/docs/developer/CHANGELOG.md
@@ -273,6 +273,10 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
   seam sensitivity is unchanged (#404).
 - Recorded as the preferred free-slip BC in the project guidance (#300);
   conda PETSc floor raised to ≥ 3.25 for FMG/rotation API consistency (#304).
+- `uw.postprocessing.zhong2008_response_from_rotated_stokes(...)` projects the
+  existing rotated-free-slip boundary traction onto the benchmark harmonic and
+  applies one Appendix A operator for geoid and self-gravity response. The pure
+  operator functions also accept precomputed topography coefficients.
 
 ### Generalized Geometric Multigrid via Custom Prolongation (July 2026)
 
diff --git a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
index e23e2e2f1..90126e2b2 100644
--- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
+++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
@@ -74,6 +74,53 @@ Combine the two rules — project the τ components with `linear_solver()` for t
 cheap linear solve, then compose `σ_rr` analytically — for accurate boundary
 stress recovery that also scales.
 
+## 3. Zhong spherical-shell geoid response
+
+Rotated free slip already exposes normal traction through
+`Stokes.boundary_normal_traction()`. The Zhong adapter projects that recovered
+traction onto the unnormalised `P_l^0` benchmark harmonic and applies the
+Appendix A geoid operator:
+
+```python
+stokes.solve()
+
+response = uw.postprocessing.zhong2008_response_from_rotated_stokes(
+    stokes=stokes,
+    radius_inner=0.55,
+    radius_outer=1.0,
+    radius_internal=0.775,
+    harmonic_degree=2,
+    load_scale=1.0,
+    include_self_gravity=True,
+)
+```
+
+The adapter delegates stress recovery to the existing rotated-free-slip API;
+it does not implement a second CBF, constrained-multiplier, or topography
+recovery path. `load_scale` must match the coefficient multiplying the model's
+internal `P_l^0` load.
+
+When surface and CMB topography coefficients are already available, call
+`zhong2008_geoid_response()` or `zhong2008_self_gravity_response()` directly.
+These functions are pure post-processing and do not require a Stokes object.
+
+The Zhong Appendix A calculation is expressed as one linear operator:
+
+```text
+N = G h + n_load
+(I - Q G) h_self_gravity = h + Q n_load
+```
+
+where `h` contains surface and CMB topography, `N` contains their geoid
+responses, and `Q` contains the two self-gravity density factors. Sharing `G`
+and `n_load` between the no-self-gravity and self-gravity paths avoids separate
+scalar implementations of the same coefficients.
+
+The rotated harmonic projector gathers boundary samples to rank zero and
+reconstructs their spherical triangulation. A future boundary-reaction
+functional could replace this step with a direct distributed finite-element
+projection without changing the Zhong operator API.
+
 ## See also
 
 - Issues [#156] (projection solver settings), [#157] (projection memory),
diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py
index 818ff2339..450e23680 100644
--- a/src/underworld3/__init__.py
+++ b/src/underworld3/__init__.py
@@ -201,6 +201,7 @@ def view():
 import underworld3.maths
 import underworld3.swarm
 import underworld3.systems
+import underworld3.postprocessing
 import underworld3.maths
 import underworld3.utilities
 import underworld3.model
diff --git a/src/underworld3/postprocessing/__init__.py b/src/underworld3/postprocessing/__init__.py
new file mode 100644
index 000000000..57997e486
--- /dev/null
+++ b/src/underworld3/postprocessing/__init__.py
@@ -0,0 +1,26 @@
+r"""
+Post-processing helpers for Underworld3 models.
+
+The post-processing package contains derived diagnostics that are useful across
+multiple model scripts but do not belong in a solver implementation.
+"""
+
+from . import geoid
+from .geoid import (
+    Zhong2008DynamicResponse,
+    Zhong2008GeoidResponse,
+    Zhong2008SelfGravityResponse,
+    zhong2008_geoid_response,
+    zhong2008_response_from_rotated_stokes,
+    zhong2008_self_gravity_response,
+)
+
+__all__ = [
+    "geoid",
+    "Zhong2008DynamicResponse",
+    "Zhong2008GeoidResponse",
+    "Zhong2008SelfGravityResponse",
+    "zhong2008_geoid_response",
+    "zhong2008_response_from_rotated_stokes",
+    "zhong2008_self_gravity_response",
+]
diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py
new file mode 100644
index 000000000..a59cb328a
--- /dev/null
+++ b/src/underworld3/postprocessing/geoid.py
@@ -0,0 +1,387 @@
+r"""Zhong et al. (2008) spherical-shell geoid response functions."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from numbers import Integral
+
+from mpi4py import MPI
+import numpy as np
+
+
+@dataclass(frozen=True)
+class Zhong2008GeoidResponse:
+    """No-self-gravity surface and CMB geoid coefficients."""
+
+    surface_geoid: float
+    cmb_geoid: float
+
+
+@dataclass(frozen=True)
+class Zhong2008SelfGravityResponse:
+    """Self-gravity-corrected topography and geoid coefficients."""
+
+    surface_topography: float
+    cmb_topography: float
+    surface_geoid: float
+    cmb_geoid: float
+    q_surface: float
+    q_cmb: float
+
+
+@dataclass(frozen=True)
+class Zhong2008DynamicResponse:
+    """Rotated-Stokes topography and corresponding geoid coefficients."""
+
+    surface_topography: float
+    cmb_topography: float
+    surface_geoid: float
+    cmb_geoid: float
+    self_gravity: Zhong2008SelfGravityResponse | None = None
+
+
+def _validate_geometry(
+    radius_inner: float,
+    radius_outer: float,
+    radius_internal: float,
+    harmonic_degree: int,
+) -> tuple[float, float, float, int]:
+    try:
+        ri = float(radius_inner)
+        ro = float(radius_outer)
+        rint = float(radius_internal)
+    except (TypeError, ValueError) as error:
+        raise TypeError("The shell radii must be real numbers.") from error
+    if not np.all(np.isfinite((ri, rint, ro))) or not 0.0 < ri < rint < ro:
+        raise ValueError(
+            "Expected finite radii ordered as "
+            "0 < radius_inner < radius_internal < radius_outer."
+        )
+    if isinstance(harmonic_degree, bool) or not isinstance(harmonic_degree, Integral):
+        raise TypeError("harmonic_degree must be an integer.")
+    degree = int(harmonic_degree)
+    if degree < 1:
+        raise ValueError("harmonic_degree must be at least one.")
+    return ri, ro, rint, degree
+
+
+def _zhong2008_geoid_operator(
+    radius_inner: float,
+    radius_outer: float,
+    radius_internal: float,
+    harmonic_degree: int,
+    load_scale: float,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Return the Zhong Appendix A topography operator and load vector."""
+
+    ri, ro, rint, degree = _validate_geometry(
+        radius_inner,
+        radius_outer,
+        radius_internal,
+        harmonic_degree,
+    )
+    load_scale = float(load_scale)
+    if not np.isfinite(load_scale):
+        raise ValueError("load_scale must be finite.")
+
+    denominator = float(2 * degree + 1)
+    operator = np.array(
+        [
+            [
+                ro / denominator,
+                ri * (ri / ro) ** (degree + 1) / denominator,
+            ],
+            [
+                ro * (ri / ro) ** degree / denominator,
+                ri / denominator,
+            ],
+        ],
+        dtype=float,
+    )
+    load = load_scale * np.array(
+        [
+            -rint * (rint / ro) ** (degree + 1) / denominator,
+            -rint * (ri / rint) ** degree / denominator,
+        ],
+        dtype=float,
+    )
+    return operator, load
+
+
+def zhong2008_geoid_response(
+    *,
+    radius_inner: float,
+    radius_outer: float,
+    radius_internal: float,
+    harmonic_degree: int,
+    surface_topography_coefficient: float,
+    cmb_topography_coefficient: float,
+    load_scale: float = 1.0,
+) -> Zhong2008GeoidResponse:
+    r"""Return the nondimensional Zhong no-self-gravity geoid coefficients.
+
+    The topography coefficients and ``load_scale`` must use the same
+    unnormalised :math:`P_l^0` forcing convention as the Stokes solve.
+    """
+
+    operator, load = _zhong2008_geoid_operator(
+        radius_inner,
+        radius_outer,
+        radius_internal,
+        harmonic_degree,
+        load_scale,
+    )
+    topography = np.array(
+        [surface_topography_coefficient, cmb_topography_coefficient],
+        dtype=float,
+    )
+    if not np.all(np.isfinite(topography)):
+        raise ValueError("The topography coefficients must be finite.")
+    geoid = operator @ topography + load
+    return Zhong2008GeoidResponse(float(geoid[0]), float(geoid[1]))
+
+
+def zhong2008_self_gravity_response(
+    *,
+    radius_inner: float,
+    radius_outer: float,
+    radius_internal: float,
+    harmonic_degree: int,
+    surface_topography_coefficient: float,
+    cmb_topography_coefficient: float,
+    load_scale: float = 1.0,
+    surface_density_contrast: float = 3300.0,
+    cmb_density_contrast: float = 5400.0,
+    planet_radius: float = 6370000.0,
+    gravity: float = 9.8,
+    gravitational_constant: float = 6.67e-11,
+) -> Zhong2008SelfGravityResponse:
+    r"""Return the Zhong self-gravity-corrected response coefficients.
+
+    Density contrasts are in kg/m³, ``planet_radius`` in metres, ``gravity`` in
+    m/s², and ``gravitational_constant`` in SI units. Shell radii and response
+    coefficients remain nondimensional.
+    """
+
+    operator, load = _zhong2008_geoid_operator(
+        radius_inner,
+        radius_outer,
+        radius_internal,
+        harmonic_degree,
+        load_scale,
+    )
+    topography = np.array(
+        [surface_topography_coefficient, cmb_topography_coefficient],
+        dtype=float,
+    )
+    physical_values = np.array(
+        [
+            surface_density_contrast,
+            cmb_density_contrast,
+            planet_radius,
+            gravity,
+            gravitational_constant,
+        ],
+        dtype=float,
+    )
+    if not np.all(np.isfinite(topography)):
+        raise ValueError("The topography coefficients must be finite.")
+    if not np.all(np.isfinite(physical_values)) or np.any(physical_values <= 0.0):
+        raise ValueError("Density contrasts and physical constants must be positive.")
+
+    gravity_scale = 4.0 * np.pi * gravitational_constant * planet_radius / gravity
+    q = gravity_scale * np.array(
+        [surface_density_contrast, cmb_density_contrast], dtype=float
+    )
+    q_matrix = np.diag(q)
+    corrected_topography = np.linalg.solve(
+        np.eye(2) - q_matrix @ operator,
+        topography + q_matrix @ load,
+    )
+    corrected_geoid = operator @ corrected_topography + load
+
+    return Zhong2008SelfGravityResponse(
+        surface_topography=float(corrected_topography[0]),
+        cmb_topography=float(corrected_topography[1]),
+        surface_geoid=float(corrected_geoid[0]),
+        cmb_geoid=float(corrected_geoid[1]),
+        q_surface=float(q[0]),
+        q_cmb=float(q[1]),
+    )
+
+
+def _spherical_triangle_area(a, b, c, radius: float) -> float:
+    determinant = abs(float(np.dot(a, np.cross(b, c))))
+    denominator = float(1.0 + np.dot(a, b) + np.dot(b, c) + np.dot(c, a))
+    return 2.0 * np.arctan2(determinant, denominator) * radius**2
+
+
+def _project_spherical_harmonic_samples(
+    coords: np.ndarray,
+    values: np.ndarray,
+    radius: float,
+    harmonic_degree: int,
+    response_sign: float,
+) -> float:
+    from scipy.spatial import ConvexHull
+
+    unit_coords = coords / np.linalg.norm(coords, axis=1)[:, None]
+    hull = ConvexHull(unit_coords, qhull_options="QJ")
+    projected_integral = 0.0
+    harmonic = np.polynomial.legendre.Legendre.basis(harmonic_degree)
+
+    for simplex in hull.simplices:
+        area = _spherical_triangle_area(
+            unit_coords[simplex[0]],
+            unit_coords[simplex[1]],
+            unit_coords[simplex[2]],
+            radius,
+        )
+        centroid = unit_coords[simplex].mean(axis=0)
+        centroid /= np.linalg.norm(centroid)
+        projected_integral += (
+            area * float(values[simplex].mean()) * harmonic(np.clip(centroid[2], -1, 1))
+        )
+
+    harmonic_norm = 4.0 * np.pi / (2 * harmonic_degree + 1)
+    return float(response_sign * projected_integral / (radius**2 * harmonic_norm))
+
+
+def _rotated_topography_coefficient(
+    stokes,
+    boundary: str,
+    radius: float,
+    harmonic_degree: int,
+    buoyancy_scale: float,
+    response_sign: float,
+) -> float:
+    buoyancy_scale = float(buoyancy_scale)
+    if not np.isfinite(buoyancy_scale) or buoyancy_scale == 0.0:
+        raise ValueError("Boundary buoyancy scales must be finite and nonzero.")
+
+    coords, sigma_nn = stokes.boundary_normal_traction(boundary, mass="auto")
+    local_rows = np.column_stack(
+        (
+            np.asarray(coords, dtype=float),
+            -np.asarray(sigma_nn, dtype=float) / buoyancy_scale,
+        )
+    )
+    gathered_rows = MPI.COMM_WORLD.gather(local_rows, root=0)
+
+    coefficient = None
+    root_error = None
+    if MPI.COMM_WORLD.rank == 0:
+        try:
+            nonempty_rows = [rows for rows in gathered_rows if rows.size]
+            if not nonempty_rows:
+                raise RuntimeError(f"No samples found on boundary {boundary!r}.")
+            merged = {}
+            for row in np.vstack(nonempty_rows):
+                key = tuple(np.round(row[:3], 12))
+                if key not in merged:
+                    merged[key] = [row[:3].copy(), 0.0, 0]
+                merged[key][1] += float(row[3])
+                merged[key][2] += 1
+            global_coords = np.array([item[0] for item in merged.values()])
+            global_values = np.array([item[1] / item[2] for item in merged.values()])
+            coefficient = _project_spherical_harmonic_samples(
+                global_coords,
+                global_values,
+                radius,
+                harmonic_degree,
+                response_sign,
+            )
+        except Exception as error:
+            root_error = f"{type(error).__name__}: {error}"
+
+    root_error = MPI.COMM_WORLD.bcast(root_error, root=0)
+    if root_error is not None:
+        raise RuntimeError(f"Spherical-harmonic projection failed: {root_error}")
+    return float(MPI.COMM_WORLD.bcast(coefficient, root=0))
+
+
+def zhong2008_response_from_rotated_stokes(
+    *,
+    stokes,
+    radius_inner: float,
+    radius_outer: float,
+    radius_internal: float,
+    harmonic_degree: int,
+    load_scale: float = 1.0,
+    surface_boundary: str = "Upper",
+    cmb_boundary: str = "Lower",
+    surface_buoyancy_scale: float = 1.0,
+    cmb_buoyancy_scale: float = 1.0,
+    include_self_gravity: bool = False,
+    surface_density_contrast: float = 3300.0,
+    cmb_density_contrast: float = 5400.0,
+    planet_radius: float = 6370000.0,
+    gravity: float = 9.8,
+    gravitational_constant: float = 6.67e-11,
+) -> Zhong2008DynamicResponse:
+    r"""Compute Zhong response coefficients from a rotated-free-slip Stokes solve.
+
+    Normal traction recovery is delegated to the existing
+    :meth:`Stokes.boundary_normal_traction` implementation. This adapter only
+    projects the two boundary responses onto the unnormalised
+    :math:`P_l^0` harmonic and applies the Zhong Appendix A operator.
+    """
+
+    ri, ro, _, degree = _validate_geometry(
+        radius_inner,
+        radius_outer,
+        radius_internal,
+        harmonic_degree,
+    )
+    if not isinstance(include_self_gravity, bool):
+        raise TypeError("include_self_gravity must be True or False.")
+
+    surface_topography = _rotated_topography_coefficient(
+        stokes,
+        surface_boundary,
+        ro,
+        degree,
+        surface_buoyancy_scale,
+        1.0,
+    )
+    cmb_topography = _rotated_topography_coefficient(
+        stokes,
+        cmb_boundary,
+        ri,
+        degree,
+        cmb_buoyancy_scale,
+        -1.0,
+    )
+    geoid = zhong2008_geoid_response(
+        radius_inner=ri,
+        radius_outer=ro,
+        radius_internal=radius_internal,
+        harmonic_degree=degree,
+        surface_topography_coefficient=surface_topography,
+        cmb_topography_coefficient=cmb_topography,
+        load_scale=load_scale,
+    )
+    self_gravity = None
+    if include_self_gravity:
+        self_gravity = zhong2008_self_gravity_response(
+            radius_inner=ri,
+            radius_outer=ro,
+            radius_internal=radius_internal,
+            harmonic_degree=degree,
+            surface_topography_coefficient=surface_topography,
+            cmb_topography_coefficient=cmb_topography,
+            load_scale=load_scale,
+            surface_density_contrast=surface_density_contrast,
+            cmb_density_contrast=cmb_density_contrast,
+            planet_radius=planet_radius,
+            gravity=gravity,
+            gravitational_constant=gravitational_constant,
+        )
+
+    return Zhong2008DynamicResponse(
+        surface_topography=surface_topography,
+        cmb_topography=cmb_topography,
+        surface_geoid=geoid.surface_geoid,
+        cmb_geoid=geoid.cmb_geoid,
+        self_gravity=self_gravity,
+    )
diff --git a/tests/parallel/test_1071_zhong2008_geoid_parallel.py b/tests/parallel/test_1071_zhong2008_geoid_parallel.py
new file mode 100644
index 000000000..a5284d97b
--- /dev/null
+++ b/tests/parallel/test_1071_zhong2008_geoid_parallel.py
@@ -0,0 +1,76 @@
+import numpy as np
+import pytest
+import sympy
+import underworld3 as uw
+
+
+pytestmark = pytest.mark.level_2
+
+
+def test_zhong2008_rotated_geoid_matches_serial_reference():
+    if uw.mpi.size == 1:
+        pytest.skip("Run this regression with at least two MPI ranks.")
+
+    radius_inner = 0.55
+    radius_outer = 1.0
+    radius_internal = 0.775
+    mesh = uw.meshing.SphericalShellInternalBoundary(
+        radiusOuter=radius_outer,
+        radiusInternal=radius_internal,
+        radiusInner=radius_inner,
+        cellSize=0.25,
+        qdegree=2,
+        degree=1,
+    )
+    velocity = uw.discretisation.MeshVariable("U_geoid_mpi", mesh, mesh.dim, degree=2)
+    pressure = uw.discretisation.MeshVariable(
+        "P_geoid_mpi", mesh, 1, degree=1, continuous=True
+    )
+    stokes = uw.systems.Stokes(
+        mesh,
+        velocityField=velocity,
+        pressureField=pressure,
+    )
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+
+    theta = mesh.CoordinateSystem.xR[1]
+    unit_r = mesh.CoordinateSystem.unit_e_0
+    load = sympy.assoc_legendre(2, 0, sympy.cos(theta)) * unit_r
+    stokes.add_natural_bc(load, "Internal")
+    stokes.add_rotated_freeslip_bc(0, "Upper", normal=unit_r)
+    stokes.add_rotated_freeslip_bc(0, "Lower", normal=-unit_r)
+    stokes.petsc_use_nullspace = True
+    stokes.petsc_options["snes_type"] = "ksponly"
+    stokes.tolerance = 1.0e-5
+    stokes.solve()
+
+    response = uw.postprocessing.zhong2008_response_from_rotated_stokes(
+        stokes=stokes,
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        radius_internal=radius_internal,
+        harmonic_degree=2,
+        load_scale=1.0,
+        include_self_gravity=True,
+    )
+    values = np.array(
+        [
+            response.surface_topography,
+            response.cmb_topography,
+            response.surface_geoid,
+            response.cmb_geoid,
+            response.self_gravity.surface_topography,
+            response.self_gravity.cmb_topography,
+            response.self_gravity.surface_geoid,
+            response.self_gravity.cmb_geoid,
+        ]
+    )
+
+    for rank_values in uw.mpi.comm.allgather(values):
+        assert np.array_equal(rank_values, values)
+
+    zhong_table_2 = np.array(
+        [0.41920, 0.77060, 0.02579, 0.03206, 0.49980, 0.93130, 0.04486, 0.05461]
+    )
+    assert np.allclose(values, zhong_table_2, rtol=0.10, atol=1.0e-12)
diff --git a/tests/test_1070_postprocessing_geoid.py b/tests/test_1070_postprocessing_geoid.py
new file mode 100644
index 000000000..9f8b794b3
--- /dev/null
+++ b/tests/test_1070_postprocessing_geoid.py
@@ -0,0 +1,171 @@
+import math
+
+import numpy as np
+import pytest
+import sympy
+import underworld3 as uw
+
+
+pytestmark = pytest.mark.level_2
+
+
+def test_zhong2008_postprocessing_is_public():
+    assert hasattr(uw, "postprocessing")
+    assert hasattr(uw.postprocessing, "zhong2008_geoid_response")
+    assert hasattr(uw.postprocessing, "zhong2008_self_gravity_response")
+    assert hasattr(uw.postprocessing, "zhong2008_response_from_rotated_stokes")
+
+
+def test_zhong2008_geoid_response_matches_direct_formula_and_load_scale():
+    radius_inner = 0.55
+    radius_outer = 1.0
+    radius_internal = 0.775
+    degree = 2
+    load_scale = 1.75
+    surface_topography = 0.3918264884592169
+    cmb_topography = 0.2653834107104151
+
+    response = uw.postprocessing.zhong2008_geoid_response(
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        radius_internal=radius_internal,
+        harmonic_degree=degree,
+        surface_topography_coefficient=surface_topography,
+        cmb_topography_coefficient=cmb_topography,
+        load_scale=load_scale,
+    )
+
+    denominator = 2 * degree + 1
+    surface_expected = (
+        radius_inner * (radius_inner / radius_outer) ** (degree + 1) * cmb_topography
+        + radius_outer * surface_topography
+        - load_scale
+        * radius_internal
+        * (radius_internal / radius_outer) ** (degree + 1)
+    ) / denominator
+    cmb_expected = (
+        radius_inner * cmb_topography
+        + radius_outer * (radius_inner / radius_outer) ** degree * surface_topography
+        - load_scale * radius_internal * (radius_inner / radius_internal) ** degree
+    ) / denominator
+
+    assert math.isclose(response.surface_geoid, surface_expected)
+    assert math.isclose(response.cmb_geoid, cmb_expected)
+
+
+def test_zhong2008_self_gravity_response_matches_direct_solve():
+    surface_topography = 0.3918264884592169
+    cmb_topography = 0.2653834107104151
+    response = uw.postprocessing.zhong2008_self_gravity_response(
+        radius_inner=0.55,
+        radius_outer=1.0,
+        radius_internal=0.775,
+        harmonic_degree=2,
+        surface_topography_coefficient=surface_topography,
+        cmb_topography_coefficient=cmb_topography,
+        load_scale=1.0,
+        surface_density_contrast=3300.0,
+        cmb_density_contrast=5400.0,
+    )
+
+    q_surface = 4.0 * np.pi * 6.67e-11 * 6370000.0 * 3300.0 / 9.8
+    q_cmb = 4.0 * np.pi * 6.67e-11 * 6370000.0 * 5400.0 / 9.8
+    denominator = 5.0
+    operator = np.array(
+        [
+            [1.0 / denominator, 0.55 * 0.55**3 / denominator],
+            [0.55**2 / denominator, 0.55 / denominator],
+        ]
+    )
+    load = np.array(
+        [
+            -0.775 * 0.775**3 / denominator,
+            -0.775 * (0.55 / 0.775) ** 2 / denominator,
+        ]
+    )
+    q_matrix = np.diag([q_surface, q_cmb])
+    expected_topography = np.linalg.solve(
+        np.eye(2) - q_matrix @ operator,
+        np.array([surface_topography, cmb_topography]) + q_matrix @ load,
+    )
+    expected_geoid = operator @ expected_topography + load
+
+    assert math.isclose(response.q_surface, q_surface)
+    assert math.isclose(response.q_cmb, q_cmb)
+    assert np.allclose(
+        [response.surface_topography, response.cmb_topography],
+        expected_topography,
+    )
+    assert np.allclose(
+        [response.surface_geoid, response.cmb_geoid],
+        expected_geoid,
+    )
+
+
+@pytest.mark.parametrize("harmonic_degree", [0, -1, 2.0, True])
+def test_zhong2008_geoid_rejects_invalid_harmonic_degree(harmonic_degree):
+    error = TypeError if isinstance(harmonic_degree, (float, bool)) else ValueError
+    with pytest.raises(error):
+        uw.postprocessing.zhong2008_geoid_response(
+            radius_inner=0.55,
+            radius_outer=1.0,
+            radius_internal=0.775,
+            harmonic_degree=harmonic_degree,
+            surface_topography_coefficient=0.4,
+            cmb_topography_coefficient=0.7,
+        )
+
+
+def test_zhong2008_rotated_stokes_adapter_matches_table_2():
+    radius_inner = 0.55
+    radius_outer = 1.0
+    radius_internal = 0.775
+    mesh = uw.meshing.SphericalShellInternalBoundary(
+        radiusOuter=radius_outer,
+        radiusInternal=radius_internal,
+        radiusInner=radius_inner,
+        cellSize=0.25,
+        qdegree=2,
+        degree=1,
+    )
+    velocity = uw.discretisation.MeshVariable("U_geoid", mesh, mesh.dim, degree=2)
+    pressure = uw.discretisation.MeshVariable(
+        "P_geoid", mesh, 1, degree=1, continuous=True
+    )
+    stokes = uw.systems.Stokes(
+        mesh,
+        velocityField=velocity,
+        pressureField=pressure,
+    )
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+
+    theta = mesh.CoordinateSystem.xR[1]
+    unit_r = mesh.CoordinateSystem.unit_e_0
+    load = sympy.assoc_legendre(2, 0, sympy.cos(theta)) * unit_r
+    stokes.add_natural_bc(load, "Internal")
+    stokes.add_rotated_freeslip_bc(0, "Upper", normal=unit_r)
+    stokes.add_rotated_freeslip_bc(0, "Lower", normal=-unit_r)
+    stokes.petsc_use_nullspace = True
+    stokes.petsc_options["snes_type"] = "ksponly"
+    stokes.tolerance = 1.0e-5
+    stokes.solve()
+
+    response = uw.postprocessing.zhong2008_response_from_rotated_stokes(
+        stokes=stokes,
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        radius_internal=radius_internal,
+        harmonic_degree=2,
+        load_scale=1.0,
+        include_self_gravity=True,
+    )
+
+    assert np.isclose(response.surface_topography, 0.41920, rtol=0.10)
+    assert np.isclose(response.cmb_topography, 0.77060, rtol=0.10)
+    assert np.isclose(response.surface_geoid, 0.02579, rtol=0.10)
+    assert np.isclose(response.cmb_geoid, 0.03206, rtol=0.10)
+    assert np.isclose(response.self_gravity.surface_topography, 0.49980, rtol=0.10)
+    assert np.isclose(response.self_gravity.cmb_topography, 0.93130, rtol=0.10)
+    assert np.isclose(response.self_gravity.surface_geoid, 0.04486, rtol=0.10)
+    assert np.isclose(response.self_gravity.cmb_geoid, 0.05461, rtol=0.10)

From d79627c07d3f3fe6f3760111ae3cbec0b5bac330 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Sun, 16 Aug 2026 20:39:15 +0530
Subject: [PATCH 2/3] Generalize spherical-shell geoid postprocessing

Move the public API under uw.postprocessing.geoid and remove Zhong-specific names from the reusable response data types and coefficient functions.

Support two-boundary shells with an optional internal load, require model-specific density and gravity inputs explicitly, and retain a focused rotated-free-slip adapter for recovering axisymmetric topography coefficients.

Rename the MPI regression, add no-load and parameter-validation coverage, and update developer documentation to distinguish generic postprocessing from a future semi-analytical propagator solver in uw.analytic.
---
 docs/developer/CHANGELOG.md                   |   9 +-
 ...ry-stress-and-projection-postprocessing.md |  45 +++-
 src/underworld3/postprocessing/__init__.py    |  18 +-
 src/underworld3/postprocessing/geoid.py       | 240 ++++++++++--------
 ...st_1071_spherical_shell_geoid_parallel.py} |  21 +-
 tests/test_1070_postprocessing_geoid.py       | 106 +++++---
 6 files changed, 265 insertions(+), 174 deletions(-)
 rename tests/parallel/{test_1071_zhong2008_geoid_parallel.py => test_1071_spherical_shell_geoid_parallel.py} (79%)

diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md
index 799e9b2a3..d417a24ad 100644
--- a/docs/developer/CHANGELOG.md
+++ b/docs/developer/CHANGELOG.md
@@ -273,10 +273,11 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
   seam sensitivity is unchanged (#404).
 - Recorded as the preferred free-slip BC in the project guidance (#300);
   conda PETSc floor raised to ≥ 3.25 for FMG/rotation API consistency (#304).
-- `uw.postprocessing.zhong2008_response_from_rotated_stokes(...)` projects the
-  existing rotated-free-slip boundary traction onto the benchmark harmonic and
-  applies one Appendix A operator for geoid and self-gravity response. The pure
-  operator functions also accept precomputed topography coefficients.
+- `uw.postprocessing.geoid` provides generic spherical-shell geoid and
+  self-gravity coefficient functions. Its rotated-Stokes adapter projects the
+  existing boundary traction onto an axisymmetric harmonic; the pure functions
+  also accept coefficients recovered by other methods and an optional internal
+  load.
 
 ### Generalized Geometric Multigrid via Custom Prolongation (July 2026)
 
diff --git a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
index 90126e2b2..cf2a1d77f 100644
--- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
+++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
@@ -74,37 +74,51 @@ Combine the two rules — project the τ components with `linear_solver()` for t
 cheap linear solve, then compose `σ_rr` analytically — for accurate boundary
 stress recovery that also scales.
 
-## 3. Zhong spherical-shell geoid response
+## 3. Spherical-shell geoid and self-gravity response
 
 Rotated free slip already exposes normal traction through
-`Stokes.boundary_normal_traction()`. The Zhong adapter projects that recovered
-traction onto the unnormalised `P_l^0` benchmark harmonic and applies the
-Appendix A geoid operator:
+`Stokes.boundary_normal_traction()`. The convenience adapter projects that
+recovered traction onto the unnormalised axisymmetric `P_l^0` harmonic and
+applies the spherical-shell geoid operator:
 
 ```python
 stokes.solve()
 
-response = uw.postprocessing.zhong2008_response_from_rotated_stokes(
+response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
     stokes=stokes,
     radius_inner=0.55,
     radius_outer=1.0,
-    radius_internal=0.775,
     harmonic_degree=2,
-    load_scale=1.0,
+    internal_load_radius=0.775,
+    internal_load_coefficient=1.0,
     include_self_gravity=True,
+    surface_density_contrast=3300.0,
+    cmb_density_contrast=5400.0,
+    planet_radius=6370000.0,
+    gravity=9.8,
+    gravitational_constant=6.67e-11,
 )
 ```
 
 The adapter delegates stress recovery to the existing rotated-free-slip API;
 it does not implement a second CBF, constrained-multiplier, or topography
-recovery path. `load_scale` must match the coefficient multiplying the model's
-internal `P_l^0` load.
+recovery path. `internal_load_coefficient` must use the same harmonic
+normalisation and sign convention as the model's internal load.
 
 When surface and CMB topography coefficients are already available, call
-`zhong2008_geoid_response()` or `zhong2008_self_gravity_response()` directly.
-These functions are pure post-processing and do not require a Stokes object.
+`uw.postprocessing.geoid.spherical_shell_geoid_response()` or
+`uw.postprocessing.geoid.spherical_shell_self_gravity_response()` directly.
+These functions are pure post-processing, work for any spherical-harmonic
+order with a consistent coefficient normalisation, and do not require a Stokes
+object. The internal load is optional.
 
-The Zhong Appendix A calculation is expressed as one linear operator:
+The density contrasts, dimensional outer-radius scale, and gravity are required
+when self-gravity is enabled. They deliberately have no Earth- or
+benchmark-specific defaults. The universal gravitational constant defaults to
+the current CODATA value and can be overridden when reproducing a paper's
+rounded constant.
+
+The calculation is expressed as one linear operator:
 
 ```text
 N = G h + n_load
@@ -116,10 +130,15 @@ responses, and `Q` contains the two self-gravity density factors. Sharing `G`
 and `n_load` between the no-self-gravity and self-gravity paths avoids separate
 scalar implementations of the same coefficients.
 
+This module does not compute a benchmark's semi-analytical Stokes solution.
+Published reference solvers, such as the Zhong et al. propagator-matrix method,
+belong in `uw.analytic`; their computed topography coefficients can be passed to
+the pure post-processing functions above.
+
 The rotated harmonic projector gathers boundary samples to rank zero and
 reconstructs their spherical triangulation. A future boundary-reaction
 functional could replace this step with a direct distributed finite-element
-projection without changing the Zhong operator API.
+projection without changing the coefficient API.
 
 ## See also
 
diff --git a/src/underworld3/postprocessing/__init__.py b/src/underworld3/postprocessing/__init__.py
index 57997e486..4e4325a1c 100644
--- a/src/underworld3/postprocessing/__init__.py
+++ b/src/underworld3/postprocessing/__init__.py
@@ -6,21 +6,5 @@
 """
 
 from . import geoid
-from .geoid import (
-    Zhong2008DynamicResponse,
-    Zhong2008GeoidResponse,
-    Zhong2008SelfGravityResponse,
-    zhong2008_geoid_response,
-    zhong2008_response_from_rotated_stokes,
-    zhong2008_self_gravity_response,
-)
 
-__all__ = [
-    "geoid",
-    "Zhong2008DynamicResponse",
-    "Zhong2008GeoidResponse",
-    "Zhong2008SelfGravityResponse",
-    "zhong2008_geoid_response",
-    "zhong2008_response_from_rotated_stokes",
-    "zhong2008_self_gravity_response",
-]
+__all__ = ["geoid"]
diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py
index a59cb328a..b37fcdcb4 100644
--- a/src/underworld3/postprocessing/geoid.py
+++ b/src/underworld3/postprocessing/geoid.py
@@ -1,4 +1,11 @@
-r"""Zhong et al. (2008) spherical-shell geoid response functions."""
+r"""Spherical-harmonic geoid and self-gravity response functions.
+
+The pure coefficient functions in this module are independent of a particular
+Stokes discretisation.  They combine surface, CMB, and optional internal-load
+coefficients through the radial Green's function for one spherical-harmonic
+degree.  A separate convenience adapter obtains the two topography coefficients
+from a completed rotated-free-slip Stokes solve.
+"""
 
 from __future__ import annotations
 
@@ -9,16 +16,26 @@
 import numpy as np
 
 
+__all__ = [
+    "GeoidResponse",
+    "SelfGravityResponse",
+    "SphericalShellResponse",
+    "spherical_shell_geoid_response",
+    "spherical_shell_self_gravity_response",
+    "spherical_shell_response_from_rotated_stokes",
+]
+
+
 @dataclass(frozen=True)
-class Zhong2008GeoidResponse:
-    """No-self-gravity surface and CMB geoid coefficients."""
+class GeoidResponse:
+    """Surface and CMB geoid coefficients without self-gravity feedback."""
 
     surface_geoid: float
     cmb_geoid: float
 
 
 @dataclass(frozen=True)
-class Zhong2008SelfGravityResponse:
+class SelfGravityResponse:
     """Self-gravity-corrected topography and geoid coefficients."""
 
     surface_topography: float
@@ -30,59 +47,66 @@ class Zhong2008SelfGravityResponse:
 
 
 @dataclass(frozen=True)
-class Zhong2008DynamicResponse:
-    """Rotated-Stokes topography and corresponding geoid coefficients."""
+class SphericalShellResponse:
+    """Spherical-shell topography and corresponding geoid coefficients."""
 
     surface_topography: float
     cmb_topography: float
     surface_geoid: float
     cmb_geoid: float
-    self_gravity: Zhong2008SelfGravityResponse | None = None
+    self_gravity: SelfGravityResponse | None = None
 
 
 def _validate_geometry(
     radius_inner: float,
     radius_outer: float,
-    radius_internal: float,
     harmonic_degree: int,
-) -> tuple[float, float, float, int]:
+) -> tuple[float, float, int]:
     try:
         ri = float(radius_inner)
         ro = float(radius_outer)
-        rint = float(radius_internal)
     except (TypeError, ValueError) as error:
         raise TypeError("The shell radii must be real numbers.") from error
-    if not np.all(np.isfinite((ri, rint, ro))) or not 0.0 < ri < rint < ro:
-        raise ValueError(
-            "Expected finite radii ordered as "
-            "0 < radius_inner < radius_internal < radius_outer."
-        )
+    if not np.all(np.isfinite((ri, ro))) or not 0.0 < ri < ro:
+        raise ValueError("Expected finite radii ordered as 0 < radius_inner < radius_outer.")
     if isinstance(harmonic_degree, bool) or not isinstance(harmonic_degree, Integral):
         raise TypeError("harmonic_degree must be an integer.")
     degree = int(harmonic_degree)
     if degree < 1:
         raise ValueError("harmonic_degree must be at least one.")
-    return ri, ro, rint, degree
+    return ri, ro, degree
 
 
-def _zhong2008_geoid_operator(
+def _spherical_shell_geoid_operator(
     radius_inner: float,
     radius_outer: float,
-    radius_internal: float,
     harmonic_degree: int,
-    load_scale: float,
+    internal_load_radius: float | None,
+    internal_load_coefficient: float,
 ) -> tuple[np.ndarray, np.ndarray]:
-    """Return the Zhong Appendix A topography operator and load vector."""
+    """Return the boundary-topography operator and fixed-load vector."""
 
-    ri, ro, rint, degree = _validate_geometry(
+    ri, ro, degree = _validate_geometry(
         radius_inner,
         radius_outer,
-        radius_internal,
         harmonic_degree,
     )
-    load_scale = float(load_scale)
-    if not np.isfinite(load_scale):
-        raise ValueError("load_scale must be finite.")
+    internal_load_coefficient = float(internal_load_coefficient)
+    if not np.isfinite(internal_load_coefficient):
+        raise ValueError("internal_load_coefficient must be finite.")
+
+    rint = None
+    if internal_load_radius is not None:
+        try:
+            rint = float(internal_load_radius)
+        except (TypeError, ValueError) as error:
+            raise TypeError("internal_load_radius must be a real number or None.") from error
+        if not np.isfinite(rint) or not ri < rint < ro:
+            raise ValueError("internal_load_radius must lie strictly between the shell radii.")
+    elif internal_load_coefficient != 0.0:
+        raise ValueError(
+            "internal_load_radius is required when internal_load_coefficient is nonzero."
+        )
 
     denominator = float(2 * degree + 1)
     operator = np.array(
@@ -98,38 +122,43 @@ def _zhong2008_geoid_operator(
         ],
         dtype=float,
     )
-    load = load_scale * np.array(
-        [
-            -rint * (rint / ro) ** (degree + 1) / denominator,
-            -rint * (ri / rint) ** degree / denominator,
-        ],
-        dtype=float,
-    )
+    load = np.zeros(2, dtype=float)
+    if rint is not None and internal_load_coefficient != 0.0:
+        load = -internal_load_coefficient * np.array(
+            [
+                rint * (rint / ro) ** (degree + 1) / denominator,
+                rint * (ri / rint) ** degree / denominator,
+            ],
+            dtype=float,
+        )
     return operator, load
 
 
-def zhong2008_geoid_response(
+def spherical_shell_geoid_response(
     *,
     radius_inner: float,
     radius_outer: float,
-    radius_internal: float,
     harmonic_degree: int,
     surface_topography_coefficient: float,
     cmb_topography_coefficient: float,
-    load_scale: float = 1.0,
-) -> Zhong2008GeoidResponse:
-    r"""Return the nondimensional Zhong no-self-gravity geoid coefficients.
-
-    The topography coefficients and ``load_scale`` must use the same
-    unnormalised :math:`P_l^0` forcing convention as the Stokes solve.
+    internal_load_radius: float | None = None,
+    internal_load_coefficient: float = 0.0,
+) -> GeoidResponse:
+    r"""Return spherical-shell geoid coefficients without self-gravity.
+
+    All source coefficients must use one consistent spherical-harmonic
+    normalisation.  The radial potential kernel depends on degree but not order,
+    so coefficients for any order :math:`m` may be used.  A positive optional
+    internal-load coefficient follows the outward radial-load convention and
+    enters the potential with the opposite sign to boundary topography.
     """
 
-    operator, load = _zhong2008_geoid_operator(
+    operator, load = _spherical_shell_geoid_operator(
         radius_inner,
         radius_outer,
-        radius_internal,
         harmonic_degree,
-        load_scale,
+        internal_load_radius,
+        internal_load_coefficient,
     )
     topography = np.array(
         [surface_topography_coefficient, cmb_topography_coefficient],
@@ -138,61 +167,53 @@ def zhong2008_geoid_response(
     if not np.all(np.isfinite(topography)):
         raise ValueError("The topography coefficients must be finite.")
     geoid = operator @ topography + load
-    return Zhong2008GeoidResponse(float(geoid[0]), float(geoid[1]))
+    return GeoidResponse(float(geoid[0]), float(geoid[1]))
 
 
-def zhong2008_self_gravity_response(
+def spherical_shell_self_gravity_response(
     *,
     radius_inner: float,
     radius_outer: float,
-    radius_internal: float,
     harmonic_degree: int,
     surface_topography_coefficient: float,
     cmb_topography_coefficient: float,
-    load_scale: float = 1.0,
-    surface_density_contrast: float = 3300.0,
-    cmb_density_contrast: float = 5400.0,
-    planet_radius: float = 6370000.0,
-    gravity: float = 9.8,
-    gravitational_constant: float = 6.67e-11,
-) -> Zhong2008SelfGravityResponse:
-    r"""Return the Zhong self-gravity-corrected response coefficients.
-
-    Density contrasts are in kg/m³, ``planet_radius`` in metres, ``gravity`` in
-    m/s², and ``gravitational_constant`` in SI units. Shell radii and response
-    coefficients remain nondimensional.
+    surface_density_contrast: float,
+    cmb_density_contrast: float,
+    planet_radius: float,
+    gravity: float,
+    internal_load_radius: float | None = None,
+    internal_load_coefficient: float = 0.0,
+    gravitational_constant: float = 6.67430e-11,
+) -> SelfGravityResponse:
+    r"""Return self-gravity-corrected spherical-shell response coefficients.
+
+    Density contrasts are signed and in kg/m³, ``planet_radius`` is in metres,
+    ``gravity`` is in m/s², and ``gravitational_constant`` is in SI units.
+    Shell radii and response coefficients remain nondimensional.
     """
 
-    operator, load = _zhong2008_geoid_operator(
+    operator, load = _spherical_shell_geoid_operator(
         radius_inner,
         radius_outer,
-        radius_internal,
         harmonic_degree,
-        load_scale,
+        internal_load_radius,
+        internal_load_coefficient,
     )
     topography = np.array(
         [surface_topography_coefficient, cmb_topography_coefficient],
         dtype=float,
     )
-    physical_values = np.array(
-        [
-            surface_density_contrast,
-            cmb_density_contrast,
-            planet_radius,
-            gravity,
-            gravitational_constant,
-        ],
-        dtype=float,
-    )
+    density_contrasts = np.array([surface_density_contrast, cmb_density_contrast], dtype=float)
+    physical_constants = np.array([planet_radius, gravity, gravitational_constant], dtype=float)
     if not np.all(np.isfinite(topography)):
         raise ValueError("The topography coefficients must be finite.")
-    if not np.all(np.isfinite(physical_values)) or np.any(physical_values <= 0.0):
-        raise ValueError("Density contrasts and physical constants must be positive.")
+    if not np.all(np.isfinite(density_contrasts)):
+        raise ValueError("Density contrasts must be finite.")
+    if not np.all(np.isfinite(physical_constants)) or np.any(physical_constants <= 0.0):
+        raise ValueError("Physical constants must be finite and positive.")
 
     gravity_scale = 4.0 * np.pi * gravitational_constant * planet_radius / gravity
-    q = gravity_scale * np.array(
-        [surface_density_contrast, cmb_density_contrast], dtype=float
-    )
+    q = gravity_scale * density_contrasts
     q_matrix = np.diag(q)
     corrected_topography = np.linalg.solve(
         np.eye(2) - q_matrix @ operator,
@@ -200,7 +221,7 @@ def zhong2008_self_gravity_response(
     )
     corrected_geoid = operator @ corrected_topography + load
 
-    return Zhong2008SelfGravityResponse(
+    return SelfGravityResponse(
         surface_topography=float(corrected_topography[0]),
         cmb_topography=float(corrected_topography[1]),
         surface_geoid=float(corrected_geoid[0]),
@@ -300,41 +321,58 @@ def _rotated_topography_coefficient(
     return float(MPI.COMM_WORLD.bcast(coefficient, root=0))
 
 
-def zhong2008_response_from_rotated_stokes(
+def spherical_shell_response_from_rotated_stokes(
     *,
     stokes,
     radius_inner: float,
     radius_outer: float,
-    radius_internal: float,
     harmonic_degree: int,
-    load_scale: float = 1.0,
+    internal_load_radius: float | None = None,
+    internal_load_coefficient: float = 0.0,
     surface_boundary: str = "Upper",
     cmb_boundary: str = "Lower",
     surface_buoyancy_scale: float = 1.0,
     cmb_buoyancy_scale: float = 1.0,
     include_self_gravity: bool = False,
-    surface_density_contrast: float = 3300.0,
-    cmb_density_contrast: float = 5400.0,
-    planet_radius: float = 6370000.0,
-    gravity: float = 9.8,
-    gravitational_constant: float = 6.67e-11,
-) -> Zhong2008DynamicResponse:
-    r"""Compute Zhong response coefficients from a rotated-free-slip Stokes solve.
+    surface_density_contrast: float | None = None,
+    cmb_density_contrast: float | None = None,
+    planet_radius: float | None = None,
+    gravity: float | None = None,
+    gravitational_constant: float = 6.67430e-11,
+) -> SphericalShellResponse:
+    r"""Compute spherical-shell response from a rotated-free-slip Stokes solve.
 
     Normal traction recovery is delegated to the existing
     :meth:`Stokes.boundary_normal_traction` implementation. This adapter only
     projects the two boundary responses onto the unnormalised
-    :math:`P_l^0` harmonic and applies the Zhong Appendix A operator.
+    axisymmetric :math:`P_l^0` harmonic.  Use the pure coefficient functions
+    directly for other harmonic orders or topography-recovery methods. Density
+    contrasts, planet radius, and gravity are required when
+    ``include_self_gravity`` is true.
     """
 
-    ri, ro, _, degree = _validate_geometry(
+    ri, ro, degree = _validate_geometry(
         radius_inner,
         radius_outer,
-        radius_internal,
         harmonic_degree,
     )
     if not isinstance(include_self_gravity, bool):
         raise TypeError("include_self_gravity must be True or False.")
+    if include_self_gravity:
+        missing = [
+            name
+            for name, value in (
+                ("surface_density_contrast", surface_density_contrast),
+                ("cmb_density_contrast", cmb_density_contrast),
+                ("planet_radius", planet_radius),
+                ("gravity", gravity),
+            )
+            if value is None
+        ]
+        if missing:
+            raise ValueError(
+                "Self-gravity requires explicit values for " + ", ".join(missing) + "."
+            )
 
     surface_topography = _rotated_topography_coefficient(
         stokes,
@@ -352,33 +390,33 @@ def zhong2008_response_from_rotated_stokes(
         cmb_buoyancy_scale,
         -1.0,
     )
-    geoid = zhong2008_geoid_response(
+    geoid = spherical_shell_geoid_response(
         radius_inner=ri,
         radius_outer=ro,
-        radius_internal=radius_internal,
         harmonic_degree=degree,
         surface_topography_coefficient=surface_topography,
         cmb_topography_coefficient=cmb_topography,
-        load_scale=load_scale,
+        internal_load_radius=internal_load_radius,
+        internal_load_coefficient=internal_load_coefficient,
     )
     self_gravity = None
     if include_self_gravity:
-        self_gravity = zhong2008_self_gravity_response(
+        self_gravity = spherical_shell_self_gravity_response(
             radius_inner=ri,
             radius_outer=ro,
-            radius_internal=radius_internal,
             harmonic_degree=degree,
             surface_topography_coefficient=surface_topography,
             cmb_topography_coefficient=cmb_topography,
-            load_scale=load_scale,
-            surface_density_contrast=surface_density_contrast,
-            cmb_density_contrast=cmb_density_contrast,
-            planet_radius=planet_radius,
-            gravity=gravity,
+            internal_load_radius=internal_load_radius,
+            internal_load_coefficient=internal_load_coefficient,
+            surface_density_contrast=float(surface_density_contrast),
+            cmb_density_contrast=float(cmb_density_contrast),
+            planet_radius=float(planet_radius),
+            gravity=float(gravity),
             gravitational_constant=gravitational_constant,
         )
 
-    return Zhong2008DynamicResponse(
+    return SphericalShellResponse(
         surface_topography=surface_topography,
         cmb_topography=cmb_topography,
         surface_geoid=geoid.surface_geoid,
diff --git a/tests/parallel/test_1071_zhong2008_geoid_parallel.py b/tests/parallel/test_1071_spherical_shell_geoid_parallel.py
similarity index 79%
rename from tests/parallel/test_1071_zhong2008_geoid_parallel.py
rename to tests/parallel/test_1071_spherical_shell_geoid_parallel.py
index a5284d97b..189d61217 100644
--- a/tests/parallel/test_1071_zhong2008_geoid_parallel.py
+++ b/tests/parallel/test_1071_spherical_shell_geoid_parallel.py
@@ -7,25 +7,23 @@
 pytestmark = pytest.mark.level_2
 
 
-def test_zhong2008_rotated_geoid_matches_serial_reference():
+def test_rotated_spherical_shell_geoid_matches_serial_reference():
     if uw.mpi.size == 1:
         pytest.skip("Run this regression with at least two MPI ranks.")
 
     radius_inner = 0.55
     radius_outer = 1.0
-    radius_internal = 0.775
+    rint = 0.775
     mesh = uw.meshing.SphericalShellInternalBoundary(
         radiusOuter=radius_outer,
-        radiusInternal=radius_internal,
+        radiusInternal=rint,
         radiusInner=radius_inner,
         cellSize=0.25,
         qdegree=2,
         degree=1,
     )
     velocity = uw.discretisation.MeshVariable("U_geoid_mpi", mesh, mesh.dim, degree=2)
-    pressure = uw.discretisation.MeshVariable(
-        "P_geoid_mpi", mesh, 1, degree=1, continuous=True
-    )
+    pressure = uw.discretisation.MeshVariable("P_geoid_mpi", mesh, 1, degree=1, continuous=True)
     stokes = uw.systems.Stokes(
         mesh,
         velocityField=velocity,
@@ -45,14 +43,19 @@ def test_zhong2008_rotated_geoid_matches_serial_reference():
     stokes.tolerance = 1.0e-5
     stokes.solve()
 
-    response = uw.postprocessing.zhong2008_response_from_rotated_stokes(
+    response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
         stokes=stokes,
         radius_inner=radius_inner,
         radius_outer=radius_outer,
-        radius_internal=radius_internal,
         harmonic_degree=2,
-        load_scale=1.0,
+        internal_load_radius=rint,
+        internal_load_coefficient=1.0,
         include_self_gravity=True,
+        surface_density_contrast=3300.0,
+        cmb_density_contrast=5400.0,
+        planet_radius=6370000.0,
+        gravity=9.8,
+        gravitational_constant=6.67e-11,
     )
     values = np.array(
         [
diff --git a/tests/test_1070_postprocessing_geoid.py b/tests/test_1070_postprocessing_geoid.py
index 9f8b794b3..90375027c 100644
--- a/tests/test_1070_postprocessing_geoid.py
+++ b/tests/test_1070_postprocessing_geoid.py
@@ -9,63 +9,84 @@
 pytestmark = pytest.mark.level_2
 
 
-def test_zhong2008_postprocessing_is_public():
+def test_spherical_shell_postprocessing_is_public():
     assert hasattr(uw, "postprocessing")
-    assert hasattr(uw.postprocessing, "zhong2008_geoid_response")
-    assert hasattr(uw.postprocessing, "zhong2008_self_gravity_response")
-    assert hasattr(uw.postprocessing, "zhong2008_response_from_rotated_stokes")
+    assert uw.postprocessing.__all__ == ["geoid"]
+    assert hasattr(uw.postprocessing.geoid, "spherical_shell_geoid_response")
+    assert hasattr(uw.postprocessing.geoid, "spherical_shell_self_gravity_response")
+    assert hasattr(uw.postprocessing.geoid, "spherical_shell_response_from_rotated_stokes")
 
 
-def test_zhong2008_geoid_response_matches_direct_formula_and_load_scale():
+def test_spherical_shell_geoid_response_matches_direct_formula_and_load():
     radius_inner = 0.55
     radius_outer = 1.0
-    radius_internal = 0.775
+    rint = 0.775
     degree = 2
-    load_scale = 1.75
+    internal_load_coefficient = 1.75
     surface_topography = 0.3918264884592169
     cmb_topography = 0.2653834107104151
 
-    response = uw.postprocessing.zhong2008_geoid_response(
+    response = uw.postprocessing.geoid.spherical_shell_geoid_response(
         radius_inner=radius_inner,
         radius_outer=radius_outer,
-        radius_internal=radius_internal,
         harmonic_degree=degree,
         surface_topography_coefficient=surface_topography,
         cmb_topography_coefficient=cmb_topography,
-        load_scale=load_scale,
+        internal_load_radius=rint,
+        internal_load_coefficient=internal_load_coefficient,
     )
 
     denominator = 2 * degree + 1
     surface_expected = (
         radius_inner * (radius_inner / radius_outer) ** (degree + 1) * cmb_topography
         + radius_outer * surface_topography
-        - load_scale
-        * radius_internal
-        * (radius_internal / radius_outer) ** (degree + 1)
+        - internal_load_coefficient * rint * (rint / radius_outer) ** (degree + 1)
     ) / denominator
     cmb_expected = (
         radius_inner * cmb_topography
         + radius_outer * (radius_inner / radius_outer) ** degree * surface_topography
-        - load_scale * radius_internal * (radius_inner / radius_internal) ** degree
+        - internal_load_coefficient * rint * (radius_inner / rint) ** degree
     ) / denominator
 
     assert math.isclose(response.surface_geoid, surface_expected)
     assert math.isclose(response.cmb_geoid, cmb_expected)
 
 
-def test_zhong2008_self_gravity_response_matches_direct_solve():
+def test_spherical_shell_geoid_response_without_internal_load():
+    response = uw.postprocessing.geoid.spherical_shell_geoid_response(
+        radius_inner=0.55,
+        radius_outer=1.0,
+        harmonic_degree=3,
+        surface_topography_coefficient=0.4,
+        cmb_topography_coefficient=-0.2,
+    )
+
+    denominator = 7.0
+    expected = np.array(
+        [
+            (0.4 + 0.55 * 0.55**4 * -0.2) / denominator,
+            (0.55**3 * 0.4 + 0.55 * -0.2) / denominator,
+        ]
+    )
+    assert np.allclose([response.surface_geoid, response.cmb_geoid], expected)
+
+
+def test_spherical_shell_self_gravity_response_matches_direct_solve():
     surface_topography = 0.3918264884592169
     cmb_topography = 0.2653834107104151
-    response = uw.postprocessing.zhong2008_self_gravity_response(
+    response = uw.postprocessing.geoid.spherical_shell_self_gravity_response(
         radius_inner=0.55,
         radius_outer=1.0,
-        radius_internal=0.775,
         harmonic_degree=2,
         surface_topography_coefficient=surface_topography,
         cmb_topography_coefficient=cmb_topography,
-        load_scale=1.0,
+        internal_load_radius=0.775,
+        internal_load_coefficient=1.0,
         surface_density_contrast=3300.0,
         cmb_density_contrast=5400.0,
+        planet_radius=6370000.0,
+        gravity=9.8,
+        gravitational_constant=6.67e-11,
     )
 
     q_surface = 4.0 * np.pi * 6.67e-11 * 6370000.0 * 3300.0 / 9.8
@@ -103,35 +124,55 @@ def test_zhong2008_self_gravity_response_matches_direct_solve():
 
 
 @pytest.mark.parametrize("harmonic_degree", [0, -1, 2.0, True])
-def test_zhong2008_geoid_rejects_invalid_harmonic_degree(harmonic_degree):
+def test_spherical_shell_geoid_rejects_invalid_harmonic_degree(harmonic_degree):
     error = TypeError if isinstance(harmonic_degree, (float, bool)) else ValueError
     with pytest.raises(error):
-        uw.postprocessing.zhong2008_geoid_response(
+        uw.postprocessing.geoid.spherical_shell_geoid_response(
             radius_inner=0.55,
             radius_outer=1.0,
-            radius_internal=0.775,
             harmonic_degree=harmonic_degree,
             surface_topography_coefficient=0.4,
             cmb_topography_coefficient=0.7,
         )
 
 
-def test_zhong2008_rotated_stokes_adapter_matches_table_2():
+def test_internal_load_requires_a_radius():
+    with pytest.raises(ValueError, match="internal_load_radius is required"):
+        uw.postprocessing.geoid.spherical_shell_geoid_response(
+            radius_inner=0.55,
+            radius_outer=1.0,
+            harmonic_degree=2,
+            surface_topography_coefficient=0.4,
+            cmb_topography_coefficient=0.7,
+            internal_load_coefficient=1.0,
+        )
+
+
+def test_rotated_adapter_requires_explicit_self_gravity_parameters():
+    with pytest.raises(ValueError, match="Self-gravity requires explicit values"):
+        uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
+            stokes=object(),
+            radius_inner=0.55,
+            radius_outer=1.0,
+            harmonic_degree=2,
+            include_self_gravity=True,
+        )
+
+
+def test_rotated_stokes_adapter_matches_zhong_table_2():
     radius_inner = 0.55
     radius_outer = 1.0
-    radius_internal = 0.775
+    rint = 0.775
     mesh = uw.meshing.SphericalShellInternalBoundary(
         radiusOuter=radius_outer,
-        radiusInternal=radius_internal,
+        radiusInternal=rint,
         radiusInner=radius_inner,
         cellSize=0.25,
         qdegree=2,
         degree=1,
     )
     velocity = uw.discretisation.MeshVariable("U_geoid", mesh, mesh.dim, degree=2)
-    pressure = uw.discretisation.MeshVariable(
-        "P_geoid", mesh, 1, degree=1, continuous=True
-    )
+    pressure = uw.discretisation.MeshVariable("P_geoid", mesh, 1, degree=1, continuous=True)
     stokes = uw.systems.Stokes(
         mesh,
         velocityField=velocity,
@@ -151,14 +192,19 @@ def test_zhong2008_rotated_stokes_adapter_matches_table_2():
     stokes.tolerance = 1.0e-5
     stokes.solve()
 
-    response = uw.postprocessing.zhong2008_response_from_rotated_stokes(
+    response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
         stokes=stokes,
         radius_inner=radius_inner,
         radius_outer=radius_outer,
-        radius_internal=radius_internal,
         harmonic_degree=2,
-        load_scale=1.0,
+        internal_load_radius=rint,
+        internal_load_coefficient=1.0,
         include_self_gravity=True,
+        surface_density_contrast=3300.0,
+        cmb_density_contrast=5400.0,
+        planet_radius=6370000.0,
+        gravity=9.8,
+        gravitational_constant=6.67e-11,
     )
 
     assert np.isclose(response.surface_topography, 0.41920, rtol=0.10)

From 47b96006f7c6a7ddbae92624274d88abb46dbddc Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Sun, 16 Aug 2026 20:47:54 +0530
Subject: [PATCH 3/3] Separate pure and rotated harmonic degree constraints

Allow the generic spherical-shell coefficient functions to evaluate the degree-zero radial potential while keeping the rotated-Stokes adapter at degree one or greater because boundary normal traction recovery removes its mean.

Add focused degree-zero formula and adapter validation tests and document the distinction.
---
 ...ry-stress-and-projection-postprocessing.md |  4 ++-
 src/underworld3/postprocessing/geoid.py       |  9 +++++--
 tests/test_1070_postprocessing_geoid.py       | 25 ++++++++++++++++++-
 3 files changed, 34 insertions(+), 4 deletions(-)

diff --git a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
index cf2a1d77f..fb007f0ba 100644
--- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
+++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
@@ -110,7 +110,9 @@ When surface and CMB topography coefficients are already available, call
 `uw.postprocessing.geoid.spherical_shell_self_gravity_response()` directly.
 These functions are pure post-processing, work for any spherical-harmonic
 order with a consistent coefficient normalisation, and do not require a Stokes
-object. The internal load is optional.
+object. They support non-negative harmonic degrees, and the internal load is
+optional. The rotated-Stokes adapter requires degree one or greater because
+normal-traction recovery removes the degree-zero mean.
 
 The density contrasts, dimensional outer-radius scale, and gravity are required
 when self-gravity is enabled. They deliberately have no Earth- or
diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py
index b37fcdcb4..537a79fee 100644
--- a/src/underworld3/postprocessing/geoid.py
+++ b/src/underworld3/postprocessing/geoid.py
@@ -72,8 +72,8 @@ def _validate_geometry(
     if isinstance(harmonic_degree, bool) or not isinstance(harmonic_degree, Integral):
         raise TypeError("harmonic_degree must be an integer.")
     degree = int(harmonic_degree)
-    if degree < 1:
-        raise ValueError("harmonic_degree must be at least one.")
+    if degree < 0:
+        raise ValueError("harmonic_degree must be non-negative.")
     return ri, ro, degree
 
 
@@ -358,6 +358,11 @@ def spherical_shell_response_from_rotated_stokes(
     )
     if not isinstance(include_self_gravity, bool):
         raise TypeError("include_self_gravity must be True or False.")
+    if degree == 0:
+        raise ValueError(
+            "The rotated-Stokes adapter requires harmonic_degree >= 1 because "
+            "boundary_normal_traction() removes the degree-zero mean."
+        )
     if include_self_gravity:
         missing = [
             name
diff --git a/tests/test_1070_postprocessing_geoid.py b/tests/test_1070_postprocessing_geoid.py
index 90375027c..dce646a86 100644
--- a/tests/test_1070_postprocessing_geoid.py
+++ b/tests/test_1070_postprocessing_geoid.py
@@ -123,7 +123,20 @@ def test_spherical_shell_self_gravity_response_matches_direct_solve():
     )
 
 
-@pytest.mark.parametrize("harmonic_degree", [0, -1, 2.0, True])
+def test_spherical_shell_geoid_supports_degree_zero():
+    response = uw.postprocessing.geoid.spherical_shell_geoid_response(
+        radius_inner=0.55,
+        radius_outer=1.0,
+        harmonic_degree=0,
+        surface_topography_coefficient=0.4,
+        cmb_topography_coefficient=-0.2,
+    )
+
+    assert math.isclose(response.surface_geoid, 0.4 - 0.55**2 * 0.2)
+    assert math.isclose(response.cmb_geoid, 0.4 - 0.55 * 0.2)
+
+
+@pytest.mark.parametrize("harmonic_degree", [-1, 2.0, True])
 def test_spherical_shell_geoid_rejects_invalid_harmonic_degree(harmonic_degree):
     error = TypeError if isinstance(harmonic_degree, (float, bool)) else ValueError
     with pytest.raises(error):
@@ -136,6 +149,16 @@ def test_spherical_shell_geoid_rejects_invalid_harmonic_degree(harmonic_degree):
         )
 
 
+def test_rotated_adapter_rejects_degree_zero():
+    with pytest.raises(ValueError, match="requires harmonic_degree >= 1"):
+        uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
+            stokes=object(),
+            radius_inner=0.55,
+            radius_outer=1.0,
+            harmonic_degree=0,
+        )
+
+
 def test_internal_load_requires_a_radius():
     with pytest.raises(ValueError, match="internal_load_radius is required"):
         uw.postprocessing.geoid.spherical_shell_geoid_response(