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.
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 invert3d.in 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.