Author: kvan637 (page 1 of 8)

Congratulations, Dr. Jami Johnson!

It’s the end of an era… Jami Johnson successfully completed her PhD thesis in photoacoustic imaging.  She is now off to The University Pierre and Marie Curie, in Paris (France) to pursue postdoctoral research in medical imaging. We wish Jami all the best, and we will follow her in what will undoubtedly be a very exciting career. The PALs want to thank Jami for all her efforts and the wonderful memories she leaves behind.


Postgraduate opportunities in our lab

Below are short descriptions of research projects available for postgraduate students in our group at the University of Auckland. If you are interested in any of these, please contact us for technical questions. For questions around enrollment at the University of Auckland, please go here. Scholarships may be available for the best applications through the University of Auckland, or the Dodd Walls Centre.

ResonanT Ultrasound Spectroscopy

Resonant ultrasonic spectroscopy fills an important gap between our ultrasonic and seismic research. Together with the PORO lab, we have rock samples to complement laser ultrasound and a host of other petrophysical data with new RUS results.

Laser ultrasonic rock physics under high pressure and/or temperature

The PAL and PORO lab join forces by combining our respective strengths in laser ultrasound and rock physics to improve data quality and quantity. In this project, we are building the expertise to do laser ultrasound in a pressure vessel with optical windows. Source/receiver locations are varied under computer control with an arduino-controlled servo rotational stage.

unraveling the mysteries of the Auckland volcanic field

Part of an active volcanic field, questions surround the nature of Auckland’s past, present and future. Using a suite of seismic techniques that range from ambient seismic noise tomography, to receiver functions and body wave tomography, we aim to build a representative model that helps us explain the geodynamics of the Auckland Volcanic Field.

quality control of timber and fruit products with laser ultrasound

Laser ultrasound can be applied to products of interest to a wide community. Current methods of testing the quality of fruit and timber, for example, can often be described by one or more of the following terms: sparse, contacting, expensive, and often destructive. In this project, we aim to explore the opportunities for laser ultrasound in estimating the quality of  fruit and timber in a non-contacting and non-destructive matter.

Our first glimpses into the Auckland Volcanic Field

Today, the New Zealand Journal of Geology and Geophysics published Josiah‘s first results using ambient seismic noise to peek deep into the lithosphere under the Auckland Volcanic Field (AVF). This noise, created by the oceans surrounding the AVF, suggests the lithosphere under Auckland is made up of oceanic and continental components. Further studies aimed at increasing the resolution of our results may help explain why this active volcanic field is where it is.

Laser ultrasonics on ice cores

Many years ago, Andrei Kurbatov asked us if we could do laser ultrasound on ice cores, during half time of an international football match we were watching in Golden, Colorado. Dylan Mikesell started these measurements in a home-made refridgerator in the Physical Acoustics Labd as a research project toward his undergraduate degree at Colorado School of Mines. We have come a very long way since, and today — with the additional help of Thomas Otheim and Hans Peter Marshall, our first paper on laser ultrasound on ice cores from Greenland and the Antarctic appeared in the journal Geosciences.

A new paper on Marchenko-equation medical imaging

Our collaboration with TU Delft has led to a new publication in the Journal of the Acoustical Society of America. Led by Joost van der Neut and Jami Johnson, this article adds the Marchenko equation to our tool set to image ultrasonic and photoacoustic data for medical applications. Perfect timing, Jami, for your PhD thesis defense tomorrow!

Lab photos

The orientations of the stations in the Auckland seismic network

The stations of the seismic network that monitors the Auckland Volcanic Field are a mix of surface and borehole seismographs. Our aim here is to figure out the orientation of the three orthogonal components of each seismometer.

The vertical component

The vertical components of the seismic network are tested with some earthquakes with an epicentre that is only a ~100 km from the AVF. These are very deep earthquakes, assuring that high-frequency signal impinges from below on the network.

The figure below shows that seismic data from these events generally results in an arrival with a negative maximum, except stations KBAZ and WTAZ:


From this example and seismic data from other earthquakes, we conclude the polarity on the vertical component is reversed for stations KBAZ and WTAZ.

The horizontal components

The majority of the seismic stations in the AVF are in boreholes between 50 to 400 m deep. Particularly for these stations,  the orientation of the horizontal components of the each station is unknown. Luckily, there are several techniques to estimate the orientations of the orthogonal horizontal components.

