Numerical Optimization for EM Inversion — a Step-by-Step Guide (1D → 2D → 3D)

efore we write a single line of code, let's start with a simple question: How do you see something you can't touch?

Let's assume you're on a boat in a thick fog. You can't see the seafloor, but you have a sonar system. You send a "ping" (a sound wave) into the water. This is the "forward problem"—using known physics to predict what will happen. A few seconds later, you hear an echo. That echo is your "data". Your brain instantly does the math: "Based on the echo's timing and strength, the seafloor must be 100 meters deep and rocky." That mental calculation—turning observed data back into an unseeable model of the world—is the essence of inversion.

A doctor uses a CT scanner for the same reason. They can't see inside your body, so they pass X-rays (physics) through you and measure what comes out (data). A powerful computer then inverts that data to build a 3D model of your insides. Numerical Inversion is our "CT scanner" for the Earth. It's the mathematical and computational process of taking our limited, noisy measurements from the surface and using them to build a reliable picture, or "model," of the Earth's invisible subsurface.

Instead of sending a sonar "ping" or X-rays, we use electromagnetism (EM). An EM inversion does this:

This isn't just an academic exercise. EM inversion is the engine behind finding:

In all these cases, drilling is expensive and risky. Inversion is what lets us "drill" a thousand times on a computer first, so we only have to drill once in the real world. Now, let's build the engine.

his is the "sonar ping" for our problem. In geophysical inversion, we first need a "Forward Model"—the physics that predicts what we should measure. For frequency-domain EM, that physics is governed by the vector wave equation.

This equation is our workhorse. It describes how an electric field $\mathbf{E}$ (our "ping") interacts with the subsurface. The key is the term $k^2(\mathbf{x})$, which contains the Earth's properties: $k^2(\mathbf{x}) = \omega^2\mu\epsilon(\mathbf{x}) - i\omega\mu\sigma(\mathbf{x})$. This means if we change the conductivity $\sigma$ (our model), the equation gives us a different predicted response (our data).

We can summarize this whole physical simulation with one simple function: $\mathbf{F}(\mathbf{m})$. You give it a model $\mathbf{m}$ (the conductivity map), and it returns the simulated data $\mathbf{d}_{\text{sim}} = \mathbf{F}(\mathbf{m})$. That's the forward problem. Now for the inverse problem. We already have our real measurements, $\mathbf{d}_{\text{obs}}$. Our goal is to find the one model $\mathbf{m}$ that is:

We enforce this balance by creating a "score," called the Objective Function ($\Phi$), that we try to make as low as possible. It has two parts:

Here, $\Phi_d$ is the Data Misfit, which measures how badly our prediction matches the real data (the error). $\Phi_m$ is the Regularization (or "model cost"), which penalizes models that are too rough or complex. The value $\beta$ is a simple slider that "trades off" the two: do we care more about a perfect data fit or a smoother model?

The entire goal of our inversion is now simple: find the single model $\mathbf{m}$ that makes the $\Phi$ score as low as possible. The next steps will show the exact algorithms (Gauss-Newton and Levenberg-Marquardt) we use to find that minimum.

Okay, let's get specific. Our "Forward Model" $\mathbf{F}(\mathbf{m})$ isn't just an abstract idea; it's a specific piece of physics. For most EM surveys, the workhorse is the diffusive form of Maxwell's equations. This is a common simplification of the full wave equation we saw in Eq. (1). At the lower frequencies we use in geophysics, the "conduction current" ($i\omega\sigma\mathbf{E}$) is dominant, and the "displacement current" ($\omega^2\epsilon\mathbf{E}$) is so tiny we can safely neglect it. This gives us our practical, workhorse PDE:

This equation is the core of our function $\mathbf{F}(\mathbf{m})$. It directly links the property we want to find—the conductivity $\sigma(\mathbf{x})$ (our model $\mathbf{m}$)—to the electric field $\mathbf{E}$ that we can measure. This "diffusive" physics gives us a crucial rule of thumb for understanding our problem: the skin depth. $\displaystyle \delta \approx \sqrt{\tfrac{2}{\mu\,\sigma\,\omega}}$

The skin depth $\delta$ tells us how far our EM "ping" can penetrate the earth. Any geological features much smaller than $\delta$ will be invisible or "smeared out," and features much deeper than $\delta$ won't send an echo back. This simple formula is our first reality check for survey design and for deciding if a simpler 1D or 2D model is good enough, or if we need to use a full 3D simulation.

