Skip to content

image_nlayer — image-charge dispatcher for n-layer soil

Physical context

The image-charge family of solvers represents the layered Green's function as a sum of point sources at mirrored positions in the soil. This representation is closed-form for two regimes:

  • n = 1 (homogeneous). A single source plus its air-mirror at \(z = -z_s\) reproduces the Neumann boundary at \(z = 0\) exactly.
  • n = 2 (two-layer). \(\Gamma_1(\lambda) \equiv K_1\) is constant in \(\lambda\), so the multiple-reflection expansion of \(1 / (1 - K_1 e^{-2\lambda h_1})\) is a geometric series whose coefficients are the well-known \(K_1^n\) Tagg / Sunde weights.

For \(n \ge 3\) the recursive reflection coefficient

\[ \Gamma_1(\lambda) \;=\; \frac{K_1 + \Gamma_2(\lambda)\, e^{-2\lambda h_2}} {1 + K_1\, \Gamma_2(\lambda)\, e^{-2\lambda h_2}} \]

becomes a non-trivial rational function of \(e^{-2\lambda h_i}\) that no longer collapses to a clean geometric series in the spatial domain. Stefanescu / Sunde (1968, ch. 3.5) describe a doubly-nested expansion that recovers a real-image series, but its implementation is fragile for hard contrasts: the convergence factor at every interface is \(|K_i|\), and the doubly-nested loop has to converge in every nested level. For practical typical contrasts this works only with non-trivial Aitken / Pade acceleration on top.

image_nlayer is the dispatcher that ties the closed-form image-charge family together. Rather than carrying a fragile \(n \ge 3\) implementation, it delegates to the engines that are designed for that regime (cim, bem, mom_sommerfeld) and raises a clear error otherwise. ADR-0002 documents this design decision.

Dispatch table

n_layers Effective backend Notes
1 image Homogeneous half-space; the soil is internally cast to HomogeneousSoil if a MultiLayerSoil with one layer is supplied.
2 image_2layer Tagg / Sunde geometric series. A MultiLayerSoil with exactly two entries is internally cast to TwoLayerSoil.
≥ 3 — (raises ValueError) The error message points to cim, mom_sommerfeld, or bem for the regime.

The returned FieldResult.backend is rewritten to "image_nlayer" so that cross-engine comparisons see one unified label per backend selection. The original delegate is preserved in FieldResult.metadata['dispatched_to'].

Why no \(n \ge 3\) implementation?

Three reasons, each independently sufficient:

  1. Convergence sensitivity. The doubly-nested Stefanescu series converges as \(\prod_i |K_i \Gamma_{i+1}|^{n_i}\). For typical parameters with \(|K_i| \approx 0.5\)–\(0.9\), the series can take \(\gtrsim 10^4\) terms in the worst-case nested combinations. Without Aitken / Pade acceleration the engine would routinely fail to converge inside max_terms.

  2. Methodological redundancy. The complex-image method (cim) solves the same problem with \(\sim 8\) complex exponentials and gives a closed-form spatial Green's function for any \(n\). The per-evaluation cost is independent of the layer count, and the matrix-pencil fit is numerically robust for the typical contrast range.

  3. Cross-validation strength. mom_sommerfeld provides the absolute reference for \(n \ge 3\) (direct numerical Sommerfeld quadrature), and bem provides an alternative weighting (CIM kernel, collocation). Adding a fragile real-image series would not increase confidence — it would just add another engine that the cross-validation matrix would have to discount when it diverges.

The trade-off is documented in ADR-0002 and revisited every time the supported use-case scope expands.

Backend ergonomics

Engine.solve automatically forwards backend="image" to:

  • image_2layer if the world holds a TwoLayerSoil,
  • image_nlayer if the world holds a MultiLayerSoil,

so client code written for the homogeneous case keeps working when only the soil model is swapped. The path through image_nlayer for an \(n \ge 3\) soil therefore terminates with a clear ValueError suggesting the appropriate alternative engine — useful safety rail for users who copy-paste a 2-layer setup and only later upgrade the soil model.

Validity envelope

Property Range / value
Soil model HomogeneousSoil, TwoLayerSoil, or MultiLayerSoil
Effective layer count 1 or 2 only — n ≥ 3 raises ValueError
Frequency inherited from delegate (quasi-static, \(f < 1\,\text{kHz}\))
Electrode placement every segment must lie inside the upper layer

