I Built a 3D Electromagnetic Ocean
The last time I wrote an update on the electromagnetic (em) induction solver, I’d finally added the ocean as a thin sheet wrapped around the Earth. That was enough to calculate how much electromagnetic power went into the ocean vs the solid Earth, but we need to know where in the ocean that power is actually dissipated.
So, that’s where we are now, in phase 3 of development, and of course, it’s much bigger than I expected.
Where does that power go?
In phase 2, we were just asking how much electromagnetic power goes into the ocean. In phase 3, we want to know where that power goes relative to three dimensions: stay on the surface? penetrate deeper? does bathymetry change it? do coastlines change it? does the conductivity of the water matter? And, the question that started ALL OF THIS: how much of it gets routed into the abyssal ocean.
So, now the solver has an ocean with depth, it has a seafloor, it has land, it have conductive rock underneath the water, and electric currents can move through the whole coupled ocean-Earth system.
The physics I’m solving is simple:
∇ × E = iωμH
Ampère / conduction
∇ × H = J = σE
And the quantity we ultimately care about is the electrical energy being turned into heat:
q = ½ σ|E|²
A changing magnetic field creates an electric field, and that electric field can drive current through conductive material. Resistance turns some of that electrical energy into heat.
The challenge is computing this around an entire planet… and that’s what I’m doing.
A: Make it faster
Ok, my first problem is, this gets computationally expensive fast. Spherical electromagnetic calculations get ugly once the resolution starts climbing, so I built some acceleration into the solver and tested a Rust backend for the expensive transforms.
I mostly wanted to track where the time was going. And, I learned something…rewriting everything in a faster language isn’t the answer. It’s true, some pieces got dramatically faster, but the full solver really didn’t.
So Python is still running the scientific workflow, and I only accelerate the parts where it actually buys me something.
B: Give the ocean depth
This was the first big physics jump, because until here, the ocean was being treated like a thin conductive sheet. Obviously the real ocean has depth, so I gave it one.
I started with the most boring ocean possible:
Ocean depth = 4,000 m
Seawater conductivity = 3.2 S/m
Yes, this means the same depth everywhere, and the same conductivity everywhere. Unrealistic, I know, but it’s the perfect validation case, because I already knew roughly what the answer should be.
That meant I could test the electric field, current density, Joule heating, continuity across the seafloor, and the total electromagnetic energy entering the system. And, guess what?, the energy accounting closed :)
C: Real ocean
We know real seawater conductivity changes with temperature and salinity, so next I let conductivity vary with depth.
Current density is
which means two water columns can see almost the same electric field and still carry different currents if their conductivity profiles are different. That changes where the heating happens.
D: Real seafloor
Next, I added actual global bathymetry from ETOPO so every ocean column gets its own depth. If there’s a deep basin, then it’s a deep water column, and with a continental shelf, it’s shallow. And, conductive Earth continues underneath the ocean.
Sounds simple, but it’s really not when implementing :(
E: Global Maxwell problem
At this point, I had depth, conductivity, bathymetry, and geography, but if I solved every ocean column independently, currents couldn’t move sideways through that volume. :(
The actual Earth is globally coupled, if you change conductivity somewhere, you can change the electromagnetic solution somewhere else. So, Phase 3E became the full volume Maxwell solve. The solver carries the tangential electric and magnetic fields through the spherical shell and reconstructs the rest of the electromagnetic state from Maxwell’s equations.
[EB, EC, HB, HC]
The shallow modeled shell extends through the ocean and into the upper solid Earth:
Depth = 12 km
Next, I ran the real global case with ETOPO bathymetry and geographic SeaSigma conductivity. And it closed!
Ocean power ≈ 1,352.68 W
Solid-Earth power ≈ 1,057.63 W
Total power ≈ 2,410.31 W
Ocean fraction ≈ 56.12%
What I found interesting is that the ocean & solid-Earth split landed really close to what I got with the much simpler thin sheet solver. That’s pretty cool because the solver actually knows the ocean has depth, and it knows where the seafloor is. So, now we can calculate heating inside the water column.. a pretty big milestone.
F: Validation
Getting an answer and proving the answer is numerically stable are not the same thing…. welcome to phase 3F 😭
This phase has been humbling, and I don’t like it. The fields are represented with spherical harmonics and low harmonic degree captures the broad planetary-scale structure while higher degree gives the solver smaller and smaller spatial features.
My degree control is ‘l_max’.
larger lmax = smaller spatial features the solver can represent
Eventually, if the solution is correctly resolved, I should be able to keep increasing l_max and have the scientific answer barely move. That’s what I’m testing now, and my early convergence tests looked bad.
Total global power was slowly converging, but the heating maps kept changing way too much. At first I thought okay... maybe the EM solution just needs a lot more resolution. But, no…every time I increased the spherical-harmonic resolution, I was also sampling the Earth onto a different physical grid.
So I was changing the electromagnetic resolution...and slightly changing the Earth underneath it. Not good, guys, not good. 😭
I stopped the compuation and rebuilt the material side from the original high-resolution ETOPO and SeaSigma data. Now, I can sample the Earth once and freeze it so it has the same bathymetry, conductivity, ocean, forcing, vertical grid. Only l_max changes now, which is much cleaner. But.. still not converging.
Weirdly, the total watts look pretty good, and this is the part I’m glad I caught. If I’d only looked at global power, I could very easily convince myself this thing was basically done, (which I really want to do), but it isn’t.
At the latest completed resolution comparison:
Total-power change ≈ 0.048%
Ocean-power change ≈ 0.232%
4–6 km ocean-power change ≈ 0.263%
Surface electric-field map change ≈ 6.25%
4–6 km heating-map change ≈ 6.34%
Heating change at 5 km depth ≈ 6.07%
So total power is barely moving (this is excellent), but the actual geographic heating pattern is still moving by ~6% (not excellent). If all I wanted was the total electromagnetic power entering the system, I’d be feeling pretty good. But, I care about where the heating goes, especially in the deep ocean, so total watts converging is not enough.
And yes... the abyssal ocean is showing up. I’m already tracking modeled Joule heating by ocean depth, I’m just not comparing it to the observed abyssal warming yet. Not until the map stops moving.
One number has stayed interesting through the current high-resolution runs:
Ocean Joule power in the 4–6 km layer ≈ 11.1% of total modeled ocean Joule power
That caught my attention, and I want to be really clear about what it means. I’m not saying electromagnetic induction explains the abyssal warming (yet), I haven’t even done the final magnitude comparison. Right now, I’m still asking whether I trust the geographic heating solution itself.
First the map has to stop moving, then we can ask whether it matches anything real.
So, I bet you’re wondering, why is convergence so slow? This is always on my mind, why is everything so slow... A lot of it looks related to the real Earth being full of sharp boundaries like coastlines, seafloor, water suddenly becoming rock.
Spectral methods don’t really love sharp discontinuities, like Fourier ringing around a hard edge. Same general problem. A decent amount of the remaining deep-ocean error is clustering around those boundaries. The fields are still changing away from the boundaries too, so I can’t just mask out the coastlines and say good enough. Unfortunately 😂
I need to keep going…
Where I am now
I’ve finally figured out why the solver has been fighting me so hard around coastlines and the seafloor. The Maxwell physics itself still looks good, it’s the sharp boundaries.
Ocean water conducts electricity, but rock conducts very differently. So at a coastline, or right at the seafloor, the conductivity can change almost instantly. And, spherical harmonics really do not like that 😂
They’re great at representing smooth fields around a sphere. So, this morning I started testing ways to handle those boundaries better without changing any of the actual physics. First I tried just using more points. I evaluated the material interaction on finer and finer grids to see if brute force resolution would make the problem disappear. It definitely helped, but not enough.
Even going from an already very fine O6 grid to O8 still changed the material operator by around 1.7–2%. That’s too much, so then I tried something more targeted. Instead of increasing resolution everywhere, I only added extra numerical resolution near the places where the material actually changes… the boundaries.
That worked really well on some of my test geometries, a simple seafloor step converged beautifully. Then I tried a rectangular ocean patch... and of course it decided to be weird. As I increased the quadrature resolution, the error didn’t decrease smoothly, it bounced around.
But here’s the interesting part, at very high resolution, it still became extremely accurate. That gave me a pretty good clue about what’s happening. I don’t think the global spherical-harmonic solver is fundamentally the problem. I think I’m still asking the numerical integration points to figure out which side of a sharp boundary they happen to land on. Basically, I’m sampling a discontinuity and hoping enough points eventually describe it correctly.
So I’m changing the approach again… Instead of continuing to sample a cell with a boundary more and more densely, I’m going to find the boundary first, split the cell along it, and integrate the two materials separately. (lightbulb moment, guys)
So if one numerical cell contains both ocean and rock, the solver won’t treat it like one mixed cell anymore. It’ll know, this part is water, this part is rock, and it’ll integrate them separately. This is phase 3H!
If this works, I should finally be able to lock the numerical resolution. And, we can finally answer… how much of this induced electromagnetic energy actually reaches the abyssal ocean?
