`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 and all $n \ge 1$ image terms are at least $2 z_s$ (or $2 h_1 - z - z_s$, etc.) away from the segment midpoint, so the point-source approximation is safe.

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.

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)
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])
    # 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
        u_list = [
            complex(np.mean(phi_per_freq[k][idxs])) for k in range(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¶