Cross-validation notes

By construction the dispatcher reproduces its delegate bit-exactly:

Stack Delegate Expected agreement
HomogeneousSoil image \(10^{-9}\) (identity)
TwoLayerSoil image_2layer \(10^{-9}\) (identity)
MultiLayerSoil, \(n = 1\) image (auto-cast) \(10^{-9}\)
MultiLayerSoil, \(n = 2\) image_2layer (auto-cast) \(10^{-9}\)
MultiLayerSoil, \(n \ge 3\) — raises

The dispatcher therefore appears in the cross-engine test matrix as a "free" extra engine for the \(n = 1, 2\) cases, with a tolerance of \(10^{-9}\) against its delegate.

References

The dispatcher itself does not introduce new mathematics. The delegate engines carry their own references on image and image_2layer. For the \(n \ge 3\) design decision see:

  • Sunde, E. D. (1968). Earth Conduction Effects in Transmission Systems, Dover, ch. 3.5. Stefanescu series for \(n\) layers and why it converges slowly for high contrasts.
  • 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 reformulations that would be required to make a real n-layer image series practical.
  • ADR-0002 — the in-repo justification of the dispatcher design.

Example

import groundfield as gf

# 2-layer soil — dispatcher delegates to image_2layer.
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, "rod", name="g1",
                    position=(0.0, 0.0, 0.0), length=1.5)
gf.create_source(world, attached_to="g1", magnitude=1.0)

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

API reference

image_nlayer

Image-charge backend dispatcher for n-layer soil (image_nlayer).

Mathematical / physical model

The image-charge family of grounding solvers represents the layered half-space Green's function as a (truncated) sum of point sources at mirrored positions in the soil. This sum is closed-form for two regimes:

  • n = 1 (homogeneous soil). A single source plus its air-mirror at \(z = -z_s\) reproduces the Neumann boundary at \(z = 0\) exactly. This is the classical image backend.

  • n = 2 (two-layer soil). The recursive reflection \(\Gamma_1(\lambda) \equiv K_1\) is constant in \(\lambda\), and the geometric expansion of the multiple reflection between the air boundary and the layer interface gives the Tagg / Sunde series $$ \varphi(s, z) \;=\; \frac{\rho_1\, I}{4\pi}\, \sum_{n=0}^{\infty} K_1^{n}\, \Bigl(\tfrac{1}{r_n^{++}} + \tfrac{1}{r_n^{+-}}

    • \tfrac{1}{r_n^{-+}} + \tfrac{1}{r_n^{--}}\Bigr). $$ For \(n \ge 3\) the upward reflection \(\Gamma_1(\lambda)\) becomes a non-trivial function of \(\lambda\) (see :func:groundfield.solver._layered.reflection_gamma). Expanding it as a closed-form real image-charge series exists in principle (Stefanescu / Sunde 1968, ch. 3.5) but the implementation is fragile for hard contrasts because of the doubly-nested geometric expansion. Within the groundfield engine family this regime is therefore covered by

  • cim — Complex Image Method (closed form, complex images).

  • mom_sommerfeld — direct numerical Sommerfeld quadrature (reference engine).
  • bem — boundary-element collocation with CIM kernel.
Dispatch behaviour

solve_image_nlayer is a small wrapper that recognises the layer count and forwards the call to the matching closed-form backend:

========= =============== =========================== n_layers actual backend note ========= =============== =========================== 1 image homogeneous case 2 image_2layer Tagg / Sunde geometric series ≥ 3 — (raises) use cim / mom_sommerfeld / bem ========= =============== ===========================

The backend tag of the returned :class:FieldResult is rewritten to "image_nlayer" so that comparisons across engines see a single identifier per backend selection.

Validity
  • Quasi-static, \(f < 1\,\mathrm{kHz}\).
  • All electrodes must lie inside the upper layer.
References
  • Sunde 1968, ch. 3.5; Tagg 1964, ch. 5.
  • Stefanescu, S. & Schlumberger, C. (1930).
  • Dawalibi, F. P. & Barbeito, N. (1991).

solve_image_nlayer 

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

Image-charge dispatcher for n-layer soil.

Routes to image for n = 1 and to image_2layer for n = 2. For n \ge 3 raises a clear :class:ValueError pointing the caller to the engines that actually support that regime.

