Category: PALNews (page 1 of 9)

News related to the PAL

Inverting teleseismic event with adaptive stacking

getting the time residuals

So far, our travel-time information for the tomographic inversion have been obtained from manually inspecting individual traces associated with the source and receiver pairs. Here we use another method of finding the arrival time residuals, by using the adaptive stacking method. This method make use of predicted times from an earth model (eg. ak135) to get an approximate alignment of trace, before iteratively improve the alignment, thus providing us with the residuals from the aforementioned model. The code that is used to perform this is named tcas, and it requires us to specify the trace data values for the earthquake.

The example earthquake we use here is a magnitude 7.6 event in the Pacific Island region. The earthquake is recorded in 17 out of 18 stations that we consider. Eight iterations are enough to reach a convergence in tcas. Taking into account the reverse polarity in KBAZ and WTAZ, the result shows good alignment between the traces. Information about the ak135 model residuals and the error estimates can be found in one of the output files. This can then be used as the starting information for the travel-time inversion with FMTOMO.


Unfortunately our inversion tests so far have not been quite promising. Judging from the development of standard-deviation and chi^2, along with the histograms of the relative residuals, we are yet to reach the ideal outcome.

std-dev and chi^2

Nine iterations of inversion process show a decent reduction in the standard-deviation of about 2 seconds. This means a variance reduction of 83%. The histograms mostly reflect the reduction in standard-deviation with the more centered final model. However, chi^2 is rather large which indicates that the data misfit is still quite sizable.

Unsurprisingly, the slices show little implication of detailed structural features. But the general idea of a positive anomaly in the Waikato region remains. The velocity model we use is specified by the extent and the number of grid points. Latitudinally it extends from 34.5 S to 41.0 S, and from 172.5 E to 178.5 E longitudinally. Vertically it reaches the depth of 305 km and 1.5 km in elevation. Depth-wise there are 28 number of grid points and 50 grid points in both the latitudinal and longitudinal direction.

Fitting the model

The fact that chi^2 is extremely high means that we haven’t quite fit the model with the data. Which is odd considering we only have one event to fit. Our subsequent attempts of using just the outer broadband stations, as well as reducing the model extent also didn’t improve chi^2 much.

To test our procedure with this technique, we perform the simplest inversion with only one source and one station. The same event is used. We took station KUZ (coordinate 36.745 S/175.721 E and elevation of 76 m) as an example with direct travel-time residual extracted from the stacking method earlier.

The starting model is based on the ak135 1D-Earth model. There are 50 by 50 grid nodes, with 28 grid nodes to depth. The model reach to the depth of 100 km. The latitudinal extent of this model is from 35.1 S to 38.8 S, and the longitudinal extent is from 173.5 E to 176.0 E. Therefore albeit it has the same amount of grid nodes, the model is smaller than in the above section. The two interfaces, one at the top and bottom of our model, have grid points following the velocity grid nodes i.e. 50 in both North-South and East-West directions. The propagation grid is twice as compacted as the velocity and interface grid.

In the invert3d inversion parameters, interfaces are not inverted for. Subspace dimension is kept at 20. And the global damping (epsilon) and smoothing (eta) are 10 and 2 respectively. In full, excluding the associated input files the file looks like this:

0.07            c: Minimum distance between interfaces (km)
0.2             c: Minimum permitted velocity (km/s)
1               c: Remove mean from predicted teleseisms (0=no,1=yes)
1 1 1           c: Velocity inversion (0=no,1=yes),damp,smooth
0 0.5 0.01      c: Interface inversion (0=no, 1=yes),damp,smooth
0 1.0 1.0       c: Source inversion (0=no, 1=yes),damp1,damp2
20              c: Subspace dimension (max=50)
10              c: Global damping factor (epsilon)
0               c: Apply second derivative smoothing (0=no,1=yes)
2               c: Global smoothing factor (eta)
6371.0          c: Earth radius in km