ur forward model (Eq. 3) predicts the full electric $\mathbf{E}$ and magnetic $\mathbf{H}$ fields everywhere in the ground. But we can't measure that. We can only measure a few components at specific receiver locations on the surface. Our instruments record these surface fields and convert them into robust, interpretable ratios. For MT/AMT/CSAMT, this gives us two key pieces of data:

Different methods (like FDEM or TDEM) might record field components directly. But in all cases, these measurements—$\mathbf{Z}$, $\mathbf{T}$, or raw fields—are all just different "views" of the same underlying physics. They are the observed data $\mathbf{d}_{\text{obs}}$ that our inversion will try to match.

Our forward model $\mathbf{F}(\mathbf{m})$ has to solve that PDE (Eq. 3). The cost and complexity of that solve depend entirely on the "dimensionality"—the assumptions we make about the Earth's geometry.

Crucially, the inversion algorithm we'll build in the next steps is agnostic to this choice. The only part that changes is the "physics engine"—the 1D, 2D, or 3D forward solver we plug into it.

We invert for $\mathbf{m}=\log\sigma$ to enforce positivity and to make steps multiplicative in $\sigma$. The numerical meshes (1D layers, 2D/3D cells) are model grids; survey-specific data (impedances, fields) are predicted at receiver locations via interpolation and functional evaluation. Sensitivities are with respect to $\mathbf{m}$: the Jacobian $\mathbf{J}=\partial\mathbf{d}/\partial\mathbf{m}$ is what couples physics to optimization.

Inversion balances two forces: fit the data you trust and stay geologically plausible. We encode this as a sum of a weighted data misfit and a model regularizer:

$\mathbf W_d$ scales each datum by its uncertainty (and optional error floor). If a standard deviation $\sigma_i$ is known, a common choice is $(\mathbf W_d)_{ii}=1/\sigma_i$, so that a “good” fit corresponds to a normalized $\chi^2$ near 1.

$\mathbf W_m$ promotes structure you expect in the earth. For smooth Tikhonov regularization, $\mathbf W_m$ stacks zero-order and first/second differences, e.g. $\mathbf W_m=\operatorname{diag}(\alpha_s\,\mathbf I,\ \alpha_x\,\mathbf D_x,\ \alpha_y\,\mathbf D_y,\ \alpha_z\,\mathbf D_z)$. You can also switch to total-variation (TV) for sharper contacts, or include a reference model $\mathbf m_{\mathrm{ref}}$ to pull solutions toward known geology.

$\beta$ trades off fit and smoothness. We typically start with a larger $\beta$ (stability first), then reduce it (“cool” it) until the data misfit reaches the noise level—this is the discrepancy principle, aiming for $\Phi_d \approx N_{\text{data}}/2$.

To use gradient-based solvers we need the derivative of $\Phi$ with respect to the model:

Here $\mathbf J=\partial \mathbf F/\partial \mathbf m$ is the sensitivity (Jacobian) and $\mathbf r$ is the residual. In 1D we can often form $\mathbf J$ explicitly; in 2D/3D we apply $\mathbf J$ and $\mathbf J^\top$ implicitly via forward/adjoint solves—critical for memory and speed.

The quadratic misfit is optimal for Gaussian noise. With spikes or static shifts, replace the squared norm by a robust penalty $\Phi_d^{\rho}=\sum_i \rho\!\left((\mathbf W_d\,\mathbf r)_i\right)$, e.g. Huber or Tukey. This preserves sensitivity near zero but limits the influence of outliers, improving convergence and model stability.

• Set reasonable error floors in $\mathbf W_d$ (e.g., 5–10% amplitude, 1–2° phase).
• Choose $\mathbf W_m$ to match geology (smooth vs. blocky), and supply $\mathbf m_{\mathrm{ref}}$ when you have trusted prior information.
• Schedule $\beta$: start high, reduce by a factor 0.5–0.8 until $\Phi_d$ hits the noise level.
• Work in $\mathbf m=\log\sigma$ to enforce positivity; apply bounds via projection if needed.

  We now turn the objective into an algorithm. The core idea is simple:   linearize the forward map around the current model, solve for a model update   that reduces the objective, and repeat. This leads to the   Gauss–Newton / Levenberg–Marquardt (GN/LM) system:

  The matrix on the left blends three ingredients: the data curvature   $\mathbf J^\top \mathbf W_d^\top \mathbf W_d\,\mathbf J$ (how changes in the model   affect data), the regularization curvature   $\beta\,\mathbf W_m^\top \mathbf W_m$ (how changes affect roughness), and the   damping $\mu\,\mathbf I$ which controls step size. When $\mu\to 0$ we recover   pure Gauss–Newton; when $\mu$ is large, the step behaves like a cautious,   well-conditioned gradient descent. After solving for $\Delta\mathbf m$, we   accept a scaled step with a line search or trust-region rule and update   $\beta$ and $\mu$ adaptively.

  In Step 2 we will derive this system carefully and implement it in   practice:

  By the end of Step 2 you will have a working GN/LM loop that hits a target   misfit (discrepancy principle), respects bounds in $\mathbf m=\log\sigma$, and   produces models that are both geologically plausible and quantitatively justified.

