`image_2layer` — Tagg / Sunde series for 2-layer soil¶

Physical context¶

The 2-layer model is the canonical typical layered-soil description: a finite upper layer of thickness $h_1$ and resistivity $\rho_1$ over a semi-infinite lower layer of resistivity $\rho_2$. It is the simplest soil that captures the dominant first-order effect of heterogeneity — a frozen, weathered, or saturated surface layer above a different sub-stratum — without introducing the identification problem that more layers create. A typical parameter range is $\rho_1 \in [50, 1000]\,\Omega\, \text{m}$, $\rho_2 \in [10, 5000]\,\Omega\,\text{m}$, $h_1 \in [0.5, 5]\,\text{m}$.

The single design degree of freedom that controls the layered behaviour is the interface reflection coefficient

\[ K_1 \;=\; \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} \in (-1, +1). \]

For $\rho_2 < \rho_1$ ($K_1 < 0$, conductive bottom) the lower layer pulls current away from the surface and the cluster impedance falls. For $\rho_2 > \rho_1$ ($K_1 > 0$, resistive bottom) the opposite holds. The limiting cases $K_1 \to \pm 1$ correspond to a perfect insulator below or a perfect conductor below.

Governing equation: Tagg / Sunde geometric series¶

A point current source $I$ at depth $z_s \in (0, h_1)$ produces a potential field that satisfies the Neumann boundary at $z = 0$ and the continuity-of-current boundary at $z = h_1$. The spectral-domain solution involves the Sommerfeld integral

\[ \varphi(s, z) \;=\; \frac{\rho_1\, I}{4\pi} \int_0^{\infty} \frac{ e^{-\lambda |z-z_s|} + e^{-\lambda (z+z_s)} + K_1\, e^{-\lambda (2 h_1 - z - z_s)} + K_1\, e^{-\lambda (2 h_1 - |z-z_s|)} } { 1 - K_1\, e^{-2\lambda h_1} } \, J_0(\lambda s)\, d\lambda. \]

Because $\Gamma_1(\lambda) \equiv K_1$ is constant in $\lambda$ (only one interface contributes) the multiple-reflection multiplier $1 / (1 - K_1 e^{-2\lambda h_1})$ expands as a geometric series in $K_1 e^{-2\lambda h_1}$, and the Sommerfeld integral collapses into a closed-form sum of point sources at mirrored positions.

After substitution and simplification, the spatial form is the classical Tagg / Sunde series:

\[ \varphi(s, z) \;=\; \frac{\rho_1\, I}{4\pi} \Biggl[ \tfrac{1}{r_0^+} + \tfrac{1}{r_0^-} + \sum_{n=1}^{\infty} K_1^n \Bigl(\tfrac{1}{r_n^{++}} + \tfrac{1}{r_n^{+-}} + \tfrac{1}{r_n^{-+}} + \tfrac{1}{r_n^{--}}\Bigr) \Biggr], \]

with the image distances

\[ r_n^{\sigma\tau} \;=\; \sqrt{(x{-}x_s)^2 + (y{-}y_s)^2 + (z - \sigma\,2 n h_1 - \tau\,z_s)^2}, \qquad \sigma, \tau \in \{+1, -1\}. \]

Per index $n$ four images are placed at $z_n^{\sigma\tau} = \sigma\,2 n h_1 + \tau\,z_s$ with weight $K_1^n$. The $n=0$ term contributes two images (the original source at $z = +z_s$ and the air mirror at $z = -z_s$, both with weight 1) — this is exactly the homogeneous image-charge solution.

Numerical strategy¶

Series truncation¶

The series is truncated at the smallest index for which the geometric tail bound falls below the tolerance,

\[ \frac{|K_1|^{n+1}}{1 - |K_1|} < \text{tol} \qquad (\text{default } \text{tol} = 10^{-6}), \]

or at max_terms (Engine.image_max_terms, default $300$), whichever happens first. The tail bound is the exact geometric remainder of the weight series $\sum_{m>n} |K_1|^m$; the historic per-term test $|K_1|^n < \text{tol}$ underestimated the neglected remainder by the factor $1/(1-|K_1|)$ — a factor 17 at $|K_1| = 0.94$. The default max_terms covers $|K_1| = 0.94$ (277 terms); harder contrasts converge more slowly and the per-segment cost grows linearly with the required number of image terms.

