pycsamt.forward.em3d#

Quasi-3D magnetotelluric forward solver.

MT3DForward approximates the full 3-D MT impedance tensor by running MT2DForward on pairs of orthogonal 2-D cross-sections (XZ and YZ planes) extracted from a Grid3D model.

Physics of the quasi-3D approximation#

For a 3-D resistivity model σ(x, y, z), the full impedance tensor is:

Z = [ Z_xx Z_xy ]

[ Z_yx Z_yy ]

The quasi-3D method splits the tensor estimation into two independent 2-D FD solves:

  1. XZ profiles — one per unique y-row of stations. Each solve treats the model as 2-D with structure in the xz-plane. From the 2-D TE mode (E_y / H_x) we get Z_xy; from TM (E_x / H_y) we get Z_yx.

  2. YZ profiles — one per unique x-column of stations. The solver uses the yz-slice resistivity with y as the horizontal axis. TE in this plane gives Z_xy from the NS structure; TM gives Z_yx from the NS structure.

The final quasi-3D tensor is the arithmetic average of contributions from both profile directions:

Z_xy ← ½ (Z_xy_xz + Z_xy_yz) Z_yx ← ½ (Z_yx_xz + Z_yx_yz) Z_xx = Z_yy = 0 (valid for structures without azimuthal asymmetry)

For a 1-D earth both profiles give the same response and the average reduces to the exact 1-D solution. For 2-D structures (variation in one horizontal direction only), one profile captures the full signal and the other gives the background, so the average is intermediate. For genuine 3-D structures, both profiles contribute and the result is an approximation that correctly captures the leading-order lateral variation — sufficient for generating AI training data.

Accuracy and limitations#

  • Off-diagonal components (Z_xx, Z_yy) are set to zero; the true 3-D Z_xx and Z_yy are typically 5–20 % of Z_xy in magnitude but are important for strike estimation.

  • The approximation breaks down for strongly 3-D structures where galvanic coupling between adjacent columns is significant (e.g. thin resistive dykes).

  • For the AI training use case, the quasi-3D data vastly improves on the fully independent-1D pseudo-3D baseline by injecting realistic lateral coupling physics.

References

Wannamaker, P.E. (1999). Affordable magnetotellurics: Interpretation in natural environments. SEG Distinguished Instructors Short Course, 7.

Mackie, R.L. et al. (1993). Three-dimensional electromagnetic modeling using finite differences. Geophysics, 58(2), 215-226.

Classes

ForwardResponse3D(freqs, stations_xy, zxy, ...)

Full (approximate) impedance tensor from the 3-D MT forward solver.

MT3DForward(freqs, grid, *[, method, verbose])

Quasi-3D magnetotelluric forward solver.
class pycsamt.forward.em3d.MT3DForward(freqs, grid, *, method='quasi3d', verbose=True)[source]#

Bases: object

Quasi-3D magnetotelluric forward solver.

Runs MT2DForward on orthogonal XZ and YZ cross-sections extracted from a Grid3D model and assembles the result into an approximate full impedance tensor.

Parameters:

  • freqs (array-like) – Frequencies [Hz].

  • grid (Grid3D) – 3-D resistivity model with surface station positions.

  • method ({'quasi3d'}) – Forward solver method. Currently only 'quasi3d' is implemented. A true 3-D FD Yee-grid solver ('fd3d') is planned for a future phase.

  • verbose (bool) – Print per-profile progress.

Examples

Halfspace — should recover the 1-D MT response:

>>> import numpy as np
>>> from pycsamt.forward.grid3d import Grid3D
>>> from pycsamt.forward.em3d import MT3DForward

>>> freqs = np.logspace(-1, 2, 8)
>>> g = Grid3D.halfspace(rho=100.0, nx=16, ny=16, nz=12,
...     x_max=5_000.0, y_max=5_000.0, z_max=3_000.0,
...     nx_stations=3, ny_stations=3)
>>> resp = MT3DForward(freqs, g).run()
>>> resp.rho_a_xy.shape
(8, 9)

