Wednesday, March 22, 2017

Geological modeling in python

Latest results will be posted here at the top as progress warrants...

==== 22 Mar 2017 ====

Now we are getting somewhere! Improvements here include long and short wavelength variation on basement and Mississippian unconformities, realistic level of noise added, two Penn limestone beds and a late monocline down to the left. The monocline is the first function for a new function called bender that will include monoclines, bumps and dip. Remember, the idea is not to reproduce this particular field line in all it's detail, but to reproduce essential characteristics and trace them back to the earth model. For example, the shallow Penn LS event on field data is a seismic thin bed as evidenced by a 90 degree phase shift. Either rock properties, thickness or both will have to be changed since this phase shift is not seen in model data. Also, was able to load a model 2D line into OpendTect, leading toward the real goal...3D!!

==== 6 Mar 2017 ====
Low impedance weathered LS zone is now partially embedded in LS at the unconformity, give better fit to amplitude variation observed in real data.

Note isolated doublet at high point of unconformity surface that might be interpreted as a sinkhole is actually a LS outlier underlain by weathered LS

--------------- original 9 Mar 2017 blog entry ---------------------------

For some time I have been considering the difficulty of making realistic geological models for input to seismic modeling. There are no good tools out there, despite half of century of clear need. Oh sure, there are some computational geometry packages that can model anything from a space ship to the Rocky Mountains given enough time and money, and geostatistical packages that rarely make geologically realistic models. One approach would be to mimic nature and start with a granite surface, deposit sediments, follow a sea level curve, subside, uplift, erode, preserve, and so on. Appealing from a 'let's completely replicate nature' viewpoint,  but would require a vast number of parameters to control everything.

The goal, for me, is to generate 3D seismic data with a level of geological realism that makes it suitable for seismic interpretation training. After all, in this scenario the answer is known -- the geological earth model in depth and time consisting of P-wave, S-wave and density at every sample point. I am more in the camp of video gamers who construct complex and realistic landscapes with a bag of tricks and simple algorithms. The simplest code that creates to most realistic result, that is the goal.

The specific object of my current attention is US mid-continent data typical of NE Oklahoma where my U Arkansas geosciences department has several donated 3D seismic surveys. You can find an image of a typical seismic profile here. From bottom up, key elements are: (1) a granitic basement unconformity with significant relief and a weathered granite conglomerate of irregular thickness, (2) Cambrian through Mississippian carbonates, thickly layered and laterally consistent, (3) a major unconformity at the Miss-Penn boundary with rugged relief and a low velocity, low density, deeply weathered, irregular thickness unit sitting on the unconformity surface, and (4) Pennsylvanian sandstone units of laterally variable thickness and properties with shale between sandstone units with slightly vertically variable parameters. Here is an example:

Of particular interest is the issue of subtle structural distortion on seismic time data due to overlying stratigraphic lateral variations. This kind of distortion cannot be fixed by prestack depth migration unless the migration velocity field is exquisitely precise, which is not currently achievable. Perhaps full waveform inversion my deliver that kind of velocity model someday. But till then, we are left with subtle structures that may be false closures and tempting drill targets.

I wanted to do all this for research interest and as a python learning exercise. Depth-to-time conversion and some display are done with SeismicUnix, everything else in python. My first examples were pretty crude (Figures 1 and 2 below). Too simplistic for interpretation exercises.

After months of head scratching, python learning and geological sketches, the result is getting to a realistic level of complexity (Figures 3-5). I learned about bumpy, nymph arrays in N dimensions, loops, conditionals, and how to write a files. Still stumped on writing binary files that can be read by SeismicUnix, so everything is text files so far but I think this will not work in 3D. So more learning ahead. Simulation of seismic vertical resolution is accomplished by convolution with a wavelet and lateral resolution by smoothing to bin size before depth conversion and convolution. Still needs improvement, still only 2D, but a much more realistic result. With a bit more progress and the jump to 3D, these earth models can be used to study issues related to subtle false structure, vertical and lateral resolution, effect of sub resolution features, interpretation methods for stratigraphic targets and depth conversion for subtle structure.

The thought of doing all this in C is daunting to the point I would not have tried it, even with the SeismicUnix C environment as a framework. Long live python.

Figure 1. Earth parameter modell in depth from early effort in Nov 2016. The Miss-Penn unconformity and chat layer on it look pretty good though. This early version did not have a basement unconformity.

Figure 2. P-wave model in time and the noise-free seismic response for earth model in Fig 1. Note the subtle structure for layers below the prominent Miss-Penn unconformity, even though the depth model above shows these beds are horizontal.
Figure 3. Current code P-wave velocity model in depth. Now we have a basement unconformity, a basal conglomerate of variable thickness and many Penn sandstone layers cut into shale or earlier sandstone units.
Figure 4. Current code earth model in P-wave time.

Figure 5. Current code noise-free seismic response for earth model in Fig 4. The Penn section may be a bit overcooked, but that can be controlled by a few parameters. At lease we now see the characteristic wormy appearance typical of complex sand-shale sequences worldwide and false structure for deep events.


Ricardo Serrano said...

Good day Prof. Liner. I'm Ricardo Serrano and I'm the founder of a platform called Geomodelr. With Geomodelr not only you can create realistic models only interpreting cross sections, but you also can query them with Python with really simple operations.
The resulting model will be fully 3D and solid, that means that at every point you query there will be one and only one formation.
I'm really interested in your work and we actually used geomodelr to create a seismic simulation of a realistic 3D model:
Effects of realistic topography on the ground motion of the Colombian Andes - A case study at the AburrĂ¡ Valley, Antioquia.
If you want, we can talk about it and I can help you to generate any realistic model you want in hours, days top.

Best Regards, Ricardo.

Prof. Christopher L. Liner said...

Ricardo, sorry for the long delay. Thank you for the comment and offer to help build geological models. The point for me is not the need to make a specific model, but to think about the minimal set of python functions I need to make something useful. Also interesting to think about the niche of my kind of geoseismic modeling versus what is out there in open source or commercial code. It is the journey that has value for me. I will be posting progress soon and hope to see more of your comments.

Thanks again, Chris