When max_terms is reached without meeting the tolerance, FieldResult.metadata['converged'] is set to False, a warning is logged and a SeriesTruncationWarning is emitted (silence with warnings.simplefilter("ignore", SeriesTruncationWarning) or raise Engine.image_max_terms).

Self-action¶

The wire-segment self-action splits cleanly:

The $n = 0$ direct contribution carries the singular line self-potential — handled by the homogeneous formula $2 \ln(L_i / a_i) / L_i$.
The $n = 0$ air-mirror term is a point image at distance $2 z_s$ and all $n \ge 1$ image terms are at least $2 h_1 - z - z_s$ away from the segment midpoint, so the point-source approximation is safe as long as the air mirror is well separated from the conductor, i.e. $2 z \gg a$ and $L_i \le 4 z_i$.

This is the same construction as in the homogeneous backend, with $\rho \to \rho_1$ as the prefactor and the geometric image series added on top — including the burial-depth guards. Since 0.15.0:

a segment midpoint with $z \le 0$ or $2|z| \le a$ raises ValueError, because its mirror image coincides with the conductor and the self-image term has no finite value (a ring at $z = 0$ used to return 164 Ω against a physical limit of ~11 Ω). Both conditions are enforced in the shared reaction-matrix kernel solver.image._self_corrected_kernel, which this backend routes its $n = 0$ term through, so mom / bem / cim / mom_sommerfeld / mutual reject the same geometries;
a segment whose extent crosses the surface ($z_i - \tfrac{1}{2} L_i |e_{z,i}| < 0$, i.e. midpoint buried but upper end in the air) also raises ValueError. This needs the segment tangent, so it is checked here and in image but not in the shared kernel;
$L_i > 4 z_i$ emits a ShallowSegmentWarning quantifying the bias with the horizontal-segment formula. It fires for plain electrodes too — a ring of radius 25 m at $z = 0.8$ m with $\Delta s = 5$ m warns "by 26 %" — and it is not emitted by the other backends of the family even though they share the bias.

See image for the derivation, the proof that the horizontal formula is the applicable one wherever the warning can still fire, and the honest statement of the residual bias (including the collocation-versus-Galerkin distinction).

Reported node potential¶

Electrode potentials, grounding_impedance and cluster_impedance are reduced with the length-weighted (Galerkin) average of the segment potentials, identical to the homogeneous backend and to the row reduction of the multi-port matrix inside the solve. This matters here in particular: the ADR-0007 interface split gives a rod crossing $h_1$ segments of unequal length, and up to 0.14.1 the reported value was a plain segment mean — measured $-10.9\,\%$ against the impedance the solver enforced for a 4.509 m rod split at $h_1 = 5$ m, and different reported potentials for two electrodes bonded into one galvanic node.

Auto-dispatch¶

Engine.solve picks image_2layer automatically when the user passed backend="image" and world.soil is a TwoLayerSoil. A client code written for the homogeneous backend therefore keeps working when only the soil model is swapped — no string change required.

Validity envelope¶

Property	Range / value
Soil model	`TwoLayerSoil` only
Frequency	quasi-static, $f < 1\,\text{kHz}$
Electrode placement	every segment must have $z < h_1$ (raises `ValueError` otherwise)
Burial depth (midpoints)	$z > 0$ and $2 z > a$ for every segment midpoint — hard `ValueError` otherwise, raised in the shared reaction-matrix kernel so `mom` / `bem` / `cim` / `mom_sommerfeld` / `mutual` reject the same geometries
Burial depth (extent)	$z_i - \tfrac{1}{2} L_i \lvert e_{z,i}\rvert \ge 0$: no segment may cross the surface. Hard `ValueError`, needs the segment tangent and therefore fires in `image` / `image_2layer` only. An upper end exactly at $z = 0$ is legal (ordinary driven rod)
Self-image bias	$L_i \le 4 z_i$ for an unbiased self-image term; `ShallowSegmentWarning` otherwise, quantified with the horizontal-segment formula. Emitted by `image` / `image_2layer` only although the bias is shared by the whole family
Reflection coefficient	$
Series truncation	adaptive tail bound $
Wire radius / segment ratio	thin-wire, $a \ll L_i$

The hard constraint that all segments lie inside the upper layer is a direct consequence of the series construction: the image positions $\pm 2 n h_1 \pm z_s$ are derived assuming the source is above the interface. A segment at $z = h_1$ would put a source on the interface and break the geometric expansion.