Ambient Noise Orientation

Stachnik et al. (2012) used a earthquake-sourced Rayleigh wave polarisation method to estimated the orientation of ocean bottom seismometers. Zha et al. (2013) used ambient seismic noise, applying a polarisation technique to the Rayleigh wave signal in CCFs to estimate the orientation of ocean-bottom-seismometers. Here we use ambient noise orientation, based on the method used by Zha et al. (2013), but we use it for borehole seismometers. We performed ambient noise orientation on all 3-component seismometers in the network, using the surface stations as a control group. The  method involves three basic stages:

This method involves three stages:

  1.  Computing the three-component cross-correlation functions for all station pairs from noise.
  2. Measurement of sensor orientation through an optimization process.
  3. Data selection and statistical analysis to obtain final orientation angles.
three-component cross-correlation functions

Cross-correlation functions (CCFs) can be computed for any pairs of seismic sensors. Computing all the possible cross-correlation functions for a pair of three-component seismometers would produce nine functions. For this technique, we use some of the CCFs that potentially contain Rayleigh wave signal. Rayleigh wave motion occurs in the vertical plane so we expect virtual Rayleigh wave signal to emerge in the vertical-vertical(C_zz), and the radial-vertical (C_rz) CCFs. For seismometers with unknown orientation, the C_rz data is distributed across the horizontal-verticals (C_1z and C_2z) cross-correlation functions.


For each pair of seismometers, the symmetric CCFs are bandpass filtered to the frequency band where Rayleigh wave with the greatest amplitude are expected to emerge (In our case 0.05–0.5Hz). We then perform a grid-search through 360-degrees, each performing as an estimate of the angle that would rotate C_1z and C_2z into the C_rz and C_tz . For each estimate of C_rz,  the zero-lag cross-correlation and coherence between the C_rz and 90-degrees phase-shifted C_zz is computed. The maximum cross-correlation marks the best estimate of C_rz. The coherence at maximum cross-correlation is used in the data selection. The back azimuth can then be calculated because the coordinates of both stations are known.


The figure below this process is illustrated for the AWAZ-ABAZ pair. On the left, there is C_zz, C_1z, and C_2z, and the phase shifted C_zz in the dashed red. In the middle panel , the correlation and coherence are plotted by angle. On the right, there is C_zz, and the C_rz, and C_tz based on the angle at maximum correlation (260 degrees in this case), and the phase shifted C_zz in the dashed red.

Data Selection and Statistical Analysis

Estimates with a coherence below 0.5 at maximum correlation tend to be scattered to a greater degree. These coherences were discarded. We then take the mean of the remaining estimates as our final estimate of the orientation. In the figure below, the estimates for the AWAZ station are presented in a polar plot. The radial axis represents the number of estimates in each 1-degree bin. Red represents estimates with coherences below 0.5 that are removed.


Stachnik, J.C., Sheehan, A.F., Zietlow, D.W., Yang, Z., Collins, J., & Ferris, A., 2012. Determination of New Zealand ocean bottom seismometer orientation via Rayleigh-wave polarization. Seismological Research Letters, 83(4).

Zha Y, Webb SC, & Menke W. Determining the orientations of ocean bottom seismometers using ambient noise correlation. Geophysical Research Letters. 2013 Jul 28;40(14).

Ambient seismic noise correlations for the AVF

Here is an example of the correlation of 5 days of noise between stations WTAZ and MBAZ. The positive times are seismic surface waves going from WTAZ to MBAZ, whereas the negative times would capture seismic waves going in the opposite direction. The reason you see a distinct wave on the positive side of this plot, but not on the negative side is that the predominant noise (from these 5 days) is from the West. It also appears that the earlier part of the wave is lower in frequency than the higher parts. This is called dispersion and stems from the fact that longer-wavelength surface waves travel through deeper parts of the Earth, which are generally faster. This information can be used to extract the seismic velocity structure in the ground between these two stations.