With 11 iterations, the outcome of the inversion in residuals.dat shows no improvement in any iteration step. The difference between initial and final model time is also less than 0.2 of a second. It is clear that we are missing a step here.

Chapter 2: Dimensional Analysis, Buckingham Pi Theorem, and the physics of flight

Chapter 2 of A Guided Tour of Mathematical Methods for the Physical Sciences deals with dimensional analysis; a great way to test the accuracy of equations that describe physical laws. It also introduces the Buckingham Pi Theorem, allowing us — under the right conditions — to find the relationship between physical parameters based on such a dimensional analysis. At the end of this post is a link to a jupyter notebook to explore the topic yourself.

The relation between weight and flying speed

One of the examples examines the relationship between physical parameters involved in flight. The book by Henk Tennekes contains a collection of the cruising speed $v$ and weight $W$ of heaps of flying objects and animals. Let us have a look at the data first:

In [ ]:
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

df = pd.read_csv(url,skiprows=5, index_col=False)
speed = df['cruising speed (m/s)'].values
weight = df['weight (newtons)'].values

As the weights $W$ and speeds $v$ of these fliers ranging from small insects to large aircraft span orders of magnitude, it makes sense to plot on log-log scale:

In [8]:
plt.plot(np.log10(speed), np.log10(weight),'k.')

These data appear to follow a line (on this log-log scale), with the largest outliers involving “exotic” fliers such as a large dynosaur and a human-powered aircraft. Nevertheless, a regression on log-log scale of the data, results in:

In [6]:
from scipy.stats import linregress
slope, intercept, r_value, p_value, std_err = linregress(np.log10(speed), np.log10(weight))

# set up an array to plot the regression line from the fitted parameters:
xfid= np.linspace(0,2.5)

plt.plot(np.log10(speed), np.log10(weight),'k.',label='fliers')
plt.plot(xfid, xfid*slope+intercept,label='regression with slope {:1.2f}'.format(slope))

The slope of the best-fitting line is approximately 5.58. Is this a coindidence? No! If you want to learn why, we encourage you to work your way through section 2.3-5 of our book. A dimensional analysis with the Buckingham Pi Theorem shows

$$W \propto v^{6}.$$

On a logarithmic scale, we find that $$\log(W) \propto \log(v^{6}) = 6 \log (v),$$ so that our experimentally determined slope of 5.58 is not half bad! Of course, the relationship of weight and speed of flight involve other variables. if you want to learn how the density of the object and the air, gravity, and the shape and size of the wings come in to play, please read Chapter 2 of our book. And if you want to have a play with these data yourself, here is the jupyter notebook.

large scale real data tomography: teleseismic events only

We perform the inversion using teleseismic events only for comparison. The stations and picks are identical to the data for teleseismic earthquake in this attempt.


Eighteen stations are used, divided into 11 in central Auckland, and 7 surrounding the former array. Low resolution channels are removed.


There are 121 events considered, though 13 are discarded due to bad data. We also have issue with the time window while we were picking.


The inversion shows negative velocity anomaly in most area in the Northland. The negative velocity in the North is a feature of our joint inversion.

There isn’t any change in travel-time residual distribution between initial and final model. The distribution is centred more in the positive direction which explains the mainly negative velocity anomaly.

AVF tomography: local events + northern illumination

Relative travel-time tomographic inversion is performed to the Auckland Volcanic Field (AVF), extending to include the northern parts of the North Island New Zealand. First P-wave arrivals are picked manually from earthquake traces that can be ranged in terms of quality. The travel-time residual between the picked arrival and the background ak135 1D-Earth model is the basis of this tomographic inversion.


Eighteen stations are used in the inversion. The stations can be differentiated into the 11 that are centred to the city of Auckland, and the remaining 7 that surround the former array. For all of the stations, we are not using the low resolution channels.


The events that we use here can be separated into local and teleseismic events. The local events make up the sources that are originated from the South of the stations array. Meanwhile the teleseismic events, mostly from the pacifics, cover the illumination from the North.

