2.20.3.5. pycsamt.forward.em2d#

2-D magnetotelluric finite-difference forward solver.

MT2DForward solves the 2-D MT boundary-value problem for both the TE and TM polarisation modes on a non-uniform staggered grid using the finite-difference (FD) method.

Physics#

For a 2-D earth with conductivity σ(x, z) and a plane electromagnetic wave incident from above, the governing scalar PDEs (time convention e^{+iωt}) are:

TE mode (electric field parallel to strike, E_y):

∂²E_y/∂x² + ∂²E_y/∂z² − iωμ₀σ(x,z) E_y = 0

TM mode (magnetic field parallel to strike, H_y):

∂/∂x[σ⁻¹ ∂H_y/∂x] + ∂/∂z[σ⁻¹ ∂H_y/∂z] − iωμ₀ H_y = 0

Grid and node numbering#

The FD grid has (nx+1) × (nz+1) nodes. Node at column j, row i (depth) has global index k = i*(nx+1) + j. The surface is row i = 0, the bottom is row i = nz.

Boundary conditions#

Dirichlet conditions are applied at all four boundaries using the analytic 1-D MT response of the edge column profiles (left/right) or the halfspace (top/bottom). For the TE mode the normalised plane-wave source E_y = 1 at the surface is the natural starting point; the BC values at depth follow from the 1-D recursion. For the TM mode the analogous 1-D H_y field is used.

Each frequency requires two independent sparse solves (one per mode). With scipy.sparse.linalg.spsolve and grids of order 50 × 40 nodes (~2000 unknowns) each solve takes < 0.1 s on a single CPU core.

References

Wannamaker, P.E., Stodt, J.A., Rijo, L. (1987). A stable finite element solution for two-dimensional magnetotelluric modelling. Geophys. J. Int. 88, 277-296.

de Groot-Hedlin, C., Constable, S. (1990). Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data. Geophysics 55, 1613-1624.

Wait, J.R. (1954). On the relation between telluric currents and the Earth’s magnetic field. Geophysics 19, 281-289.

Classes

ForwardResponse2D(freqs, stations_x, zxy, ...)

Output of the 2-D MT forward solver.

MT2DForward(freqs, grid, *[, verbose])

2-D magnetotelluric finite-difference forward solver.
class pycsamt.forward.em2d.MT2DForward(freqs, grid, *, verbose=True)[source]

Bases: object

2-D magnetotelluric finite-difference forward solver.

Solves both TE and TM modes for a list of frequencies on the provided Grid2D and returns a ForwardResponse2D with apparent resistivity and phase at all station positions.

Parameters:

  • freqs (array-like) – Frequencies [Hz] at which to evaluate the response.

  • grid (Grid2D) – 2-D resistivity model and station layout.

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

Examples

Uniform halfspace — should recover the 1-D MT response:

>>> import numpy as np
>>> from pycsamt.forward.grid2d import Grid2D
>>> from pycsamt.forward.em2d import MT2DForward

>>> freqs = np.logspace(-2, 3, 10)
>>> g = Grid2D.halfspace(rho=100.0, nx=30, nz=20,
...                       x_max=5_000.0, z_max=3_000.0,
...                       n_stations=5)
>>> fwd = MT2DForward(freqs, g)
>>> resp = fwd.run()
>>> resp.rho_a_te.shape
(10, 5)

2-D model with conductive anomaly:

>>> g2 = Grid2D.with_anomaly(bg_rho=100.0, anomaly_rho=2.0,
...     anomaly_bounds=(1500.0, 3500.0, 300.0, 900.0),
...     nx=40, nz=25, x_max=6000.0, z_max=3000.0, n_stations=10)
>>> resp2 = MT2DForward(freqs, g2).run()
run()[source]

Run the forward solver for all frequencies.

Return type:

ForwardResponse2D

class pycsamt.forward.em2d.ForwardResponse2D(freqs, stations_x, zxy, zyx, rho_a_te, phase_te, rho_a_tm, phase_tm, residual_te, residual_tm, grid)[source]

Bases: object

Output of the 2-D MT forward solver.

All impedance and apparent-resistivity arrays have shape (n_freqs, n_stations).

Parameters:

  • freqs (ndarray, shape (n_freqs,)) – Evaluation frequencies [Hz].

  • stations_x (ndarray, shape (n_stations,)) – Station x-positions [m].

  • zxy (ndarray, shape (n_freqs, n_stations), complex) – TE-mode surface impedance Z_xy = E_y / H_x [V/A].

  • zyx (ndarray, shape (n_freqs, n_stations), complex) – TM-mode surface impedance Z_yx = E_x / H_y [V/A] (negative by convention: Z_yx = −Z_xy for a 1-D earth).

  • rho_a_te (ndarray, shape (n_freqs, n_stations)) – TE apparent resistivity [Ω·m].

  • phase_te (ndarray, shape (n_freqs, n_stations)) – TE impedance phase [degrees, 0–90° for a normal model].

  • rho_a_tm (ndarray, shape (n_freqs, n_stations)) – TM apparent resistivity [Ω·m].

  • phase_tm (ndarray, shape (n_freqs, n_stations)) – TM impedance phase [degrees].

  • residual_te (ndarray, shape (n_freqs,)) – Relative linear-system residual norm(A @ x - b) / norm(b) of the direct sparse solve for each mode, one value per frequency. This measures the achieved numerical accuracy of scipy.sparse.linalg.spsolve; it is not an iteration count, since the solver is direct.

  • residual_tm (ndarray, shape (n_freqs,)) – Relative linear-system residual norm(A @ x - b) / norm(b) of the direct sparse solve for each mode, one value per frequency. This measures the achieved numerical accuracy of scipy.sparse.linalg.spsolve; it is not an iteration count, since the solver is direct.

  • grid (Grid2D) – The model grid used for the forward run.

freqs: ndarray
stations_x: ndarray
zxy: ndarray
zyx: ndarray
rho_a_te: ndarray
phase_te: ndarray
rho_a_tm: ndarray
phase_tm: ndarray
residual_te: ndarray
residual_tm: ndarray
grid: Grid2D
property n_freqs: int[source]
property n_stations: int[source]
property periods: ndarray[source]

Periods [s].

station_response(station)[source]

Return all response arrays for one station as a dict.

Parameters:

station (int)

Return type:

dict

to_feature_array(*, mode='both', log_rho=True, include_phase=True)[source]

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

Parameters:

  • mode ({'te', 'tm', 'both'}) – Which mode(s) to include.

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

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

Return type:

ndarray, shape (n_stations, n_features)

plot(station=0, ax=None, *, figsize=(9, 5))[source]

Quick diagnostic plot of ρ_a and phase for one station.

Parameters:

  • station (int) – Station index.

  • ax (array of Axes or None) – Two axes (rho, phase). Created when not provided.

  • figsize (tuple)

Returns:

axes

Return type:

ndarray of Axes, shape (2,)