Convergence and cost¶

Series convergence. Geometric in $|K_1|$ (terms for the $10^{-6}$ tail bound):
$|K_1| = 0.5$ → 20 terms.
$|K_1| = 0.8$ → 69 terms.
$|K_1| = 0.9$ → 153 terms.
$|K_1| = 0.94$ → 277 terms (just inside the default max_terms = 300).
$|K_1| = 0.95$ → 328 terms (above the default — raises a SeriesTruncationWarning; increase Engine.image_max_terms).
Per-segment cost. $O(N^2 \cdot M)$ kernel evaluations with $M$ the truncated series length. For typical cases with $N \le 1000$ and $M \le 200$ the engine still runs in seconds.
Reduction to homogeneous. At $K_1 = 0$ the engine collapses to the homogeneous backend bit-exactly — the only series term with non-zero weight is $n = 0$ with weight $1$, identical to the image backend with $\rho = \rho_1$.

Cross-validation notes¶

Counterpart	Expected agreement	What is checked
`image` at $K_1 = 0$	bit-exact	series collapse
`mom` on the same world	$\le 2\,\%$	same kernel, different resolution scheme
`cim` on the same world	bit-exact (CIM falls back on this kernel for $n = 2$)	`cim` deliberately reuses the Tagg / Sunde self-kernel for $n = 2$
`bem` on the same world	bit-exact (BEM falls back on this kernel for $n = 2$)	same reason
`mom_sommerfeld` on the same world	$\le 5\,\%$	quadrature reproduces the geometric series
Sunde / Dwight closed forms	$\le 10\,\%$	reference values for canonical geometries (rods, rings)

The engine's role inside the cross-validation matrix is to be the closed-form anchor: every other engine is checked against it for $n = 2$.

References¶

Tagg, G. F. (1964). Earth Resistances, Pitman, ch. 5. The classical statement of the geometric image series.
Sunde, E. D. (1968). Earth Conduction Effects in Transmission Systems, Dover, sect. 3.5. Derivation of the multi-layer Sommerfeld integral and the 2-layer reduction.
Stefanescu, S. & Schlumberger, C. (1930). Sur la distribution électrique potentielle autour d'une prise de terre ponctuelle, Journal de Physique. Original geophysical formulation of the image-charge series.
Dawalibi, F. P. & Barbeito, N. (1991). Measurements and computations of the performance of grounding systems buried in multilayer soils. IEEE PWRD 6(4) — extension to $n$ layers and modern computational practice.
Colominas, I., París, J., Navarrina, F. & Casteleiro, M. (2012). Improvement of computer methods for grounding analysis in layered soils by using high-efficient convergence acceleration techniques. Adv. Eng. Soft. 44 — Aitken / Pade acceleration of the series for $|K_1| \to 1$.

Example¶

import groundfield as gf

soil = gf.TwoLayerSoil(rho_1=100.0, rho_2=500.0, h_1=2.0)
world = gf.create_world(soil=soil)
gf.create_electrode(world, "ring", name="g1",
                    center=(0.0, 0.0, 0.8), radius=5.0,
                    wire_radius=0.005)
gf.create_source(world, attached_to="g1", magnitude=1.0)

engine = gf.create_engine(backend="image_2layer",
                          segment_length=0.1,
                          frequencies=[50.0])
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])

API reference¶

image_2layer ¶

Image-charge backend for 2-layer soil (Tagg / Sunde).

Computes the potential field in a layered half-space (upper layer $\rho_1$ of thickness $h_1$, semi-infinite lower layer $\rho_2$) using the classical image-charge series.

Theory (short)

With the reflection coefficient at the layer interface $$ K \;=\; \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} $$ a point current source $I$ at $(x_s, y_s, z_s)$ with $0 < z_s < h_1$ produces, at any field point in the upper layer, the potential $$ \varphi(x, y, z) \;=\; \frac{\rho_1\, I}{4\pi}\,\Bigl[ \tfrac{1}{r_0^+} + \tfrac{1}{r_0^-} + \sum_{n=1}^{\infty} K^n \Bigl(\tfrac{1}{r_n^{++}} + \tfrac{1}{r_n^{+-}} + \tfrac{1}{r_n^{-+}} + \tfrac{1}{r_n^{--}}\Bigr)\Bigr], $$ with the image-source distances $$ r_n^{\sigma\tau} \;=\; \sqrt{(x-x_s)^2 + (y-y_s)^2 + (z - \sigma\,2 n h_1 - \tau z_s)^2}, \qquad \sigma, \tau \in {+1,-1}. $$ The $n = 0$ term (two images, weight 1) is exactly the homogeneous image-charge backend. The 2-layer backend therefore reduces to the homogeneous result for $\rho_2 = \rho_1$ (i.e. $K = 0$) — usable as a consistency test.