There are 116 events sourced locally that we looked at. Most of these events have good trace data, but there are 4 events in which picking is not applicable due to the noise. From the remaining 112, some stations don’t show earthquake occurrence in the traces. This happened to the stations on the outside periphery, and the most times with station OUZ. We think that this may be related to the location of the stations. Station OUZ is the furthest station away from the sources, and thus the time window in which the traces are shown could be out of range here. All data that exhibit these situations are not considered for the inversion. Other traces that are regarded too noisy are also not included.

local events

Similarly with the teleseismic events, there are 121 events taken into account, 13 of which are discarded due to bad data. The same situation with the local data is also observed here. Station HIZ and OUZ at the two ends of the array are the most often disregarded from the inversion because of this. Again on top of all of this, if a trace data is too noisy it would also be discarded.

teleseismic events


To assess the inversion performance, we plot the evolution of standard deviation (in seconds) and chi^2 throughout the iterative steps. We see that after 3 iterations, we won’t get any more significant improvement to the model.The graphs show that the overall improvement from the initial model is less than 10%. Relaxing the damping parameter does reduce the std deviation and chi^2 but only slightly. However doing so might introduce redundant features.

standard deviation and chi^2

inversion result

The following result is obtained after performing 3 iteration steps. We see that the Auckland Volcanic Field is at the edge of a predominantly negative velocity anomaly (in km/s). We are still unsure on what might have caused this observation. But the negative anomaly in Auckland is a characteristic for areas of active volcanic activity. For example in the southern extent of the model, the negative anomalies reflect the locations of Mt. Taranaki, and the Taupo Volcanic Zone.

The red boxes signify the extent of our model that we’re most confident on. However it doesn’t mean that features outside of this box are all artefacts, they are simply less accurate due to the configuration of our sources and stations.

We add here the travel-time residual distribution for the initial and final model after 3 iterations. There is little change as explained by the standard-deviation analysis. The histograms are centred near zero, although we do have travel-times that are very late or very early.

large scale real data tomography: 100+ local events

In this travel-time tomographic inversion we plan to use more events to get a significant result in terms of AVF subsurface velocity field.


Eighteen stations are used in the inversion. These are divided into 11 stations located around the central area, and 7 outer surrounding stations. We removed the low resolution channels associated with the station.


events and picks

We have 116 earthquakes around the region that we plan on using.

Almost all events have good trace data, which means that we have a good quality source input that we can work with. Our next step is to determine the first arrival pick with impeccable accuracy. The automatic trigger does a good job with picking the arrivals to a couple of seconds or so, but to consistently get an accuracy to the second or even less, picking by hand is still the go to method.

While we are doing the picking, some number of traces don’t show the event. We presume that this is mostly due to the location of the related station. The relative distances can make the time window in which the traces are shown becoming too short. For example, the northernmost station OUZ is the one that mostly show this condition. For the following results, the corresponding picks are excluded from the inversion, but we can revisit these picks and include them later on.

Even though the number of traces that don’t show an earthquake is small, it is still small things that we can improve on.


The split in the positive and negative anomaly (in km/s) which we also saw in the inversion with less events can be seen more distinctly this time around. Positive anomaly is prominent to the North and similarly, negative in the South. Another highlight is the negative anomaly right below Auckland (36.6-37.2 S, 174.6-175.2 E) at 50+ km.

We can even see a dip in the anomaly although we cannot be sure how much influence does the smearing effect has here. It is consistent with the primary paths direction from the South-East. However unlike any smearing that we have seen before in previous results, that are usually constricted to the outer and deep part of the model, here it is rather central. For example, the deep lineaments in the NS slice are likely to be smearing, and are in agreement with the three deepest events. Furthermore, our previous analysis from the checker-board test shows that at least for the first 70 km or so, the smearing is not too dominant. Therefore, we may be seeing a jump or a continuation of structure to depth commencing somewhere in the south Auckland.

