ENH: Implementing current density as output of empymod

[prisae]

Message from @efinden

To understand more of the free and polarization currents in the ground it would be useful to adjust the output of empymod simply by multiplying the electric field by the eta H-function, in order to obtain the total current density in the frequency domain and then use model.tem to transform it to the time-domain.

In other words, just to use Ohms law in the frequency domain: J(w) = E(w)*etaH(w).

[prisae]

Thanks @efinden for this suggestion. I made a very rough proof-of-concept implementation in #220 - implemented, but NOT TESTED yet. Please have a go and try it out, I'd appreciate feedback!

You can install it via

python -m pip install git+https://github.com/emsig/empymod@ecurrent 

Note: No idea if this is anywhere close to correct. I will have to read into it more, and also compare it to other results. It was just a quick try!

[efinden]

Thank you very much for the new implementation!

I did some testing today and have a question on which locations in space that etaH are computed.

Is the etaH working on the output of the electric field for every receiver-location?

I tested one case with resisivity=200 Ohm-m, for a halfspace over air of resistivity 1e40. The air layer has a nearly zero conductivity, and thus etaH for the air layer should be very low (I used res=1e-40).

The configuration is a rectangular bipole source in the air layer at z=-0.1 m. Electric source. Electric receiver. Receivers are spaced out between x,y,z = [-5,5]. The current density of the air layer is very low, as expected, but the current density is also very low for the 200 Ohm-m halfspace.

I would perhaps expect the y-component of the electric field multiplied with the etaH-function in the ground of 200 Ohm-m to be much larger. Maybe it is something in my missing understanding of the implementation that is causing this results.

[prisae]

Good point. I multiplied it with the etaH from the source layer (as that is what is used in loop with zetaH), but it does not make sense here.

Then I have to put it into the loop where it goes over receiver depths. I will push an update.

[prisae]

I pushed a change, could you test again? Now it should multiply it with the etaH from the receiver layer

[efinden]

Thanks! It is fine, then I understand it as well!

A loop over the receiver depths will make it!

I will run a couple of tests, it takes some time with the filters I currently use and a lot of rec. positions, so perhaps I will not report here before tomorrow!

Looking forward to try!

[efinden]

A test of a pure resistive 200 Ohm-m halfspace, in the diffusional model (which means setting the permittivity to zero) is done. When asking for receivers to show the y-component of the electric field, and with the ecurrent = True, the modified code seems to reproduce exactly the current density as obtained when the conductivity is multiplied with the electric field outside empymod (on the output of empymod).

Thank you very much!

I look forward to investigate the polarization currents now for an IP case. Also, I will try to model a larger domain to see if the "smoke ring" of currents caused by the TEM system behaves as expected with depth.

[prisae]

Great! As I said, I only implemented according to your description - I haven't done any background research/read-up yet or thought deeper about it, so it is good to keep being curious about the feasibility of the outcome.

[efinden]

Thank you! I did a bit of background check on the physics before I asked and it seems all right, but other considerations are always welcome.

I post below an image of the current density due to a TEM switch off from a 2x2 m loop lying 10 cm above the surface for a 200 Ohm-m halfspace. Only some points in space are plotted, and the surrounding points are given the same intensity in the pixels, therefore it looks like the surface layer has zero current, this is only due to the sparse plotting.

In this case I use the quasi-static approximation and force a constant permittivity to zero, therefore there is only free current shown, despite the title says total current density.