Parameters:

Name Type Description Default
world 'World'

World whose soil is any layered soil model.

 required
engine 'Engine'

Engine configuration; engine.segment_length controls the discretisation.

 required
max_terms int

Forwarded to :func:solve_image_2layer for n = 2; ignored otherwise.

 200
tol int

Forwarded to :func:solve_image_2layer for n = 2; ignored otherwise.

 200

Returns:

Type Description
FieldResult

backend is rewritten to "image_nlayer" so that cross-engine reports use one unified label.
Source code in src/groundfield/solver/image_nlayer.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def solve_image_nlayer(
    world: "World",
    engine: "Engine",
    *,
    max_terms: int = 200,
    tol: float = 1e-6,
) -> FieldResult:
    """Image-charge dispatcher for n-layer soil.

    Routes to ``image`` for ``n = 1`` and to ``image_2layer`` for
    ``n = 2``. For ``n \\ge 3`` raises a clear :class:`ValueError`
    pointing the caller to the engines that actually support that
    regime.

    Parameters
    ----------
    world
        World whose soil is any layered soil model.
    engine
        Engine configuration; ``engine.segment_length`` controls the
        discretisation.
    max_terms, tol
        Forwarded to :func:`solve_image_2layer` for ``n = 2``;
        ignored otherwise.

    Returns
    -------
    FieldResult
        ``backend`` is rewritten to ``"image_nlayer"`` so that
        cross-engine reports use one unified label.
    """
    if not isinstance(world.soil, (HomogeneousSoil, TwoLayerSoil, MultiLayerSoil)):
        raise TypeError(
            "Backend 'image_nlayer' supports HomogeneousSoil, "
            "TwoLayerSoil, and MultiLayerSoil. "
            f"Got: {type(world.soil).__name__}."
        )
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")

    stack = as_layer_stack(world.soil)
    n = stack.n_layers
    _log.info("image_nlayer: dispatching for n_layers=%d", n)

    if n == 1:
        # Forward to the homogeneous backend; HomogeneousSoil is the
        # only soil it accepts. If the user supplied a degenerate
        # MultiLayerSoil with one layer, build an equivalent
        # HomogeneousSoil on the fly.
        if isinstance(world.soil, HomogeneousSoil):
            res = solve_image(world, engine)
        else:
            from groundfield.world import World

            tmp = World(
                name=world.name,
                soil=HomogeneousSoil(resistivity=float(stack.rhos[0])),
                electrodes=list(world.electrodes),
                conductors=list(world.conductors),
                sources=list(world.sources),
                boundary=world.boundary,
            )
            res = solve_image(tmp, engine)
    elif n == 2:
        # Forward to the closed-form Tagg/Sunde series.
        if isinstance(world.soil, TwoLayerSoil):
            res = solve_image_2layer(
                world, engine, max_terms=max_terms, tol=tol
            )
        else:
            # MultiLayerSoil with exactly two layers: cast to TwoLayerSoil.
            from groundfield.world import World

            two = TwoLayerSoil(
                rho_1=float(stack.rhos[0]),
                rho_2=float(stack.rhos[1]),
                h_1=float(stack.h[0]),
            )
            tmp = World(
                name=world.name,
                soil=two,
                electrodes=list(world.electrodes),
                conductors=list(world.conductors),
                sources=list(world.sources),
                boundary=world.boundary,
            )
            res = solve_image_2layer(
                tmp, engine, max_terms=max_terms, tol=tol
            )
    else:
        raise ValueError(
            f"image_nlayer: layer count {n} ≥ 3 is not supported by "
            "the real image-charge series (Γ_1(λ) is no longer "
            "constant in λ). Use 'mom_sommerfeld' (direct Sommerfeld "
            "quadrature of the full layered Green's function) or "
            "'fem' — the 'cim'/'bem' n ≥ 3 kernels were found "
            "structurally incomplete in the 2026-07-08 audit and "
            "are rejected as well."
        )

    # Rewrite backend tag so that cross-engine reports use the unified label.
    res = res.model_copy(
        update={
            "backend": "image_nlayer",
            "metadata": {**res.metadata, "n_layers": int(n), "dispatched_to": res.backend},
        }
    )
    return res
  • ADR-0002 — engine selection heuristic.