Thursday, April 27, 2017

Python, wells, LAS and info

Those of you who deal with wireline well log data know that the standard format is call LAS. A nice site from the Canadian Well Logging Society that explains the format can be found here. In the old days, a well would be logged for gamma ray, density and maybe a couple of electric logs. In today's world, a full logging assault can yield a vast collection of data types (called curves). Not only that, but these might be split across many LAS files. One well in our data base has 17 LAS files! Each individual curve has a short abbreviation (a mnemonic) such as GR (gamma ray), RHOB (bulk density) and DT (compressional sonic). The mnemonics are somewhat standardized, but in fact there is broad variability between logging companies. Each curve also has units and a description field in the LAS file.

How do you figure out what curves are contained in such a collection of LAS files? Well, one approach is to convert each LAS to an excel file, but that gets messy and brings in the data as well as the curve information. The Agile Geoscience folks have discussed this topic and python in detail; here we have a less ambitious goal that allows a simple solution.

To deal with this problem, I wrote a code in python that I call - here it is

import lasio
import os

n = 1
for file in os.listdir("."):
    if file.endswith(".las") or file.endswith(".LAS"):
        var =
        print "\n las"+str(n), file
        print var.well
        print var.other
        for curve in var.curves:
            print("%s\t[%s]\t%s\t%s" % (curve.mnemonic, curve.unit, curve.value, curve.descr))
        n += 1

This little code is magic, mainly because it draws on lasio, a library for reading and writing LAS files. Stepping through the code line-by-line goes like this
  1. import the lasio library
  2. import the operating system library
  3. set counter n to 1
  4. for each file in the current directory...
  5. if the file name ends in las or LAS...
  6. read the file
  7. print 'las' with n appended, then the file name
  8. print the 'well' data section of the file
  9. print the 'other' data section of the file
  10. for each curve in the file print mnemonic, units, value and description
  11. increase n by 1 and continue till all files are read
My early versions of this code required the file names to be pasted into the code, which was ok for a few LAS files. But that got very clunky with more than four files.

I work on a mac which has a command line terminal and python build in. Python has lots of default libraries, but others are always needed. I install them with pip, but first you must install pip itself

sudo easy_install pip

Then this command line installs lasio

sudo pip install lasio

Now you are ready to go. The code above can be copied into a text editor and saved as (or another name you like better). Now go to a directory full of LAS files from a well and copy into that folder. On the mac, open a terminal window and cd to the LAS directory, then type


or, if you want to capture the output into a file

python > info_jones1a.txt

The result will look something like the listing at the end of this post. The beauty of this little code fragment is the amount of information it makes readily available in complex well situations with many, or redundant, or overlapping LAS files. 

Hope you enjoy the code and it helps you organize a few tough LAS situations.

las1 Jones_1A_run1.las

Mnemonic  Unit  Value                       Description      
--------  ----  -----                       -----------      
STRT      F      437.0                      START DEPTH      
STOP      F     4000.5                      STOP DEPTH       
STEP      F     0.5                         STEP             
NULL            -999.25                     NULL VALUE       
COMP            Jones                       COMPANY          
WELL            1A                          WELL             
FLD             WILDCAT                     FIELD            
LOC             SHL: 1320' FNL & 2210' FWL  LOCATION         
CNTY                                        COUNTY           
STAT                                        STATE            
CTRY                                        COUNTRY          
SRVC                                        SERVICE COMPANY  
API                                         API NUMBER       
DATE            8-Jun-2012                  LOG DATE         
UWI                                         UNIQUE WELL ID   

