The code is for adjoint-state traveltime tomography (ATT) in the Cartesian coordinates. This simple version is mainly for training purpose.
The details of the algorithm can be found at Tong, P. (2021). Adjoint-state traveltime tomography: Eikonal equation-based methods and application to the Anza area in southern California. Journal of Geophysical Research: Solid Earth, 126, e2021JB021818. https://doi.org/10.1029/2021JB021818.
In order to run the codes, you have to install GNU Fortran [Operating system setup tutorials in Chinese]:
# Fedora
$ sudo dnf install gcc-gfortran
# CentOS
$ sudo yum install gfortran
# Ubuntu/Debian
$ sudo apt install gfortran
# macOS
$ brew install gfortran
Generic Mapping Tools (GMT) is used to plot figures. Refer to [GMT reference book in Chinese] for GMT installation and usage.
$ cd data
Prepare your data in dataInFormat_D_xyz
following the designed format.
I use syntheticData.f90
to generate synthetic earthquakes and seismic stations.
$ cd dataSelection/0_commandCenter/
Set the parameters in parametersData.F90
.
You can select the regions for earthquakes and seismic stations.
The two regions can be different.
$ ./workflowStep.sh
$ cd ../../../commandCenter
You should have the following files ready at ../data/
:
sources
receivers
nevtstaray
evtstaminmax
traveltimeReceiverGathers
Set the parameters in parametersGenerator.F90
. The dimension of the forward grid,
number of multiple grids, iteration number and some others can be defined here.
You can compile and execute parametersGenerator.F90
at this time.
$ gfortran -o xpara parametersGenerator.F90
$ ./xpara
$ rm xpara
cd ../model
If this is for a recovery test, you can set up the target velocity model by editing velocity3d_true.F90
.
$ gfortran -o xvel velocity3d_true.F90
$ ./xvel
$ rm xvel
Otherwise, just define the initial model for seismic tomography by editing velocity3d.F90
.
$ gfortran -o xvel velocity3d.F90
$ ./xvel
$ rm xvel
cd ../mesh
Edit memeshgenerator.F90
to adjust the size of the inversion grid.
The following three parameters should be properly set.
invx = 12
invy = 3
invz = 11
In general, we sample one-wavelength anomaly by 5 grid points or the grid interval is about one fourth of the anomaly wavelength.
$ gfortran -o xmesh meshgenerator.F90
$ ./xmesh
$ rm xmesh regmesh.mod
cd ../commandCenter
# Run it if this is a recovery test
$ ./workflow_obstime.sh
$ ./workflow_inversion.sh
The obtained velocity model is located at ../inversion/
as velocity3d015
.
You can display the results along the cross-section set by lineEnds
:
$ ./plot-cross-section.sh
The Bash script will call zCutVelocity.f90
and the gmt script vcut.gmt
.