search

swung-research / 3d-csem-open-source-landscape

Werthmüller, D., R. Rochlitz, O. Castillo-Reyes, and L. Heagy, 2021, Towards an open-source landscape for 3-D CSEM modelling: Geophysical Journal International, 227(1), 644--659
https://doi.org/10.1093/gji/ggab238
Creative Commons Attribution Share Alike 4.0 International
11 stars 5 forks source link

3D Model #10

Closed prisae closed 4 years ago

prisae commented 4 years ago

3D Model

You can find the 3D model in the notebook notebooks/MarlimR3D.ipynb.

Questions:

Any inputs on this model @Moouuzi, @ocastilloreyes, @lheagy ?

That is how the horizontal resistivities look like:

Selection_001

Marlim R3D

Papers

DATA

Model and survey lines

... TODO Model and Survey...

3D block model

... TODO Results ...

(Related to #4.)

prisae commented 4 years ago

@Moouuzi thanks for bringing this up (regarding mesh). We'll soon have a new model, I just received another email from Paulo:

We generated two new seg-y files with the 100 x 100 x 20 grid. I'll upload them soon.

prisae commented 4 years ago

As my grid-to-grid interpolation routine just continuous with the last values I didn't even notice this, so very glad you caught that!

Moouuzi commented 4 years ago

I've recalculated the previous solution with updated resistivities because I discovered some important artifacts during the resistivity interpolation -> Of course, the interpolation approach from the regular grid on my tetrahedral mesh produces some mismatches, those are most problematic at the sea-bottom and should be removed. Furthermore, I used the correct water resitivity of 0.32 Ohmm (instead of 0.3) :-)

I've also added another solution using reciprocity - Interesting, that the results are in agreement, but show different boundary effects.

I forgot to uncomment the "append boundary mesh" line. Thus, I'm pretty sure that the boundary effects will become way smaller in the next runs with the new model. Disregarding the explanation for E_z below, the only strange difference that remains is: what's wrong with the inline - E-y component?

results_p1_new

What I meant with the importance of "rx/tx" positions relative to the layer interfaces in the mesh:

Here we can clearly see the effects I often observed before: It's obvious that in land-based applications, 'E_z' varies significantly if received 1, 10 or 100 cm below or above the earth's surface due to the field discontinuity. I estimated these effects (and compared against empymod solutions), E_z varies over several orders of magnitudes if received from different heights between 0.001 and 1 m below or above the surface. This behavior is also very frequency-dependent.

The same holds for the water/subsurface interface, but the effect is of course smaller. Anyway, if I incorporate the digital elevation model via interpolation from the original dem, the provided "real" Rx or Tx coordinates are located either slightly above or slightly below the interface in the mesh due to the DEM interpolation - the finer the surface-mesh, the better the approximation, but never exact.

On this account, I always slighty shift Rx coordinates relative to my triangulated-interfaces, ensuring that Rx results are interpolated from the elements below the interface (or above). In the figure above, the black curves shows results from -849.7 m, the blue curve results from z = -849.8 m (original coordinate was z=-879.78m and was in the sea bottom leading to the E_z mismatch in the first plot). Even this 10 cm difference can lead to critical misinterpretations. If there are real slopes in the topography or bathymetry, the discontinuity is not aligned with E_z anymore but visible in all components! - I discovered this from a modeling of the Stromboli volcano.

In the reciproce test model, it becomes even more interesting: My transmitter is placed "exactly" on the interface, which means (considering the reciprocity and regarding it as Rx) - there is an ambiguity caused by the discontinuity: The field can be approximated either as limit from above or limit from below. This is exactly seen in the reciproce E_z fields, some values belong to results obtained above the sea-bottom, some from below. I assume the choice depends on numerical tolerances.

prisae commented 4 years ago

These are fantastic results. Yes, I agree with what you say regarding $E_z$. They are also the least important, as they are not actually measured in the field (AFAIK). And I don't think they will ever agree between FD and FE codes with topography. However, everything you said could be interesting to put into the discussion of the results.

prisae commented 4 years ago

This is in agreement with my first results I posted on Slack once. Ex agreed quite well, Ez was out by a factor of 10, which could be similar to the issue you describe, but Ey was also off. Funny. After my holidays I will model it also on a much denser grid, and if we have two results that show the same I will get back to them...

prisae commented 4 years ago

That is the result I posted on 23rd December on Slack:

index

The model was rather coarse, just containing 3 million cells.

prisae commented 4 years ago

Also, inline Ey is rarely used, it is a very weak field (but broadside Ey is used). To be honest, I don't trust their Ey, because their step-like behaviour is very suspicious. So if our results agree I think we have the good result there...

Having said that, my 'real data experience' is very limitied and old, from PGS time, where everything was time-domain and a completely different story, so I am learning as we go and I might be wrong ;-) @ocastilloreyes might have more experience with actually measured node-based CSEM data?

