3D forward modelling

MTGeophysics.jl has no 3D solver of its own. For 3D it uses ModEM, a well-tested MPI-parallel solver, and handles everything around it: reading and writing the files, running ModEM, and comparing its predictions with the data. In this section we look at each of these steps.

Read and write 3D files

using MTGeophysics

d = load_data_modem("data.dat")          # impedance tensor and tipper
m = load_ws3d_model("model.rho")         # WS3D model, log10 ρ in memory
write_ws3d_model("new.rho", m.dx, m.dy, m.dz, m.A, m.origin)

WS3D is the model format of ModEM 3D. load_ws3d_model marks air cells as NaN, so topography is easy to recognise in memory.

Run ModEM

A ModEM forward run takes a model and a data file and writes the predicted data:

mpirun -n 21 Mod3DMT -F model.rho data.dat pred.dat

-F means forward, and -n sets the number of MPI processes. The 3D VFSA inversion runs this command for you at every step.

Compare with the data

fit = chi2_and_rms("data.dat", "pred.dat")
fit.rms

The RMS is computed from the impedances and, when present, the tipper.

Some sites suffer from galvanic distortion: small bodies near the surface that scale the electric field. chi2_and_rms_distorted fits a real 2 × 2 distortion matrix to each site before computing the misfit, which shows how much of the misfit is due to distortion:

fit = chi2_and_rms_distorted("data.dat", "pred.dat"; damping = 0.01)

damping pulls each matrix towards the identity: 0 lets it fit freely, larger values allow less correction.

Apparent resistivity and phase

load_data_modem also computes the apparent resistivity and phase of every component, with their errors, as d.ρ, d.φ, d.ρerr and d.φerr. Their size is (periods, components, sites).