C1 [IN] Caliper 1 {F11.4}
C2 [IN] Caliper 2 {F11.4}
CDF [LBF] Calibrated Downhole Force {F11.4}
CTEM [DEGF] Cartridge Temperature {F11.4}
DCAL [IN] Differential Caliper {F11.4}
DEVI [DEG] Hole Deviation {F11.4}
GR_EDTC [GAPI] EDTC Gamma Ray {F11.4}
GTEM [DEGF] Generalized Borehole Temperature {F11.4}
HAZI [DEG] Hole Azimuth {F11.4}
HAZIM [DEG] Hole Azimuth {F11.4}
P1AZ [DEG] Pad 1 Azimuth {F11.4}
RB [DEG] Relative Bearing {F11.4}
SDEV [DEG] Sonde Deviation {F11.4}
SDEVM [DEG] Sonde Deviation {F11.4}
STIT [F] Stuck Tool Indicator, Total {F11.4}
TENS [LBF] Cable Tension {F11.4}

las2 Jones_1A_run2.las

Mnemonic  Unit  Value                       Description      
--------  ----  -----                       -----------      
STRT      F      437.0                      START DEPTH      
STOP      F     4000.5                      STOP DEPTH       
STEP      F     0.5                         STEP             
NULL            -999.25                     NULL VALUE       
COMP            Jones                       COMPANY          
WELL            1A                          WELL             
FLD             WILDCAT                     FIELD            
LOC             SHL: 1320' FNL & 2210' FWL  LOCATION         
CNTY                                        COUNTY           
STAT                                        STATE            
CTRY                                        COUNTRY          
SRVC                                        SERVICE COMPANY  
API                                         API NUMBER       
DATE            8-Jun-2012                  LOG DATE         
UWI                                         UNIQUE WELL ID   

