Issue with prism_magnetic - unwanted asymmetry in solution

[lruepke]

Thanks so much for adding prism_magnetic!

Description of the problem:
I am playing with the new functions and fail to reproduce the classic 2D solutions for a burried rod in Heitzler et al. 1962 using Harmonica's new prism_magnetic function.

Setup:
Infinite rod of 10 km width and 15km vertical extent buried at 2 km depth that causes an induced field anomaly.

Expected behavior:
For an inclination of 45degree (in the northern hemisphere) and an angle of 90degree between strike of the rod and geographic north, I should get a symmetric total field anomaly with a low to the North and a high to the South. For a magnetization of 0.001emu (unit in heitzler, I guess it means emu/cm^3 equal to 1 A/m), the peaks are at +145 and -145-ish gamma in Heitzler.

Harmonica result:
I am using the development version of Harmonica from github. Version: v0.6.0.post33+gc677dfa
I use a rod that is aligned with the easting direction. A setup of declination=0 should reproduce the strike=90 of Heitzler.

I get total field anomalies as functions of strike and inclination that are very similar to the Heitzler results. But, at inclination=45 and strike=90 (declination 0 in the harmonica setup) the predicted anomaly is not fully symmetric and the magnitude is also different at +303nT and -292 nT. For other symmetric combinations of strike and inclination, I also get asymmetries.

Possible solutions:
Most likely I am doing something wrong - any help would be highly appreciated! Might also point to an issue in the kernel functions or wrappers - that's why it made it a bug report. Hope that's ok.
I am unsure about the unit conversions but as far as I understand code, the 4s and pis are taken care of inside choclo and harmonica.
I am also unsure about the projections of b_e, b_n, b_v to the total field but can't see anything wrong there.

Thanks!
lars

Full code that generated the error

import numpy as np
import harmonica as hm # this is the development version build from source
import matplotlib.pyplot as plt

easting = np.linspace(-13e3, 13e3, 201)
northing = np.linspace(-13e3, 13e3, 201)
easting, northing = np.meshgrid(easting, northing)

# in radians
inclination = np.deg2rad(45)  # in radians
declination = np.deg2rad(0)  # in radians

# Compute induced magnetization
magnetization_strength = 1 #A/m
magnetization = [
    magnetization_strength * np.cos(inclination) * np.sin(declination),
    magnetization_strength * np.cos(inclination) * np.cos(declination),
    -magnetization_strength * np.sin(inclination) # minus because vertical axis is positive upwards in harmonica
]

prisms = np.zeros((1,6))
prisms[0, 0] = -50e3 # west - make large horizontal extent
prisms[0, 1] = 50e3  # east
prisms[0, 2] = -5e3  # south
prisms[0, 3] = 5e3  # north
prisms[0, 4] = -17e3  # bottom
prisms[0, 5] = -2e3  # top

# Compute magnetic anomaly
coordinates = [easting.ravel(), northing.ravel(), np.full_like(easting, 0).ravel()]
anomaly = hm.prism_magnetic(coordinates, prisms, magnetization)
shape = easting.shape

# get total field anomaly from b_e, b_n, b_v
anomaly_bu =  anomaly[0]*np.cos(inclination) * np.sin(declination) + \
   anomaly[1]*np.cos(inclination) * np.cos(declination) - anomaly[2]*np.sin(inclination)

# Reshape the anomaly to match the original grid shape
anomaly_reshaped = anomaly_bu.reshape(shape)

# Plot
fig, axs = plt.subplots(1, 2, figsize=(16, 6))

pc = axs[0].pcolormesh(easting, northing, anomaly_reshaped, shading='auto', cmap='RdBu_r')
fig.colorbar(pc, ax=axs[0], label='Magnetic Anomaly (nT)')
axs[0].set_xlabel('Easting (m)')
axs[0].set_ylabel('Northing (m)')
axs[0].set_title('Computed total Magnetic Anomaly')

cross_section_index = len(easting) // 2
vertical_slice = anomaly_reshaped[:, cross_section_index]  

axs[1].plot(vertical_slice, northing[:, cross_section_index])
axs[1].set_xlabel('Magnetic Anomaly (nT)')
axs[1].set_ylabel('Northing (m)')
axs[1].set_title('Max: {:.1f} nT Min: {:.1f} nT'.format(max(vertical_slice), min(vertical_slice)))

plt.show()

Full error message

System information

• Operating system: MacOS 13.6
• Python installation (Anaconda, system, ETS): anaconda,
• Version of Python: 3.11.5
• Version of this package: v0.6.0.post33+gc677dfa
• If using conda, paste the output of conda list below:

Output of conda list

PASTE OUTPUT OF CONDA LIST HERE

[santisoler]

Hi @lruepke! Thanks for opening this issue, it's very appreciated!

Thanks so much for adding prism_magnetic!

You're welcome! We are also excited about it, and we are glad that some users are already testing it. That shows that we should be making a release soon also.

Regarding the issue you are raising. I don't have Heitzler et al. 1962 at hand, and couldn't find the right article for this reference. Does it have a DOI that would make finding it easier? Or maybe a link to it (even if it's behind a paywall, I can get it from my University's library)?

By reading your setup and your code, I was a bit confused by the term "rod". I always think of a rod as a cylindrical body with a length much larger than the dimensions of its circular section. In your case, it seems more like a vertical dike, right?

Regardless of that, I think your point is should still be considered.

After doing some math, I think the problem here is to assume that the total field anomaly should be antisymmetric in this case. I don't think it is.

For simplicity, I assumed a horizontal infinite cylinder (rod) whose axis is parallel to the easting and with a arbitrary magnetization $\mathbf{M}$. Blakely (1995, p. 96) gives the solution of the magnetic field $\mathbf{B}$ on any point $\mathbf{p}$ outside the cylinder:

$$ \mathbf{B}(\mathbf{p}) = \frac{m}{2\pi \mu_o r^2} \left[ 2 ( \hat{\mathbf{m}} \cdot \hat{\mathbf{r}} ) \hat{\mathbf{r}} - \hat{\mathbf{m}} \right] $$

where $\mathbf{m} = \pi a^2 \mathbf{M}$ ($a$ being the area of the cross section of the rod), $\mathbf{r}$ is the vector that goes from the axis of the rod to the observation point $\mathbf{p}$ and it's perpendicular to the axis. The "hat" version of these vectors are unit vectors, and $m$ and $r$ are their magnitudes.

The total field anomaly can be computed as:

$$ \Delta T = \hat{\mathbf{F}} \cdot \mathbf{B} $$

where $\hat{\mathbf{F}}$ is the unit vector related to the Earth's background magnetic field.

If we assume that the rod has only induced magnetization, then its magnetization vector is parallel to $\mathbf{F}$, so:

$$ \hat{\mathbf{F}} = \hat{\mathbf{m}} $$

So we can express the total field anomaly as:

$$ \Delta T(\mathbf{p}) = \frac{m}{2\pi \mu_o r^2} \left[ 2 ( \hat{\mathbf{m}} \cdot \hat{\mathbf{r}} ) \hat{\mathbf{m}} \cdot \hat{\mathbf{r}} - \hat{\mathbf{m}} \cdot \hat{\mathbf{m}} \right] = \frac{m}{2\pi \mu_o r^2} \left[ 2 ( \hat{\mathbf{m}} \cdot \hat{\mathbf{r}} )^2 - 1 \right] $$

Now, let's consider that the rod is parallel to the east direction, located at northing = 0 and at a depth of $h$. And consider an observation point $\mathbf{p}$ located at easting = 0, at zeroth height and at a northing coordinate of $y$.
In this scenario we have that:

$$ \mathbf{p} = (0, y, 0) $$

$$ \mathbf{r} = (0, y, h) $$

And consider that $\hat{\mathbf{m}} = (0, m_y, m_z)$.

In this scenario, the total field anomaly on $\mathbf{p}$ can be calculated as:

$$ \Delta T(y) = \frac{m}{2\pi\mu_0 r^2} \left[ 2 (m_y y + m_z h)^2 - 1 \right] $$

It's not hard to proof that this expression of the total field anomaly is not antisymmetric with respect to $y$ if $m_y \ne 0$. It does shows an antisymmetric behaviour, but it's not antisymmetric across $y=0$, so $\Delta T(y) \ne -\Delta T(-y)$. It does exists a value $\alpha$ that satisfies that $\Delta T(y + \alpha) = -\Delta T(-y + \alpha)$, I haven't done the math to get an expression for it though.

By extending this analytical result to your model, I think the code is right: the total field anomaly should not be antisymmetric across $y=0$ for the 45 degree inclination case (and not for any inclination different than 90 or -90).

Let me know if this is clear enough and if there's any more questions related to this case. Looking forward to your reply.

[lruepke]

...thanks for the detailed answer! I need to "digest" it and will then write again. Just for reference, the Heitzler report has this doi: https://doi.org/10.7916/d8-gx40-qf60 and can be found here: https://academiccommons.columbia.edu/doi/10.7916/d8-gx40-qf60. I find it an accessible paper on the Talwani-type 2D solutions.

[lruepke]

I agree with your analysis of the Blakely solution that there is no antisymmetric behavior for I=45 (and strike of the rod of 90deg -> infinite in east/x direction).

Interestingly, the Heirtzler 1962 equations (2.1 and 2.2) for the horizontal and vertical anomalies induced by an infinite buried prism (is that a rod?) give an antisymmetric solution at I=45 (eqn. 2.9) - and thereby models based on that formulation. I am puzzled by this; maybe a simplification was introduced in the Heirtzler solution that I am not seeing?

[santisoler]

Thanks @lruepke for sharing the links to Heitzler's article.

I quickly read the first part of Section A of Chapter II (Derivation of the Basic Formulas) and I have the following observations:

• It assumes an arbitrary prism to start the derivations (fig 2-1). An arbitrary prism is not a rod. Usually a rod is considered to have circular section and its length is much larger than the radius of its section.
• When defining the magnetic potential of the prism (the P magnitude in the second equation of the section), they are applying some simplifications: they approximate the prism with magnetization $\bar{\mathbf{M}}$ as a dipole with magnetic moment $\bar{\mu} = \bar{\mathbf{M}} \Delta x \Delta y \Delta z$. Analytic solutions for the magnetic field of prisms with arbitrary magnetization exists though (see https://www.fatiando.org/choclo/latest/api/generated/choclo.prism.magnetic_field.html#choclo.prism.magnetic_field)
• I haven't done a complete analysis, but I cannot see that equations 2-1 and 2-2 are antisymmetric. In both cases, when evaluating them on $-x$, only the rectangular term ($2xzM_z$) changes sign, while the quadratic terms ($x^2$ and $z^2$) don't. Similar analysis could be made on eq 2-3 and 2-4.
• In eq 2-9 the $H$ and $V$ components are also functions of the location of the receiver, and they are not antisymmetric (as pointed out in the previous point). So, except in the case at which the magnetization is perfectly vertical or perfectly horizontal, $T$ is doesn't show symmetries around $x=0$.

I arrived to these observations after a quick read, so let me know if you think any of them is wrong.

Thanks again for raising this question!

[santisoler]

I'm removing the bug label from this issue, but I'll keep it open to continue the conversation.