prisae commented 4 years ago

Here you can download the extended calculation model: https://cloud.werthmuller.org/index.php/s/pCWytp9nw6XmYPR

prisae commented 4 years ago

No wonder our large-offset results look different, given their strange 'extrapolation': index Anyway, our goal is not to judge their model, but to reproduce it. With this it should be better possible, I guess...

I updated the notebook https://github.com/prisae/3d-csem-open-source-landscape/blob/master/notebooks/MarlimR3D.ipynb on master to include both the core domain and the modelling domain.

(I'll be off on holiday from this afternoon onwards till 27/01.)

prisae commented 4 years ago

Yep, with this model it looks fine. Here is a very quick run (5 million cells, 250x250x100 m cell size), and it looks already very similar. Doing it with a finer grid should yield very comparable results.

index

Moouuzi commented 4 years ago

Thanks for the updates!

My p1 results with a boundary mesh and a finer discretizaion look even better for all frequencies. Regarding their conductivity model, I assume it is based on one of these extrapolation algorithms (Thomas has shown me an application with an algorithm tested in pyGIMLi ). Sometimes, they work quite well, but often, strange distributions are generated.

Anyway, I'm glad with the progress for now.

Dieter, have a nice holiday!

Moouuzi commented 4 years ago

Results with p1, calcualtion of all Tx and recirocity approach, for all frequencies, but not considering the new Extended conductivity model.

Computation times for all Tx: 2.5h / freq. Computation times using reciprocitiy apporach: 1h / freq.

results_p1_all_freq.pdf

Now, as a final step, I will add the new resistivity model and try to reduce the number of nodes for a reasonable p2 model to see, how much a more or less exact Interpolation of conductivities influences the results.

Moouuzi commented 4 years ago

The new resitivity model appears to be sorted in another way, at least it does not match anymore if I switch axes as before. To project in my right-handed coordinate system, i used, with the old segy data:

np.transpose(res_h[::-1, :, ::-1], (1, 0, 2)) and the x0 and y0 coordinates as derived in the notebook.

Now, I Imported the new data, used the new nx, dx, etc..., as in the notebook, but the resitivity distribution does not match to the layers anymore, and obviously, the results are worse now. Looking at the notebook, for me it seems that x/y resitivities might be switched on the plotted vertical slices. At least the new extended resistivity plot does not match the "yellowish" colors from the previous model anymore but rather the ones on the northing slice. I tried to switch x/y but that does not resolve the problem - x0/y0 might be also not correct anymore.

Maybe I'm mistaken and the bug is somewhere else, please correct me in this case.

Anyway, I prefer the fine resitivity model for interpolation reasons, so I tested to extent this model on my own in each direction, just using the values from the outermost slice in each dimension. That works quite well and leads to promising results. Might be even more realistic than their extrapolation results … plots follow.

PS: @Dieter, If you read this already in your vacation time, don't Reply:-) Not that urgent...

Moouuzi commented 4 years ago

Very good news: I used my own resitivity extension as described above and (log-transformed resistivites) linear 3D regulargrid-interpolation from scipy. My (final) Marlim version with p2, 30 min/freq, 1TB RAM @ 64 cores with reciprocity:

p2_rzp.pdf p2_rzp_errs.pdf

I think, overall ~10 % misfit is quite excellent considering all the differences in codes, discretization, conductivity distribution. Now I could also adapt a p1 version with probably even less resources, but for sure slightly unsmoother results.

@ Octavio: Do you prefer a mesh for p1 or p2? Do you need the mesh only or also some ideas for the resistivity interpolation?

prisae commented 4 years ago

Really great results @Moouuzi, looks awesome that we can reproduce the EMGS-code result so nicely, and also at a pretty great speed! Awesome! I will look into it in more detail by the end of next week (just back from holiday, and next week I have the consortium meeting).

I do think it is OK to do our own extrapolation if the results agree. Looking at the notebook https://github.com/prisae/3d-csem-open-source-landscape/blob/master/notebooks/MarlimR3D.ipynb they indeed seem to be different. To me it appears that the Northing-axis is flipped, judging from the bathymetry. Did you try that? (Hence, not swapping x/y, just inverting y.)

ocastilloreyes commented 4 years ago

Amazing results @Moouuzi. I think I should use the same conductivity model that you have worked on. In this way, the comparison would be fairer (not only among our codes but also with those we are using as a reference). How do you define your conductivity model? I only need a list of conductivity values for each mesh's element. Is this possible?

I will try to work on it in the first week of February. In the next few days, I will be in a project workshop ...

Moouuzi commented 4 years ago

@ Octavio: Unfortunately, I will not be able to prepare the input for you this week due to several other commitments. Probably, I will send you relevant files between next Wednesday and Friday.

ocastilloreyes commented 4 years ago

@Moouuzi no problem. This week I have returned to Barcelona and I am finishing some pending tasks, therefore I can wait

Good weekend!