AF10 [OHMM] Array Induction Four Foot Resistivity A10 {F11.4}
AF20 [OHMM] Array Induction Four Foot Resistivity A20 {F11.4}
AF30 [OHMM] Array Induction Four Foot Resistivity A30 {F11.4}
AF60 [OHMM] Array Induction Four Foot Resistivity A60 {F11.4}
AF90 [OHMM] Array Induction Four Foot Resistivity A90 {F11.4}
AO10 [OHMM] Array Induction One Foot Resistivity A10 {F11.4}
AO20 [OHMM] Array Induction One Foot Resistivity A20 {F11.4}
AO30 [OHMM] Array Induction One Foot Resistivity A30 {F11.4}
AO60 [OHMM] Array Induction One Foot Resistivity A60 {F11.4}
AO90 [OHMM] Array Induction One Foot Resistivity A90 {F11.4}
AT10 [OHMM] Array Induction Two Foot Resistivity A10 {F11.4}
AT20 [OHMM] Array Induction Two Foot Resistivity A20 {F11.4}
AT30 [OHMM] Array Induction Two Foot Resistivity A30 {F11.4}
AT60 [OHMM] Array Induction Two Foot Resistivity A60 {F11.4}
AT90 [OHMM] Array Induction Two Foot Resistivity A90 {F11.4}
AHFCO60 [MM/M] Array Induction Four Foot Conductivity A60 {F11.4}
AHMF [OHMM] Array Induction Mud Resistivity Fully Calibrated {F11.4}
AHORT [OHMM] Array Induction One Foot Rt {F11.4}
AHORX [OHMM] Array Induction One Foot Rxo {F11.4}
AHSCA [MV] Array Induction SPA Calibrated {F11.4}
AHTCO10 [MM/M] Array Induction Two Foot Conductivity A10 {F11.4}
AHTCO20 [MM/M] Array Induction Two Foot Conductivity A20 {F11.4}
AHTCO30 [MM/M] Array Induction Two Foot Conductivity A30 {F11.4}
AHTCO60 [MM/M] Array Induction Two Foot Conductivity A60 {F11.4}
AHTCO90 [MM/M] Array Induction Two Foot Conductivity A90 {F11.4}
AHTD1 [IN] Array Induction Two Foot Inner Diameter of Invasion(D1) 
AHTD2 [IN] Array Induction Two Foot Outer Diameter of Invasion (D2) 
AHTRT [OHMM] Array Induction Two Foot Rt {F11.4}
AHTRX [OHMM] Array Induction Two Foot Rxo {F11.4}
CDF [LBF] Calibrated Downhole Force {F11.4}
CFTC [HZ] Corrected Far Thermal Counting Rate {F11.4}
CNTC [HZ] Corrected Near Thermal Counting Rate {F11.4}
CTEM [DEGF] Cartridge Temperature {F11.4}
DCAL [IN] Differential Caliper {F11.4}
DNPH [CFCF] Delta Thermal Neutron Porosity {F11.4}
DPHZ [CFCF] HRDD Standard Resolution Density Porosity {F11.4}
DSOZ [IN] HRDD Standard Resolution Density Standoff {F11.4}
ECGR [GAPI] Environmentally Corrected Gamma-Ray {F11.4}
GDEV [DEG] HGNS Deviation {F11.4}
GR [GAPI] Gamma-Ray {F11.4}
GTEM [DEGF] Generalized Borehole Temperature {F11.4}
HCAL [IN] HRCC Cal. Caliper {F11.4}
HDRA [G/C3] HRDD Density Correction {F11.4}
HDRB [G/C3] HRDD Backscatter Delta Rho {F11.4}
HGR [GAPI] HiRes Gamma-Ray {F11.4}
HMIN [OHMM] MCFL Micro Inverse Resistivity {F11.4}
HMNO [OHMM] MCFL Micro Normal Resistivity {F11.4}
HNPO [CFCF] HiRes Enhanced Thermal Neutron Porosity {F11.4}
HPRA [] HRDD Photoelectric Factor Correction {F11.4}
HTNP [CFCF] HiRes Thermal Neutron Porosity {F11.4}
NPHI [CFCF] Thermal Neutron Porosity (Ratio Method) {F11.4}
NPOR [CFCF] Enhanced Thermal Neutron Porosity {F11.4}
PEFZ [] HRDD Standard Resolution Formation Photoelectric Factor 
PXND_HILT[CFCF] HILT Porosity CrossPlot {F11.4}
RHOZ [G/C3] HRDD Standard Resolution Formation Density {F11.4}
RSOZ [IN] MCFL Standard Resolution Resistivity Standoff {F11.4}
RWA_HILT[OHMM] HILT Apparent water resistivity {F11.4}
RXO8 [OHMM] MCFL High Resolution Invaded Zone Resistivity {F11.4}
RXOZ [OHMM] MCFL Standard Resolution Invaded Zone Resistivity {F11.4}
SP [MV] Spontaneous Potential {F11.4}
SPAR [MV] SP Armor Return {F11.4}
STIT [F] Stuck Tool Indicator, Total {F11.4}
TENS [LBF] Cable Tension {F11.4}
TNPH [CFCF] Thermal Neutron Porosity {F11.4}

Tuesday, April 25, 2017

Carbonate Essentials Webinar - Day 1

Had a very good session with several questions. We had 28 attendees, but some of them were group logins. Total registration is 50. Tomorrow is part 2!

Screen shot of the webinar in progress. Attendee names redacted to protect the innocent. 
You can always depend on Mike Graul for a good dose of humor.

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.

Monday, January 30, 2017

KFSM Earthquake Insurance Interview

Just finished an interview with Charlie Hanna of KFSM channel 5 in Ft. Smith, AR. Topic was earthquake insurance. It will air sometime in the next couple of days. Below is my whiteboard drawing made for the interview.

Explanation of the figure: In N and NE Oklahoma horizontal wells are drilled into unconventional reservoirs, typically the Mississippian Lime. These wells produce vast quantities of formation (salt) water mixed with a small percentage of oil. The oil is separated and the water needs to be disposed. For this purpose a deeper nearby vertical salt water disposal well pumps the water down into the 500 million year old Arbuckle formation that sits on granite basement rocks. The basement and Arbuckle often have ancient inactive faults from continental collisions that formed the North American continent. Salt water injection that is too rapid, at too high a pressure or from wells too close together can reactive the faults and be felt as earthquakes. The earthquake waves radiate out from the fault and, sometimes, can be felt 100 miles away in Fayetteville, AR where I live. The hilly terrain of NW Arkansas scatters and focusses the earthquake waves so that the felt effect varies widely from one neighborhood to the next.