Convergence

Series terms decay as $K^n$ (geometric). For $|K| < 1$, typically $n \lesssim 50$ is enough for a relative accuracy of $10^{-6}$. The implementation truncates after max_terms terms or as soon as $|K|^n < \text{tol}$. If max_terms is hit without reaching the tolerance, a warning is recorded in FieldResult.metadata.

Preconditions

All electrodes must lie completely in the upper layer ($z_\text{seg} < h_1$); otherwise the backend raises a clear ValueError.
For $|K| \to 1$ (extreme contrasts, e.g. wet soil over rock) convergence slows down significantly. The solver issues a log warning in that case.

References

Tagg, G. F. (1964). Earth Resistances. Pitman, ch. 5.
Sunde, E. D. (1968). Earth Conduction Effects in Transmission Systems. Dover, sect. 3.5.
ADR-0001 docs/adr/0001-two-layer-method.md.

SeriesTruncationWarning ¶

Bases: UserWarning

The Tagg/Sunde image series hit max_terms before tol.

Emitted (in addition to the logger record) whenever a backend truncates the two-layer image-charge series without reaching the requested tail tolerance — typically at high layer contrast (:math:|K| \to 1). Silence with warnings.simplefilter("ignore", SeriesTruncationWarning) or raise Engine.image_max_terms.

solve_image_2layer ¶

solve_image_2layer(
    world: "World",
    engine: "Engine",
    *,
    max_terms: int = 100,
    tol: float = 1e-06
) -> FieldResult

Image-charge solver for 2-layer soil (Tagg / Sunde).

Parameters:

Name	Type	Description	Default
`world`	`'World'`	World whose `soil` is a :class:`TwoLayerSoil`. All electrodes must lie inside the upper layer.	required
`engine`	`'Engine'`	Engine configuration. `engine.segment_length` controls the discretisation as in the homogeneous backend.	required
`max_terms`	`int`	Maximum number of series terms.	`100`
`tol`	`float`	Series truncation: stop as soon as $\|K\|^n < \text{tol}$.	`1e-06`

Returns:

Type	Description
`FieldResult`	Result object. :attr:`FieldResult.metadata` exposes diagnostic fields (`backend`, `K`, `n_terms_used`, `converged`).

Source code in src/groundfield/solver/image_2layer.py

