[Question] Interpolation approach for stratigraphic units of varying thickness

[davenquinn]

I am evaluating Loop Structural on a well log dataset of the Williston Basin. This is a relatively uncomplicated basin structurally but the study area is large so sedimentation dynamics come into play, and some formations change dramatically in thickness across the basin.

I am trying to create a structural surfaces from a set of formation tops (so value constraints, no gradient information). The tricky bit is that I would like to treat formation thickness as a free parameter that varies across the basin. Instead, it appears that the finite difference interpolator (at least, with defaults) fairly rigidly enforces a constant gradient in the val parameter when a conformable stratigraphy is being modeled. Is there a way to relax this constraint?

Here is the result of Loop3D modeling:

If I do a naïve interpolation of each surface separately with SciPy, I get something much closer to the data, but where stratigraphic superposition is not obeyed:

I've tried using unconformities and interfaces, but I am not sure I was using them correctly in any case.

My code is here for reference.

My model data frame looks like this:

                                    X             Y            Z   val  interface feature_name
WELL_ID unit_name
711     TOP070_Pierre   -756901.954423  2.764497e+06   431.292863     0          0         main
779     TOP070_Pierre   -632910.888494  2.791280e+06   180.137160     0          0         main
784     TOP070_Pierre   -632515.256795  2.795109e+06   175.260351     0          0         main
786     TOP070_Pierre   -620580.894934  2.803847e+06   121.005842     0          0         main
789     TOP070_Pierre   -632946.220886  2.795394e+06   183.794768     0          0         main
...                                ...           ...          ...   ...        ...          ...
13919   TOP542_Basement -516713.285205  2.833342e+06 -3780.746761 -28           28         main
14009   TOP542_Basement -521056.245299  2.827089e+06 -3716.129032 -28           28         main
14016   TOP542_Basement -515309.978975  2.848040e+06 -3643.891288 -28           28         main
14029   TOP542_Basement -595338.113582  2.840625e+06 -3616.459233 -28           28         main
16098   TOP542_Basement -588366.037256  2.557178e+06 -2058.318517 -28           28         main

and my stratigraphy has this general structure:

{
'main': {
  'TOP070_Pierre': {'max': 0, 'min': -1, 'id': 0},
  'TOP097_Mowry': {'max': -1, 'min': -2, 'id': -1},
    ...
}

As you can see, I'm using arbitrary integer indices to track the order of surfaces; I'd expect the final interpolation to have a varying spacing between these value contours.

I apologize if I am missing some fundamentals of this package's approach, as I only just began to use it today. I'd appreciate any recommendations on the right approach to take for this particular problem.

[lachlangrose]

The interpolators generally will all minimise variation in thickness. If you reduce the regularisation this may be less.

We haven't really focused on modelling stratigraphic complexities as the tools for these are fairly well developed in commercial code. We have tried to add tools for incorporating more complex structural geometries e.g fault displacement informed by kinematics and folding.

Using unconformities will work, happy to have a look but won't be at a computer for 3-4 weeks.

[davenquinn]

Hi @lachlangrose thank you for the perspective. I suspect your framework could be adapted to this different style of interpolation, which might be an interesting problem in its own right (and not quite the same as going in a full stratigraphic modeling direction). But it's not urgent to solve.

I'll probably retreat to more well-supported use cases for now, but may reach back at some point if I poke more at this.

This module has been a pleasure to get familiar with overall!

[lachlangrose]

Hi @davenquinn,

I'm back from my holidays so happy to have a chat about these ideas.

Looking at your pictures again, the current way that I would build a model like this would be to combine multiple scalar fields interpolated independently (as you've done with scipy) and then use unconformities to control the overprinting. Effectively this is just a series of boolean relationships between them. This should give you comparable results to the scipy version.

If you wanted to take this even further you could use inequality values to constrain the geometries between the different layers but for the examples you have shown the contact location points are probably adequate.

Cheers