Sketch of NE OK earthquake as felt in Fayetteville, AR

Sunday, January 29, 2017

Seismic interpretation bootcamp

Yesterday I taught a 1-day seismic interpretation bootcamp short course for the Tulsa Geological Society. Everyone worked on their own laptop running OpendTect software and the Netherland's offshore F3 Demo data from the dGB website.

Kudos to Shane Matson for organizing and the 25 who attended. Lots of familiar faces in the room, including Dick Banks, Jamie Woolsey, and Eric Gross. The OSU-Tulsa campus as an excellent venue. It was a great day!

Bootcampers using OpendTect in Tulsa.
Bootcamp advertisement

Thursday, January 26, 2017

Lyell and Darwin

As some Seismos readers know, a passion of mine is collecting rare books. I have been at it since the late 1970s and now have about 120 titles that one could term rare.

The photo above is a very nice part of the collection, completed only this week with he addition of Charles Darwin's The Origin of Species. The 1876 version of the 6th edition of Origin is the basis for all modern reprinting of his most famous book. It was the last edition edited by Darwin himself before his death in 1882. Only 1250 copies were printed, making 1876 the smallest print run except the 1859 first edition of Origin that also numbered 1250.

Charles Lyell was an older contemporary of Darwin and a close colleague. They are together again on my bookshelf.

From left-to-right in the photo, we have:

Lyell, 1834, Principles of Geology, 4 vol (3rd ed)
Lyell, 1885, Students Elements of Geology (4th ed)
Lyell, 1865, Elements of Geology (6th ed)
Lyell, 1863, Geological Evidences of the Antiquity of Man (1st ed)
Darwin, 1876, The Origin of Species (6th ed)
Darwin, 1871, The Descent of Man (1st ed, 2nd issue)

A rather shabby first edition of Origin sold at Sotheby's in December 2016 for $102K.

Wednesday, January 11, 2017

Seismic Sherlock

When one works with students teaching seismic interpretation it can be a sinking feeling. All the basic elements get covered in the right order. Faults are picked, horizons are tracked, amplitude is extracted, attributes are applied. High tech, lots of workflow, parameters and button pushing. Then a question comes up and the student answer makes it clear something is missing. Perhaps after picking 20 faults the student has no clue of the structural style. Or noise is closely interpreted to have geologic meaning. Of course every mistaken concept is discussed and corrected, but it leaves the teaching like someone standing before a crashing surf having explained a few pebbles.

I've been thinking there may be a better way to jump start the learning process about seismic interpretation. No cubes, or N-level software applications, not even a colored pencil. Perhaps we need to think of a seismic section like a dead body in the morgue and the teacher is like one of those Victorian doctors who lectures to a classroom in the round as his assistant dissects the corpse, or  Sherlock explaining the geological implications and geophysical limitations of a seismic section...

Here is the patient (click to enlarge).
Data fade at the top of the section is due to shallow event muting and subsequent loss of CMP fold in the gather that has been summed to make each poststack trace we are interpreting. The mute cut is uneven, as evidenced by the irregular depth (that is, time) of the fade. Any interpretation in the mute zone is extremely hazardous, amplitudes are unreliable, structure is uncertain because missing offsets in the gathers make velocity analysis highly suspect. Only fools and lunatics interpret in the mute zone. Note the strong event at 250 ms on the left end of the line is undisturbed, no obvious lateral time or amplitude disruptions. At the same depth on the right side are lateral amplitude bursts with a 'sting of pearls' look. This is characteristic of spatial aliasing that arrises from poststack traces that are too far apart for the frequency and velocity at that location. The spatial aliasing effect in this area of the data is doubtless due to shallow (and therefore low) velocities, the effect largely subsiding by 400 ms.