Here are the plots showing how the root-mean-square travel-time residuals, standard deviation, and chi^2 change through each iteration. The first iteration is when the misfit drops drastically, afterward the plots quickly becoming stable. These plots justify the number of iterations we did as after three or four iterations we already have a convergence.

We have here the histogram of the travel-time residuals for initial and final model. There is no change between the two, but that is consistent with the standard-deviation plot. Most residuals are within a couple of seconds, but we do have some that can be more than 20 s late.

Jami got engaged!

Jami Johnson sent us word of her engagement to her Kiwi  boyfriend Sam! Here is a picture of the happy couple on a recent trip to Norway. We wish Sam and Jami all the best (and for them to move back to New Zealand, of course).

Apple Seismology

Today, our paper on Apple Seismology appeared in Physics Today. It was a follow up on Sam and Zoe’s experiments, highlighting the similarities between seismic waves and normal modes in the Earth, and their acoustic equivalents in a Braeburn apple.

In the panel on the left, seismic (surface) waves traverse Earth, while in the right panel laser-generated and detected surface waves circle the apple.

Joint Inversion – finer resolution checker-board

Our previous joint inversion of assessing the resolution of tomography in the Auckland Volcanic Field gives us promising results for subsurface features at the scale of ~60 km (see this post).

The latest checker-board test to improve from this finding involve using smaller sized checker-board, around half of the previous test. The same receivers configuration and local and teleseismic sources are still used here – 1330 cataloged events (681 local + 649 teleseismic) onto 25 receiver stations.

In comparison to our previous result, this recovery indicates that yes smaller scale features can plausibly be resolved, but at the cost of having the extent of the well recovered model being diminished. Some regions outside the box may also show good checker-board recovery. But these are not always nicely interconnected (ex. to the north and northwest of lake Taupo).

More specifically, this result now implies that:

  • The region where resolution is commendable, is now mostly confined around the Auckland area only. Up to about 100 km deep, EW span of 260 km, and NS extend of 167 km.
  • Within this region, ~30 km (roughly 35 in depth 30 in width and 33 in length) sized subsurface features can be resolved.

Resolution test on the AVF – Joint Inversion

Joint inversion of local and teleseismic events

Promising results are obtained from performing a joint inversion of local and teleseismic sources in a checker-board synthetic test. We proceed with using 1330 cataloged events (681 local + 649 teleseismic) onto 25 receiver stations.

The local events made up the sources coming in from the southeast of the receivers array. Whereas teleseismic sources complement this with incoming rays from other directions but most notably from the Northwest.

The extent of the model encapsulate from Northland reaching to Hawke’s Bay. The implications from this recovered patterns are:

  • Dimension of the region where resolution is most exemplary is signified by the red box (322 x 231 x 270 km). Outside of it, the smearing artefact is too prominent.
  • Within the box, the size of the features that are resolvable, corresponds to the size of each checker-board, is in the order of ~60 km (roughly 60 in depth 62 in width and 67 in length).

Comparing resolution results – Local vs Teleseismic vs Joint Inversion

We compare three sets of results of resolution test: using local sources only, teleseismic sources only, and a joint local and teleseismic. The source locations are mainly situated so that the local sources consist of earthquakes south-southeast of the receivers array, and the teleseismic sources are north-northwest of the array. The raypath plots illustrate the statement.

Initial model is the same in all three tests. Comparing the recovered checker-board for the horizontal slices, we expected to see recovery mainly in the Southeast for local sources and to the Northwest for teleseismic sources. Therefore we see a success in joint inversion, as in the joint model, we can see an improved coverage width extending from Northland to Waikato.

In vertical slices, local sources provide better depth recovery as the results from the teleseismic sources are more smeared. However a joint inversion albeit with added smearing can resolve more feature, for example the shallow positive anomaly at 175.8 E.

The next step is to add the remainder of the sources both local and teleseismic, to be inverted jointly. Once cleared, we need to assess the extent of the resolution i.e. how small of a feature we can resolve in the inversion, which is linked to the size of our checker-board.

Skip to toolbar