Borehole Constraints

[geojoeyphillips]

Hey Matt,

This is a great package and I think I got the inversion to work for the oceans reflecting a crustal thickness variation consistent with our tomography. Working on another project on a smaller scale (~100m lateral resolution) with borehole information on the layer of interest at an elevation of ~400 m. For the high resolution gravity survey there's only 1 borehole and I was curious if there's a way to extremely force the inverted topography to fit the constraints? I tried added a more tapered function instead of the minimum distance one:

`def tapered_mindist(
points: pd.DataFrame,
grid: xr.DataArray,
radius: float,
region: list[float] | None = None,
) -> xr.DataArray:

grid = copy.deepcopy(grid)
original_dims = list(grid.sizes.keys())
constraint_points = points.copy()
min_dist: xr.DataArray = dist_nearest_points(
    targets=constraint_points,
    data=grid,
    coord_names=(str(original_dims[1]),str(original_dims[0])),

).min_dist
tapered_mask = xr.where(min_dist < radius, min_dist / radius, 1)

if region is not None:
    df = vd.grid_to_table(tapered_mask)
    df["are_inside"] = vd.inside(
        (df[original_dims[1]],df[original_dims[0]]),
        region=region,
    )
    new_tapered_mask = df.set_index([original_dims[0],original_dims[1]]).to_xarray()
    new_tapered_mask = xr.where(new_tapered_mask.are_inside, new_tapered_mask, 1)
    new_tapered_mask = xr.where(tapered_mask.isnull(), np.nan, new_tapered_mask)
    tapered_mask = new_tapered_mask

return tapered_mask`

That didn't quite work either so curious on your thoughts. Also, the topography of interest outcrops at the surface. I tried adding those surface values to the constraints but I think the density contrast I provide (sediments over crystalline rock 2700 - 1800 kg/m^3) doesn't reflect the surface points. Maybe that's why I can't get the inverted topography to match the outcrop? Looking forward to your input and thanks again for sharing this package!

[mdtanker]

Glad it's working out for you! So the weighting grid is not the ideal solution to constraining the inversion, but it was the only thing I was able to come up with for now. There are a few better alternatives I think to try out first.

For the outcropping parts, I'd try using the parameter upper_confining_layer within run_inversion. This takes a data array of topography and doesn't let the inverted topo cross above this. So this should probably be your surface topography.

For making the inversion stick to the borehole data, I would use the regional separation tools, as shown in this user guide. Ideally, the residual misfit that is input to the inversion is 0 at the places where you know the depth to the layer. This is the best way of stopping the inversion from moving the topography at these points. I'd try the DC shift method, passing a constraints_df. This should calculate a flat regional field resulting in a residual of 0 at the borehole. (The DC shift was only added in v.0.5.0)

Is your starting topography a flat layer at the depth of the 1 constraint?

Let me know if any of this helps!

[geojoeyphillips]

Interesting I totally missed the upper confining layer bit, thank you!

Regarding setting the regional grav as the borehole, I don't fully get how that works if I already removed a regional trend? I used the verde trend feature for the regional separation (I will test the 4 options you have as well!). By setting the regional grav to a point the forward problem will keep the borehole fixed but everything else will vary then?

[mdtanker]

That good to know you missed it, makes me realize I don't think it was in any of the examples.

The regional geologic signal can be removed earlier on in processing, as it seems you have done, but the regional I'm referring to is the regional component of the misfit. This inversion is based on a gravity misfit, where the misfit is defined as the difference between the gravity effect of the true topography and that of the starting model. This misfit can be separated into a regional misfit and a residual misfit. The residual misfit is what the inversion uses, so any areas with positive residual misfit values will try and raise the topography and vice-versa. Therefore, the residual misfit, by definition, should be 0 where we know the true topography (at the borehole), and you can use the regional to make sure this is true.

Hope that helps.

[geojoeyphillips]

You do have an example of upper confining layer in the regional separation example! I did not fully appreciate what it meant the first time I went through.

If you'd like for me to test the edit to regional.regional_dc_shift I can load in the function individually.

For a bit of clarification about setting grav_df["reg"] = constraints_df.sampled_grav I previously set grav_df["mistfit"]= grav_df["observed_grav"] and grav_df["res"] = grav_df.misfit should I not just set the grav_df["res"] at the borehole to be 0 then and not grav_df["reg"]?

[mdtanker]

Yes that would be great if you could test it! I have just implemented the fix in PR #74, which should be available now as v0.6.0 on pip, and should be on conda in the next day or 2.

This new version has several breaking changes, so make sure to look at the changelog to see how you'll need to update your code.

The gravity misfit is defined as misfit = observed - forward and misfit = region_misfit + residual_misfit. By setting the regional equal to the sampled misfit value at the borehole, you are setting regional_misfit = misfit (at the borehole) and therefore residual_misfit = 0 at the borehole. I hope this makes some sense.

Please ask me for more clarification if you need it, it's really helpful to hear your opinions and how you're interpreting the package so I can make the docs clearly in the future!

[geojoeyphillips]

Thanks for the update! The regional_dc_shift worked to force the borehole at that depth. Thanks again for doing that so quickly!

I plan to keep using this package for our various projects so I am excited things are working out so well.

I have been curious about one other thing. In the case of sedimentary basins where you could have layers pinching out or in the case of alluvium covering a basement unit in one part of the target region and then a lava flow in the other. Do you think there's a way to set up a grid with laterally varying density contrast to find a target topography layer? Maybe that's too complicated but I'm curious if you have thought of this!

Thanks again!

[mdtanker]

Glad to here! If you had any figures of the inversion outputs or anything that you can share I'd love to see how its working.

I've also been thinking about laterally varying density contrasts. I think it should work fine, you'll just need to use a slightly more complex xarray.where() function for defining the density grid, especially if you uses a zref value that is not entire above or below your surface. Let me know if it works alright.

[mdtanker]

Following up, Invert4Geom now supports spatially variable density contrasts values. See this example for how it works.

[geojoeyphillips]

Woah that's great! I wrote a complicated xr.where() but this looks much simpler. The shared example says that the density contrast will be 1799 or 2999 if upward less than or equal to 0? So this is in place of zref then? Thank you!!