Full 3D Monopole & Dipole Open Hole Sonic Logging in Fast, Slow and Anisotropic Formation


Acoustic wave propagation in multilayered structures is used in different industries for numerous measurements, evaluation and imaging applications. Technology based on this method is used in the Oil and Gas industry to characterise underground reservoirs during their entire life cycle; starting from exploration to their abandonment. This technology is referred to as Borehole Acoustic Measurement (Logging) and is further classified in to Sonic and Ultrasonic Logging

Sonic Logging Principle:

Borehole Sonic Logging is based on the principle of measuring the travel time of Compressional (P), Shear (S) and Stoneley waves propagating along a borehole with different radial depths into the formation. The operating frequency range for these measurements is typically from 0.5 to 20 kHz.  This helps to probe several meters into the formation. In this tutorial we will show you how to set up a sonic logging model and run it on the cloud.

The Model:


The model we will show you how to build today will be a fluid-filled borehole surrounded by a fast or slow formation. A formation is characterised as fast or slow depending on whether its shear wave velocity is slower or faster than the compressional wave velocity in the borehole fluid. 

The sonic measurement tool or in this case the pressure wave excitation will be located at the centre of borehole, the receivers will also be located along the centre of the borehole the first receiver will be placed 1.219m from the transmitter this will be referred to as the Tx-Rx spacing and the receivers are spaced periodically ever 0.152m the Rx-Rx spacing.


The model will be excited with both a monopole and a dipole source. Due to the size of this model as well we will need to utilise OnScales MPI capabilities on the cloud to run these models in a reasonable time. We will show the current workflow needed to setup an MPI model!

Download: 3D Sonic Logging Model

Step 1: Define variables 

We need to define the following variables:

  • Borehole, tool & formation radii 
  • Height of borehole and sonic tool 
  • Tx-Rx and Rx-Rx spacing 
  • Mesh settings

Using the height of the sonic tool and the spacing variables you can calculate the number of axial receivers. The following inputs that we use will yield a total of 24 receivers.

Below is the code - it can be copy and pasted directly into a new Analyst script.

c ****************************
c ****************************