3-D conductive block:

>>> g = Grid3D.block_anomaly(bg_rho=500.0, anomaly_rho=5.0, nx=20, ny=20, nz=15,
...     x_max=6_000.0, y_max=6_000.0, z_max=4_000.0,
...     nx_stations=4, ny_stations=4)
>>> resp = MT3DForward(freqs, g, verbose=False).run()
>>> resp.rho_a_xy.std(axis=1).max() > 1.0   # lateral variation present
True
run()[source]#

Run the forward solver and return the full impedance tensor.

Return type:

ForwardResponse3D

class pycsamt.forward.em3d.ForwardResponse3D(freqs, stations_xy, zxy, zyx, zxx, zyy, rho_a_xy, phase_xy, rho_a_yx, phase_yx, rho_a_xx, phase_xx, rho_a_yy, phase_yy, method='quasi3d', grid=None)[source]#

Bases: object

Full (approximate) impedance tensor from the 3-D MT forward solver.

All arrays have shape (n_freqs, n_stations).

Parameters:

  • freqs (ndarray, shape (n_freqs,))

  • stations_xy (ndarray, shape (n_stations, 2)) – Station (x, y) positions [m].

  • zxy (ndarray of complex) – Impedance tensor components [V/A].

  • zyx (ndarray of complex) – Impedance tensor components [V/A].

  • zxx (ndarray of complex) – Impedance tensor components [V/A].

  • zyy (ndarray of complex) – Impedance tensor components [V/A].

  • rho_a_xy (ndarray) – Apparent resistivities [Ω·m].

  • rho_a_yx (ndarray) – Apparent resistivities [Ω·m].

  • rho_a_xx (ndarray) – Apparent resistivities [Ω·m].

  • rho_a_yy (ndarray) – Apparent resistivities [Ω·m].

  • phase_xy (ndarray) – Impedance phases [degrees].

  • phase_yx (ndarray) – Impedance phases [degrees].

  • phase_xx (ndarray) – Impedance phases [degrees].

  • phase_yy (ndarray) – Impedance phases [degrees].

  • method (str) – Solver method ('quasi3d').

  • grid (Grid3D) – Source model.

freqs: ndarray#
stations_xy: ndarray#
zxy: ndarray#
zyx: ndarray#
zxx: ndarray#
zyy: ndarray#
rho_a_xy: ndarray#
phase_xy: ndarray#
rho_a_yx: ndarray#
phase_yx: ndarray#
rho_a_xx: ndarray#
phase_xx: ndarray#
rho_a_yy: ndarray#
phase_yy: ndarray#
method: str = 'quasi3d'#
grid: Grid3D = None#
property n_freqs: int[source]#
property n_stations: int[source]#
property periods: ndarray[source]#
to_feature_array(*, components='xy_yx', log_rho=True, include_phase=True)[source]#

Flatten to a 2-D feature matrix for ML input.

Parameters:

  • components (str) – Comma/underscore-separated component names. E.g. "xy_yx" (default) uses Z_xy and Z_yx; "xy" uses only Z_xy; "all" uses all four.

  • log_rho (bool) – Return log₁₀(ρ_a) instead of ρ_a.

  • include_phase (bool) – Concatenate phase alongside ρ_a.

Returns:

X

Return type:

ndarray, shape (n_stations, n_features)

to_survey_dataset(y_models=None, *, components='xy_yx')[source]#

Wrap this response as a single-survey SurveyDataset3D.

The result can be concatenated with other surveys to build a training dataset for GCNInverter3D.

Parameters:

  • y_models (ndarray or None) – Per-station model target vectors, shape (n_stations, n_params). When None a zero placeholder is used.

  • components (str) – Passed to to_feature_array().

Returns:

Single-survey dataset (n_surveys = 1).

Return type:

SurveyDataset3D

station_response(station)[source]#

Return all tensor components for one station as a dict.

Parameters:

station (int)

Return type:

dict