e turn the objective into an algorithm. Linearize the forward map around the current model: $\mathbf F(\mathbf m+\Delta\mathbf m)\approx \mathbf F(\mathbf m)+\mathbf J\,\Delta\mathbf m$. With this quadratic approximation, minimizing $\Phi$ yields the normal equations:

The right-hand side is the gradient $\nabla\Phi(\mathbf m)=\mathbf J^\top\mathbf W_d^\top\mathbf W_d\,\mathbf r +\beta\,\mathbf W_m^\top\mathbf W_m(\mathbf m-\mathbf m_{\mathrm{ref}})$ with residual $\mathbf r=\mathbf F(\mathbf m)-\mathbf d$. The three terms on the left are the data curvature, the regularization curvature, and the damping $\mu$. Small $\mu$ gives a Gauss–Newton step; large $\mu$ behaves like a cautious gradient step.

Solve for $\Delta\mathbf m$ (usually with CG/LSQR, matrix-free). Then choose a step length $\alpha\in(0,1]$ by an Armijo backtracking line search, enforcing bounds in $\mathbf m=\log\sigma$ if present. Trust-region LM variants adapt $\mu$ using a gain ratio:

If $\rho$ is high (the quadratic model predicted the decrease well), decrease $\mu$ and allow bolder steps; if $\rho$ is poor, increase $\mu$ and retry with a smaller step. This keeps the iteration stable even when the linearization is imperfect.

Start with a larger $\beta$ to stabilize early iterations, then reduce $\beta\leftarrow c\,\beta$ with $c\in[0.5,0.8]$ once the misfit drops. The target is the noise level: $\Phi_d \approx N_{\text{data}}/2$ for a good fit.

This 1D loop is intentionally compact: it mirrors the exact logic we will reuse in 2D/3D, where matrix-free products replace explicit matrices and a preconditioned CG solves the normal equations. Next, we will plug this engine into a 1D layered forward model to recover a simple conductive basin—then scale the same algorithm to 2D and 3D.

et's put the algorithm from Step 2 to the test. We'll define a simple 3-layer "true" model: a resistive layer (100 $\Omega\cdot$m), followed by a conductive basin (10 $\Omega\cdot$m), over a resistive basement (1000 $\Omega\cdot$m). We generate synthetic data from this model using our 1D forward solver, add 5% Gaussian noise to simulate a real survey, and then feed these "observed" data to our inv1d_gnlm function. We'll start the inversion from a simple, uninformed 100 $\Omega\cdot$m half-space.

The solver converges in just a few iterations. It starts with a high $\beta$, finding the coarse structure first (as seen in the early iteration models). As the data misfit $\Phi_d$ approaches its target (the discrepancy principle), the algorithm automatically reduces $\beta$, allowing the inversion to fit the finer details. The final recovered model clearly resolves the conductive basin, demonstrating that our optimization loop is stable and effective.

ere is the most powerful idea in numerical inversion: the optimization algorithm we wrote for 1D does not change when we move to 2D or 3D. The inv1d_gnlm.py pseudocode, with its while loop, Armijo search, and $\beta$ schedule, is the exact same engine we use for 3D inversions with millions of parameters.

The only thing that changes is how we compute the four key components:

e've built a complete mental model for EM inversion. We started with the forward physics (Eq. 2), defined a stable, weighted objective function (Eq. 5), and derived the workhorse Gauss-Newton / Levenberg-Marquardt algorithm (Eq. 7) to minimize it.

The key insight is that the optimization logic is separate from the physics implementation. By using adjoints for gradients and matrix-free solvers for the update step, the same algorithm scales from a 1-second 1D model to a 10-hour 3D inversion. This modularity—separating the minimize(Objective) from the Objective.compute(gradient, residual)—is the core of all modern geophysical inversion frameworks.