gravity forward modelling: block model, save grid and gravity

[brianchen123]

I am new to Harmonica and I am trying to build a tilted block model for gravity forward modeling. Could anyone provide some sample codes? The codes that I used seem to be outdated and are producing errors. Also, where can I find the latest documentation for Harmonica?

Thank you,
Brian

import numpy as np
import verde as vd
import matplotlib.pyplot as plt
import harmonica as hm

Define the grid of the desired size and resolution

grid = vd.grid_coordinates(region=[-100, 100, -100, 100], spacing=1)

Define the geometry of the tilted block

block_model = hm.utils.make_tilted_block(grid, slope=20, direction=-45, thickness=20, density=2.5)

Create a 3D array for the final model and initialize it with the density of the block

density_model = np.ones_like(block_model) * 2.5

Define the indices of the cells below and above the tilted block

below_indices = np.where(grid[:, 1] < 20 * np.tan(np.deg2rad(20) * 0.5))[0]
above_indices = np.where(grid[:, 1] > 20 * np.tan(np.deg2rad(20) * 0.5))[0]

Set the density of the cells below and above the block to 2.7 and 2.5 g/cm^3, respectively

density_model[below_indices] = 2.7
density_model[above_indices] = 2.5

Create a Harmonica Gridder object and add the density model to it

grd = hm.Gridder(grid, data=density_model)

[MGomezN]

Hi @brianchen123! The latest Harmonica version is Version 0.6.0, Released on: 2023/03/01. This and the previous documentation are available here check the "change log" section. Also, under the "gallery" section you will find multiple examples.

[brianchen123]

Thanks

[LL-Geo]

Hey @brianchen123,

I don't think there is a function called utils.make_tilted_block and Gridder in Harmonica...

Sadly, we can't directly handle the 3D block model at the moment. We do have a discussion last year to get a 3D block function, but it not happening atm. One thing you could try is to build multi prism_layers from top to bottom. Then you can define the density based on its 3D locations. In the end, you can calculate the gravity response from these layers using prism_layer.gravity and add them together.

An alternative solution is to use Simpeg which is easier to generate and control the mesh, see the example in here.

[brianchen123]

Thanks Lu. This is helpful. I'll create three prism layers, which should be a common case for a thin block structure. It's a simple case but not easy to build using either Harmonica or SimPEG. It would be nice to use a GUI for model creation, if one exists.