Revisiting the strong 250 ms event, is shows as a peak-over-trough waveform, assuming black is positive amplitude and polarity is SEG normal. If it is further assumed that the seismic processor has done the job correctly and the wavelet is zero phase, it follows that the 250 ms event must be a seismic thin bed. A thin bed with zero phase wavelet, it is to be recalled, has the appearance of a 90 degree phase shifted wavelet, or a positive-over-negative as observed here. To quantify a bit with 10000 ft/s velocity and 55 Hz dominant frequency in this zone, the vertical resolution is about 45 ft, so the waveform of this event indicates it is less than 45 ft thick. Yet this event is very bright. With some basic investigation, one finds the surface geology is Pennsylvanian sandstone and shale, in short clastics. The possible lithologies that could give a strong reflection for a thin bed in an otherwise clastic section are anhydrite or carbonate (limestone or dolomite). Most likely a limestone, but well verification is needed.

At this point it is useful to jump to the deeper section and work back up the 250 ms event. This line is from the US mid-continent so we can safely assume the basement rocks are igneous, specifically granite. Not all granites are created equal with respect to velocity and density, but a fair guess is 21000 ft/s (6400 m/s) for P-wave velocity and 2.65 g/cc for density. Granitic basement is often layered and this is the case here. An autocorrelation test would quickly revel absence of multiples, so the basement layering is real. You can on the fact that the top of the basement is an unconformity surface.

The basement reflection event here is at 800 ms on the left, represented by a strong weaker trailing peak. If we are to have a negative amplitude event then what overlies the basement must have a higher impedance than granite. A tall order, but dolomite has just this property. Also, lateral variability of the basement reflection amplitude makes one suspect there is a zone of weathered granite or conglomerate at the basement contact. Notice in the center of the line there is an irregular group of hills representing erosional remnants or later uplift. In any case, we see the basement reflector fade to nearly nothing in this region, indicating unaltered granite overlain by tight dolomite. The impedance of these two rock types being almost equal, a case of matched impedance making the granite interface almost invisible.

The peak (dark) event on the LHS at 680 ms is likely a shale contact with tight dolomite below. Lateral variability of the event amplitude is probably related to thickness changes in the shale. On the left side, the disruption of event 680 may not be a geological reality. Immediately above we see a deep sinkhole indicating a limestone unconformity (limestone because it is far more soluble than dolomite); this unconformity extends across the entire line. If the sinkhole is filled with low-velocity rubble it will have the effect of a velocity sag. This is a false structure. On the other hand, if it were a velocity sag all deeper events should sag equally. But they do not, suggesting it is a real dissolution feature at the limestone unconformity and slowly healing with depth. How deep is this sinkhole? It looks to be about 30 ms

And speaking of the limestone unconformity at 610 ms (on the left), it is a doozy. In fact, the Missippippian-Pennsylvanian boundary that occurs in most places around the world. Typically, the overlying Penn rocks are clastics (shale and sandstone)Across the line the amplitude is not only irregular, but often changes polarity. It is also structurally irregular, pocked with hills and valleys. Here we have classic paleokarst on a limestone unconformity surface. Karst in the modern world is a landscape characterized by subsurface drainage through fracture networks and collapse features. Dissolution occurs due to atmospheric and soil biota CO2 forming acidic ground water that dissolves the limestone pretty effectively. The same process that rounds out the crisp lettering in old limestone cemetery markers. Once buried, the altered limestone near the unconformity  translates into highly variable acoustic impedance. For example, on the far right the unconformity pops as a strong negative (white) event, surely an area of deeply weathered limestone with associated high porosity and low acoustic impedance. But in the center of the line, we see the unconformity as a tight positive (black) event. Here we suspect overlying sandstone and shale in direct contact with hard, unaltered limestone.