The cross-correlation of 5 days of seismic noise between WTAZ and MBAZ, representing surface seismic waves travelling between these stations. provides a suite of tools in Python to calculate the wave fields of these “virtual earthquakes” between each combination of station pairs of the AK network. We repeat this procedure each day, and look for small differences in the calculated waveforms to inform ourselves of possible changes in the subsurface. For example, on Mount Merapi, Indonesia, changes in the waveforms correlate with changes in precipitation saturating the near-surface. However, in principle magmatic changes could be detected this way. Areas of active research include how to do the calculations so that we can exclude causes for changes in the waveforms not related to the Earth. For example, the distribution of noise sources may vary from day to day, affecting the calculated “virtual earthquakes.”


First estimates of changes (or the lack thereof) in the seismic wave speed in the AVF. Disclaimer: these results are raw data that need further study!


This blog will describe the running of FMTT on the AVF.


The installation was done by following the instructions obtained from

In fact, all the necessary files are attainable in an all in one download as a compressed Unix gzipped and tarred file specified in that link. In which it includes the instruction manual, a research paper of its application, and a couple of examples.

In a nutshell, to install all we need to do is to run the shell script compileall after the compressed file has been unpacked. But for things to run smoothly, prior to running the shell script what we did was:

  • Edit line 22 of  compileall and line 5 of source/aktimes/Makefile into the appropriate Fortran compiler. In my case this is gfortran.
  • Still in compileall, add ./ before remodl and setbrn in line 82 and 83 respectively so that these can be found during compilation.
  • Edit line and line 19 and line 181 of syntht.f90 source file. Delete (KIND=i10) in both lines to avoid mismatch of function error.

Compilation complete

FMTT itself consist of several individual executable programs. Each one requires an input file which name is identical to the executable except for the .in extension. In these input files we can specify and modify the parameters so that the desirable inversion routine is performed. More details on each program can be found inside the instruction manual.

An executable worth to note is tomo3dt, which perform the whole tomographic inversion as opposed to run the necessary executable programs one after the other. Another notable mention is gmtslicet. This program is used to help plot and visualize the inversion result as depth, E-W, and N-S slices. The complete plotting procedure is based on GMT scripts (GMT stands for Generic Mapping Tools).

Running the Example

Two examples were provided in the compressed file. Prior to do any of the examples, ak135.hed and ak135.tbl need to be copied into the respective example directory. These files can be found in source/aktimes after compilation. Once the files are copied, simply execute tomo3dt from the directory. It should returns the following lines:

gugi@sc-9020-hs:~/AVF/fmtt_v1.0/example1$ tomo3dt
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!
Program fm3dt has finished successfully!

The number of times the line “Program fm3dt has finished successfully!” is repeated depends on the number of iteration specified in the, the input file for tomo3dt. This shows that we have managed to perform a successful tomographic inversion.


Plotting the Examples

After tomo3dt has been successfully run and the rounding issue in gridi.vtx and gridc.vtx has been adjusted, the results are visualized using gmtslicet and GMT. Go to the subdirectory /gmtplot and execute gmtslicet. Prior to executing gmtslicet, the files gridi.vtx and gridc.vtx might need to be adjusted. In the top lines of both files, the numbers round off differently to each other which could result in an error while running gmtslicet. Equalizing these values will solve the error. Nick Rawlinson kindly provides us with an updated gmtslicet.f90 source file. It meant to tackle the rounding problem. It works in some computer but somehow it could be machine specific.

The execution of gmtslicet will produce three bound .gmt files (among others) for depth, E-W, and N-S strike. These are used by the plotting scripts plotd, plotew, and plotns. Prior to running these scripts, some modification were made. #!/bin/csh was added on the top of each script file. Then, depending on how GMT is called in your computer, gmt may need to be added before each GMT command namely gmtset, xyz2grd, grdimage, psscale, psxy, pscoast. Output file can also be modified.

The results of running the plotting scripts are of .ps file type which can be viewed. Depth or E-W or N-S slice is depicted by their respective script file.

Jonathan wins the S.J. Hastie Scholarship

Our youngest PAL, Jonathan Simpson, has been working this summer on performing laser ultrasonics on rocks in our pressure vessel. This semester, Jonathan will embark on a BScHons thesis in our group. Today, he was awarded the S. J. Hastie Scholarship by the Geoscience Society of New Zealand. We congratulate Jonathan, and wish him good luck in his BScHons project on leaky modes in the Hikurangi Trench!

Jonathan receiving the Hastie Scholarship from Jennifer Eccles, vice-president of the GSNZ.

Skip to toolbar