c All inputs should be in this section
symb r1 = 0.065 /* tool inner radius
symb r2 = 0.100 /* tool outer radius
symb r3 = 0.13 /* open hole radius
symb rform = 6 * $r1 /* formation radius estimate
symb h1 = 6 /* height of borehole
symb h2 = 5 /* height of sonic logging tool
symb hs = 0.5 * ( $h1 - $h2 ) /* start of the tool (equidistant from top & bottom)
symb hr = 0.285 - 0.152 /* location of last receiver on the tool
symb tr = 1.219 /* transmitter to 1st receiver spacing
symb rr = 0.152 /* receiver to receiver spacing

symb nr = ( $h2 - $tr ) / $rr /* no. of axial receivers
symb ncopy = $nr - 1

c read material properties
symb freq = 5e3
symb freqdamp = $freq

c Mesh Settings (element sizes)
symb boxx = 3.25e-3 /* element size in X & Y
symb boxz = 0.01 /* element size in Z

 Step 2: Keypoints and Grid 

It is beneficial to set up the grid using keypoints to allow for rapid parameterisation and specific control of the model. We won't go into great detail in this tutorial about keypoints but check out this article for more information on keypoints and why you should be using them.

The symbol code has a number of #key options that allow for a quick method of generating XYZ and IJK model keypoints these commands also help to parameterise the model grid generation to allow for quick changes to the model geometry to be made. We will use the #key commands #keyindx and #keycord you can find more information on these commands in this article. We will also make use of the #key coommand #keycopy which allows for quick replication of already defined keypoints we will use this to create keypoints for each of the receivers along the borehole axis.

This is exactly the same as the 2D model except we have added a 3rd dimension.

This code can be copy and pasted directly into the analyst script 

c Keypoints
/* #keycordt allows XYZ keypoints to be generated from absolute dimensions instead of the distance between keypoints like #keycord
symb #keycordt x 1 -$rform -$r3 -$r2 -$r1 0. $r1 $r2 $r3 $rform
symb #get { idx } rootmax x /* store interger suffix of last keypoint
symb #keyindx i 1 $idx 1 $boxx 1

symb #keycordt y 1 0. $r1 $r2 $r3 $rform
symb #get { jdx } rootmax y /* store interger suffix of last keypoint
symb #keyindx j 1 $jdx 1 $boxx 1

symb #keycord z 1 0 $hs $hr $rr /* bottom of borehole to first receiver
symb #keycopy z 3 4 $ncopy /* copy keypoints 3 & 4 out ncopy times
symb #get { krend } rootmax z /* store interger suffix of last keypoint
symb #keycord z $krend $z$krend $tr $hs /* finish keypoints in y starting from end of keycopy
symb #get { kdx } rootmax z /* store interger suffix of last keypoint
symb #keyindx k 1 $kdx 1 $boxy 1

/* integer suffix to another variable
symb indgrd = $idx
symb jndgrd = $jdx
symb kndgrd = $kdx

/* define the number of nodes in the computational grid, the number of dimensions for the model, and the appropriate constraint relations.
grid $i$idx $j$jdx $k$kdx 

/* To define nodal coordinates for the model.
geom keypnt $idx $jdx $kdx

Step 3: MPI Settings 

When we hear people talking about MPI, in most cases, this is a implementation associated with an application using MPI. It is widely used by most legacy simulation packages and aimed at those lucky enough to have High Performance Computing (HPC) hardware within their company to run their MPI application. Luckily, we can utilise the thousands of HPCs on the cloud. We are going to show you the necessary code to run an MPI model.

This code can be copy and pasted directly into the analyst script. For this code to work it needs to be inserted after the GRID command and before the GEOM command 

symb nptot = 500    /* total number of parts to partition the model into. parts = cores
symb npi = 10 /* total number of partitions in X
symb npj = $nptot / $npi /* calculate total number of parts in Y
symb xpi = $npi
symb xpj = $npj
symb xistep = $indgrd / $xpi
symb xjstep = $jndgrd / $xpj

max $nptot
symb ifr = 1
/* Loop through parts and partition model
do loopI I 1 $npi 1
symb jfr = 1
do loopJ J 1 $npj 1
/* Calculate part number
symb ip = $J + ( $I - 1 ) * $npj
/* Calculate and j indices
symb xito = $I * $xistep
symb ito = $xito
symb xjto = $J * $xjstep
symb jto = $xjto
/* Check for last part
if ( $I eq $npi ) then
symb ito = $indgrd
if ( $J eq $npj ) then
symb jto = $jndgrd
/* Define part
sdiv $ip $ifr $ito $jfr $jto $k1 $kndgrd
symb jfr = $jto
end$ loopJ
symb ifr = $ito
end$ loopI

Step 4: Defining Materials and assigning them to the grid 

OnScale has a material database with a number of common materials readily available for users to use however for this tutorial we will need to manually define the fast and slow formations using the material properties shown in the image above. First though we will use the database to save water to a file and we will edit that file. Save the file to the same working directory as the flxinp file 

  1. Expand fluid section 
  2. Double click watr to add this material to the Project Materials 
  3. Right click watr and select Copy - Call the new material watr2  
  4. Repeat steps 1 & to add Steel from the Metal section
  5. Click Save to File -  save the file as sonic_logging_mat.prjmat 

We have created a copy of watr to act as an interface to define a pressure load when it gets to that point

Once you have saved the file open it in the OnScale application add the following lines of code starting on line 58 of that file. We have also included the anisotropic material properties for this model 

 /* fast formation
wvsp on
type elas
prop form_fast 2350 3658 2032 0.01

/* slow formation
wvsp on
type elas
prop form_slow 2054 1890 508 0.01

/* anisotropic material properties
wvsp off
type lean
prop form_aniso 2350
23.87 15.33 9.79 0 0 0 23.87
9.79 0 0 0 15.33 0 0
0 2.77 0 0 2.77 0 4.27

Now that materials have been defined, we can now assign those materials to elements of our grid to do this we use the SITE command. This is where keypoints come in handy because now we can define regions between start and end keypoints in a given direction to assign material to all the elements encapsulated by that window.

For this model instead of assigning a cylindrical region of nodes to assign the fast formation material to we will assign the formation material to the full grid and assign water and steel to the appropriate cylindrical regions. Why? Because it is much easier to apply the appropriate boundary conditions to the model if the materials lie on the extents.

The model we are setting up here will take advantage of symmetry conditions on the YMIN boundaries

The code can be copy and pasted directly into the analyst script you are working with 

c **************************************************
c **************************************************
symb #read sonic_logging_mat.prjmat

regn form_slow /* All materials must be assigned a material initalise them all with void
cyln watr stnd z z1 $z29 $x5 $x5 $x8 $x8 * * 0 360
cyln steel stnd z $z1 $z29 $x5 $x5 $x7 $x7 $x6 $x6 0 360
regn watr2 $i5 $i5+1 $j1 $j1+1 $k28-1 $k28+1

plot matr /* plot geometry

Step 5: Boundary Conditions and Loads 

It is best to simplify the model as much as possible. If you are modelling a symmetrical object and the same conditions are applied to both sides of the model, there is no reason to model the entire structure. It is more efficient to model half (or less) of the structure and achieve the same results for a fraction of the computing time. We will apply a symmetry boundary on the minimum X extent of the grid. For more information on the topic of boundary conditions check out the following articles:

Boundary Conditions & Loading

How to Apply your Own Boundary Conditions

IMPD (Impedance) vs. ABSR (Absorbing) Boundary Conditions

We will also apply a pressure load on the interface of the materials watr and watr2 this will act as a monopole pulse we will also provide you with the text file that contained the time and amplitude values of the pulse we used. 

For this model we will excite it with a monopole pulse and a dipole pulse. The below code snippet contains both options, but you will need to comment out the ones you don't want to use.

c External boundary conditions
side xmin impd 2054 1890 508
side xmax impd 2054 1890 508
side ymin symm
side ymax impd 2054 1890 508
side zmin absr
side zmax absr

c Define time varying funciton
data hist pulse * pulse_5kHz.txt /* open the text file and save it to an array called pulse
func hist pulse /* assign array to a time history func

/* Pressure Loads (monopole pulse)
pdef pld1 func
matr mat1 watr2 out
sdef pld1 mat1 $i5 $i5+1 $j1+1 $j1+1 $k28-1 $k28+1

/* Pressure Loads (dipole pulse)
pdef pld1 func
vctr vct1 1 0 0
sdef pld1 vct1 $i5 $i5 $j1 $j1+1 $k28-2 $k28+2

Step 6: Calculations and Outputs 

The receivers situated on the borehole axis measure the time it takes for the acoustic pulse to travel in the formation in various forms while undergoing dispersion and attenuation. When the sound energy arrives at the receiver, having passed through the formation, it does so at different times in the form of different types of waves. This is because these waves travel with different velocities we want to capture these waveforms so we need to set up outputs to capture the acoustic pressure measured on each receiver 

c Request calculation of additional arrays
calc pres aprs /* open an additional APRS array - acoustic pressure

c Request time domain outputs saved to a flxhst file
hist func /* write function to a output file
symb jrx = 3 /* keypoint suffix of first receiver
do loop_rx J 1 $nr 1 /* open do loop to loop through receivers
hist aprs $i5 $i5 $j$jrx $j$jrx /* loop through receiver locations to output acoustic pressure at these points
symb jrx = $jrx + 1 /* add one to the keypoint suffix
end$ loop_rx /* close loop

Step 7: Process and Model Execution 

Before running the model, you should always issue the prcs command. In preparation for processing the time-step information, it computes the model time step for stability and opens the necessary arrays for field-data storage.

The basic information required prior to executing the model:

  • Model timestep (step)
  • Simulation runtime (simtime)
  • Simulation runtime in terms of number of timesteps (nexec)
c ***********************
c ***********************
c Time settings
time * * 0.8

c Process model

set imag mpeg /* set to generate a video
nvew 1
mirr x on

c ***********************
c ***********************
symbx ncycles = 250
symb simtime = $ncycles / $freq /* 0.01
symb nloops = 100
symb #get { step } timestep /* extract calculated timestep from prcs command
symb nexec_loop = ( $simtime / $step ) / $nloops /* partition full simulation into smaller segments for plotting

c Runtime calculations
proc exec_loop save /* name procedure exec_loop and save the code
exec $nexec_loop /* run model partially
mirr x on
plot aprs rang 0 /* plot acoustic pressure - shows wave propagation

end$ proc /* end of procedure

c run model by calling 'plot' procedure 'nloops' times
proc exec_loop $nloops

c Save symbols to file for later use
symb #get { labl } jobname /* find name of run
symb #save '$labl.symb' /* save in symb file

c End of input file

Step 8: Running the model on the cloud and downloading the results 

Before we run the simulations on the cloud, we need to enable the Advanced Cloud Job Submission option in the settings to open the settings window click mceclip0.png in the top right hand corner of the application window. Once the settings are open click Cloud Account Configuration and enable the setting 


Now that the models have been set up they can now  be ran on the cloud - To do this simply open the cloud scheduler (Run on Cloud mceclip1.png those familiar with OnScale will notice that there is no estimate option when you open the scheduler, this is because we are running an MPI simulation. In the estimators place is the Estimate Override options all you need to do as a user is request an appropriate amount of RAM for the simulation, for these simulations a RAM request of 50GB will be sufficient, so request 50GB and click Run to submit the files to the cloud and the simulation will start.

To access the storage, click mceclip1.png locate the job you have just run and download the results. 

Step 8: Post-Processing the Results 

Similar to the 2D Borehole Sonic Logging Simulations we have post-processed all of our results in MATLAB if you would like to see what the results can look like if done properly see the linked 2D tutorial there you can also download the MATLAB script that we used to process these results