def solve_image_2layer(
    world: "World",
    engine: "Engine",
    *,
    max_terms: int = 100,
    tol: float = 1e-6,
) -> FieldResult:
    """Image-charge solver for 2-layer soil (Tagg / Sunde).

    Parameters
    ----------
    world
        World whose ``soil`` is a :class:`TwoLayerSoil`. All electrodes
        must lie inside the upper layer.
    engine
        Engine configuration. ``engine.segment_length`` controls the
        discretisation as in the homogeneous backend.
    max_terms
        Maximum number of series terms.
    tol
        Series truncation: stop as soon as $|K|^n < \\text{tol}$.

    Returns
    -------
    FieldResult
        Result object. :attr:`FieldResult.metadata` exposes diagnostic
        fields (``backend``, ``K``, ``n_terms_used``, ``converged``).
    """
    if not isinstance(world.soil, TwoLayerSoil):
        raise TypeError(
            "Backend 'image_2layer' requires TwoLayerSoil. "
            f"Got: {type(world.soil).__name__}."
        )
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")

    soil: TwoLayerSoil = world.soil
    h_1 = soil.h_1
    K = soil.reflection_coefficient
    ds = engine.segment_length

    _log.info(
        "image_2layer: rho_1=%.1f, rho_2=%.1f, h_1=%.2f, K=%.4f",
        soil.rho_1, soil.rho_2, h_1, K,
    )

    # 1) Discretisation — identical to the homogeneous backend, but
    #    with segments split at the layer interface (ADR-0007 step 1,
    #    audit 2026-07-08 WP-D2) so no segment straddles h_1.
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds, layer_interfaces=(h_1,))
        elec_to_segidx[e.name] = list(range(len(all_segments),
                                            len(all_segments) + len(segs)))
        all_segments.extend(segs)

    # 2) Per-electrode input currents from sources
    elec_input_current: dict[str, complex] = {
        e.name: 0j for e in world.electrodes
    }
    n_ignored_sources = 0
    for src in world.sources:
        if src.kind != "current":
            n_ignored_sources += 1
            continue
        i_complex = src.magnitude * np.exp(1j * np.deg2rad(src.phase_deg))
        if src.attached_to in elec_input_current:
            elec_input_current[src.attached_to] += i_complex
    if n_ignored_sources:
        warnings.warn(
            f"image_2layer: {n_ignored_sources} non-current source(s) "
            "ignored — only CurrentSource is supported.",
            UserWarning,
            stacklevel=2,
        )

    # 3) Cluster building (ideal conductors only) and finite-impedance
    #    branch list (passed into the nodal-analysis solver).
    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)

    # 3b) Distributed-conductor topology (ADR-0003) + ADR-0004
    #     inductive coupling assembly.
    cond_segs, distributed_branches_objs, interior_nodes = _build_distributed_topology(
        world.conductors, cluster_id
    )
    pseudo_owners: list[str] = []
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
        pseudo_owners.append(pn)
    for n_ in interior_nodes:
        if n_ not in cluster_id:
            cluster_id[n_] = n_
            pseudo_owners.append(n_)
            elec_to_segidx[n_] = []
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [
        (db.node_a, db.node_b, db.R) for db in distributed_branches_objs
    ]
    finite_branches = list(finite_branches) + distributed_branch_tuples
    earth_inductive_model = getattr(
        engine, "earth_inductive_model", "perfect_mirror"
    )
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity

        sigma_earth_for_carson = resolve_earth_conductivity(world.soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers

        layered_earth_for_sommerfeld = resolve_earth_layers(world.soil)
    inductance_matrix_full, has_inductance, carson_builder = _assemble_inductance_matrix(
        distributed_branches_objs,
        n_lumped_branches=n_lumped_branches,
        n_total_branches=len(finite_branches),
        earth_model=earth_inductive_model,
        sigma_earth=sigma_earth_for_carson,
        layered_earth=layered_earth_for_sommerfeld,
    )

    seg_points = np.array([s.midpoint for s in all_segments])
    seg_lengths = np.array([s.length for s in all_segments])
    wire_radii = np.array([s.wire_radius for s in all_segments])

    # Burial-depth validation (audit 2026-07-28, F26/F29). The air
    # mirror of the Tagg/Sunde series is the same perfect mirror as in
    # the homogeneous backend, so the same finite-image-separation
    # requirement 2 z >> a applies to every leakage segment.
    _check_segment_depths(all_segments, where="solve_image_2layer")
    # ADR-0012 V2: per-segment radial shell coefficient. Zero for
    # any segment that does not belong to a concrete-encased
    # foundation electrode — historic case preserved bit-exact.
    seg_shell_coeffs = np.array(
        [s.concrete_shell_coefficient_ohm_m for s in all_segments],
        dtype=float,
    )

    # ADR-0007: detect cross-layer geometry. When all segments live
    # inside the upper layer, the historic Tagg/Sunde image series
    # is used (bit-exact regression). When any segment is at or
    # below the interface, the kernel automatically dispatches to
    # the rigorous Sommerfeld path.
    z_max = seg_points[:, 2].max()
    cross_layer = bool(z_max >= h_1)
    if cross_layer:
        _log.info(
            "image_2layer: cross-layer geometry detected "
            "(z_max = %.3f m, h_1 = %.3f m). Using ADR-0007 "
            "Sommerfeld fallback for all source/observer pairs.",
            z_max, h_1,
        )
    elif h_1 - z_max < 2.0 * ds:
        # Audit 2026-07-08 WP-D2: within ~one segment of the
        # interface the point-source treatment of the K-weighted
        # nearest image (distance 2(h_1 − z)) degrades. Surface the
        # margin so the user can refine ds or accept the bias.
        warnings.warn(
            f"image_2layer: deepest segment ends {h_1 - z_max:.2f} m "
            f"above the layer interface (h_1 = {h_1:.2f} m), which is "
            f"less than 2×segment_length ({2.0 * ds:.2f} m). The "
            "image-series self-terms near the interface carry an "
            "elevated discretisation error there — consider a finer "
            "segment_length.",
            UserWarning,
            stacklevel=2,
        )
    # Tag every segment with its layer index (0 = upper, 1 = lower).
    for s in all_segments:
        s.layer_index = 0 if s.midpoint[2] < h_1 else 1

    # 4) Self-kernel closure and frequency loop
    self_kernel = _two_layer_self_kernel_factory(
        soil, max_terms, tol, allow_cross_layer=cross_layer,
    )
    n_segments = len(all_segments)
    n_freq = len(engine.frequencies)
    omegas = [2.0 * np.pi * float(f) for f in engine.frequencies]
    real_electrode_names = {e.name for e in world.electrodes}

    # ADR-0010 Tier 1 (WP-F): the multi-port grounding matrix Z is
    # frequency-independent — share it across the per-frequency calls.
    _mp_cache: dict = {}

    def _solve_at(omega: float) -> tuple[dict[str, complex], np.ndarray]:
        carson_dz = (
            carson_builder(omega) if (has_inductance and carson_builder is not None)
            else None
        )
        elec_total = _solve_cluster_currents(
            electrodes=world.electrodes,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            seg_points=seg_points,
            seg_lengths=seg_lengths,
            wire_radii=wire_radii,
            elec_to_segidx=elec_to_segidx,
            self_kernel=self_kernel,
            finite_branches=finite_branches,
            pseudo_owners=pseudo_owners,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
            shell_coefficients=seg_shell_coeffs,
            multiport_cache=_mp_cache,
        )
        sc = np.zeros(n_segments, dtype=complex)
        for ename, idxs in elec_to_segidx.items():
            if not idxs:
                continue
            I_total = elec_total.get(ename, 0j)
            if I_total == 0j:
                continue
            L_total = seg_lengths[idxs].sum()
            sc[idxs] = I_total * seg_lengths[idxs] / L_total
        return elec_total, sc

    def _phi_batch(sc_list: list[np.ndarray]) -> list[np.ndarray]:
        """Segment-midpoint potentials for a list of current vectors.

        Stacks all real/imaginary parts into one excitation matrix so
        the kernel's O(N²·M_series) tensors are built exactly once
        for the whole frequency set (ADR-0010 Tier 1).
        """
        k = len(sc_list)
        stacked = np.zeros((n_segments, 2 * k))
        for m, sc in enumerate(sc_list):
            stacked[:, m] = sc.real
            stacked[:, k + m] = sc.imag
        if not stacked.any():
            return [np.zeros(n_segments, dtype=complex)] * k
        phi = self_kernel(seg_points, seg_lengths, wire_radii, stacked)
        # ADR-0012 V2: post-solve potential evaluation must carry the
        # same shell augmentation as the inverse-problem assembly
        # inside _solve_cluster_currents. Without this step the
        # electrode_potentials and cluster_impedance ignore the shell
        # drop and report the bare-soil values, hiding the V2 effect.
        if np.any(seg_shell_coeffs > 0.0):
            with np.errstate(divide="ignore", invalid="ignore"):
                shell_diag = np.where(
                    seg_lengths > 0.0,
                    seg_shell_coeffs / seg_lengths,
                    0.0,
                )
            phi = phi + shell_diag[:, None] * stacked
        return [phi[:, m] + 1j * phi[:, k + m] for m in range(k)]

    elec_per_freq: list[dict[str, complex]] = []
    sc_per_freq: list[np.ndarray] = []
    phi_per_freq: list[np.ndarray] = []
    if has_inductance:
        for omega in omegas:
            et, sc = _solve_at(omega)
            elec_per_freq.append(et)
            sc_per_freq.append(sc)
        phi_per_freq = _phi_batch(sc_per_freq)
    else:
        et, sc = _solve_at(0.0)
        ph = _phi_batch([sc])[0]
        elec_per_freq = [et] * n_freq
        sc_per_freq = [sc] * n_freq
        phi_per_freq = [ph] * n_freq

    n_terms_used_self = 0
    if any(np.any(sc) for sc in sc_per_freq):
        _, _, _, n_terms_used_self = _two_layer_image_offsets(
            K, h_1, max_terms, tol
        )

    electrode_potentials: dict[str, list[complex]] = {}
    electrode_currents: dict[str, list[complex]] = {}
    conductor_currents: dict[str, list[complex]] = {}
    conductor_potentials: dict[str, list[complex]] = {}
    for ename, idxs in elec_to_segidx.items():
        if not idxs:
            continue
        # Length-weighted (Galerkin) node potential — the same
        # reduction the multi-port matrix enforces inside
        # _solve_cluster_currents (audit 2026-07-28, F10/F11). The
        # rod-split-at-the-interface case of ADR-0007 is exactly where
        # the historic unweighted mean diverged from the solved value.
        u_list = _weighted_node_potential(
            phi_per_freq, idxs, seg_lengths, n_freq
        )
        i_list = [elec_per_freq[k][ename] for k in range(n_freq)]
        if ename in real_electrode_names:
            electrode_potentials[ename] = u_list
            electrode_currents[ename] = i_list
        else:
            conductor_potentials[ename] = u_list
            conductor_currents[ename] = i_list

    # 7) Point-source list for post-processing
    point_sources = [
        PointSource(
            position=tuple(seg_points[i].tolist()),
            current=[complex(sc_per_freq[k][i]) for k in range(n_freq)],
            electrode_name=all_segments[i].electrode_name,
            length=float(seg_lengths[i]),
        )
        for i in range(n_segments)
    ]

    cluster_members: dict[str, list[str]] = {}
    for ename in real_electrode_names:
        cluster_members[ename] = sorted(
            n for n in cluster_id
            if cluster_id[n] == cluster_id[ename] and n in real_electrode_names
        )

    converged = (
        _series_tail_bound(abs(K), n_terms_used_self) < tol
        if n_terms_used_self else True
    )
    if not converged:
        tail = _series_tail_bound(abs(K), n_terms_used_self)
        _log.warning(
            "image_2layer: max_terms=%d reached, series tail bound "
            "|K|^(n+1)/(1-|K|) = %.2e > tol=%.2e. Result may be "
            "inaccurate.",
            max_terms, tail, tol,
        )
        warnings.warn(
            f"image_2layer: the Tagg/Sunde image series was truncated "
            f"at max_terms={max_terms} with a remaining tail bound of "
            f"{tail:.2e} > tol={tol:.2e} (|K| = {abs(K):.3f}). Raise "
            "Engine.image_max_terms for high layer contrasts.",
            SeriesTruncationWarning,
            stacklevel=2,
        )

    metadata = {
        "world_name": world.name,
        "n_segments": n_segments,
        "segment_length": ds,
        "K": float(K),
        "rho_1": float(soil.rho_1),
        "rho_2": float(soil.rho_2),
        "h_1": float(h_1),
        "n_terms_used": n_terms_used_self,
        "converged": converged,
        "stub": False,
        "earth_inductive_model": earth_inductive_model,
    }
    if has_inductance and sigma_earth_for_carson is not None:
        from groundfield.coupling.carson import skin_depth

        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_earth_for_carson)
            for f in engine.frequencies
        }
    elif has_inductance:
        from groundfield.coupling.carson import skin_depth

        sigma_ref = 1.0 / float(soil.rho_1)
        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_ref)
            for f in engine.frequencies
        }
    if conductor_currents:
        metadata["conductor_node_currents"] = conductor_currents
        metadata["conductor_node_potentials"] = conductor_potentials

    return FieldResult(
        backend="image_2layer",
        frequencies=list(engine.frequencies),
        electrode_potentials=electrode_potentials,
        electrode_currents=electrode_currents,
        point_sources=point_sources,
        soil_resistivity=float(soil.rho_1),
        soil=soil,
        clusters=cluster_members,
        metadata=metadata,
    )

ADR-0001 — original methodology decision.

Counterpart	Expected agreement	What is checked
`image` at \(K_1 = 0\)	bit-exact	series collapse
`mom` on the same world	\(\le 2\,\%\)	same kernel, different resolution scheme
`cim` on the same world	bit-exact (CIM falls back on this kernel for \(n = 2\))	`cim` deliberately reuses the Tagg / Sunde self-kernel for \(n = 2\)
`bem` on the same world	bit-exact (BEM falls back on this kernel for \(n = 2\))	same reason
`mom_sommerfeld` on the same world	\(\le 5\,\%\)	quadrature reproduces the geometric series
Sunde / Dwight closed forms	\(\le 10\,\%\)	reference values for canonical geometries (rods, rings)

image_2layer — Tagg / Sunde series for 2-layer soil¶