Chapter 3 Green’s Functions
3.1 Introduction
Green’s functions (GFs) are commonly used to represent force excitation in wave propagation and electromagnetic, mechanical, and thermal problems. Green’s function poses a significant error factor in seismic waveform inversion (Yagi and Fukahata 2011), and, therefore, the selection of the most suitable model is extremely important.
For Green’s function calculations, we used a 1D velocity model derived from CRUST1.0 (Jerkins et al. 2023). The computation was performed using the Fomosto QSEIS code (Wang (1999), Heimann et al. (2018)).
Figure 3.1: Concept of the Greens function stores.
Figure 3.2: Velocity model used for Green’s function computations (see supplement material in Jerkins et al., 2023).
The number of Green’s function components required to represent the effects of a source at a receiver site is determined by the physics of the modelled process and the symmetry of the medium. There are several Green tensor component schemes and configuration types for GF stores available in Fomosto (see Heimann et al. (2019)). For our computations we used QSEIS backend, elastic10 component scheme, which means that the number of Green’s function components is 10. Since the maximum distances to the stations did not exceed ~3000 km, we can use QSEIS code for synthetic seismograms computations in a cartesian, layered half-space model space. For the larger distances, the QSSP backend is recommended since QSSP includes the spherical Earth. More information regarding the available Fomosto backends can be found here.
| Component scheme | Configuration type | Source type | Receiver distances |
|---|---|---|---|
| elastic10 | A, B | moment tensor | all |
| elastic8 | A, B | moment tensor | far field only |
| elastic5 | A, B | single force | all |
| elastic2 | A, B | isotropic source | all |
| poroelastic10 | A, B | flow rate, pressure | near field (quasi-static) |
| elastic18 | C | moment tensor | all |
3.2 How to make a config file Fomosto
- For Green’s function calculations, the model should have the following information: depth [km], Vp [km/s], Vs [km/s], density [g/cm3], Qp and Qs.
Figure 3.3: Velocity model used for Greens function computations (see supplement material in Jerkins et al., 2023.).
- Create a config file for Fomosto:
Activate gmtobsrocko virtual environment:
conda activate gmtobsrocko
Open make_config_for_fomosto.ipynb file and run the code cell by cell.
First, we import the required packages.
import pandas as pd
import numpy as np
from scipy import stats
import os
After that, we set the paths to the input and output files and read the velocity model in xlsx format using Pandas:
current_dir = os.path.abspath(os.getcwd())
#Insert the name for the velocity model file
input_file_name = "NorthSea1D_CRUSTaverage_model_annie.xlsx"
#Set the output config file name
output_file_name ="config_northsea_GFs"
velocity_model_filepath = os.path.join(current_dir,
input_file_name)
model_df = pd.read_excel(velocity_model_filepath)
model_df.drop(model_df.filter(regex="Unname"),axis=1, inplace=True)
config_output = os.path.join(current_dir,
"config_files_4_fomosto",
output_file_name)
Adjust the parameters of the config file, such as Green’s function store name (config_id), sample rate, receiver_depth, source_depth_min, source_depth_max, distance_min, distance_max, distance_delta if needed.
# Name your Green's function (GF) store
config_id = "northsea"
# Set the parameters for GF store
sample_rate = 5 # Sampling rate in [Hz]
receiver_depth = 0.0 # [m]
source_depth_min = 0.0 # [m]
source_depth_max = 45000.0 # [m]
source_depth_delta = 1000.0 # [m]
distance_min = 0.0 # [m]
distance_max = 2000000.0 # [m]
distance_delta = 5000.0 # [m]
To display currently selected velocity model use:
# display the current model
model_df
Make the config file by running the following cell:
depth = model_df['depth [km]']
p_vel = model_df['Vp[km/s]']
s_vel = model_df['Vs[km/s]']
rho = model_df['density[g/cm^3]']
qp = model_df['qp']
qs = model_df['qs']
cond, cond2, cond3, cond4 = True, True, True, True
with open(config_output, 'w+') as f:
f.write("--- !pf.ConfigTypeA\n")
f.write("id: " + config_id + "\n")
f.write("modelling_code_id: qseis.2006a\n")
f.write("earthmodel_1d: |2\n")
f.write(" " + str(depth[0])[:6])
f.write(" " + str(p_vel[0])[:6])
f.write(" " + str(s_vel[0])[:6])
f.write(" " + str(rho[0])[:6])
f.write(" " + str(qp[0])[:6])
f.write(" " + str(qs[0])[:6])
f.write("\n")
for i in range(1, len(depth)):
if depth[i-1] >=4 and cond:
f.write(" transition S\n")
cond=False
elif depth[i-1] >=5.310 and cond2:
f.write(" transition P\n")
cond2=False
elif depth[i-1] >=6.521 and cond3:
f.write(" baselocal\n")
cond3=False
elif depth[i-1] >=31 and cond4:
f.write(" mantle\n")
cond4=False
f.write(" " + str(depth[i])[:6])
f.write(" " + str(p_vel[i])[:6])
f.write(" " + str(s_vel[i])[:6])
f.write(" " + str(rho[i])[:6])
f.write(" " + str(qp[i])[:6])
f.write(" " + str(qs[i])[:6])
f.write("\n")
f.write("sample_rate: " + str(sample_rate) + "\n")
f.write("component_scheme: elastic10\n")
f.write("tabulated_phases:\n")
f.write("- !pf.TPDef\n")
f.write(" id: begin\n")
f.write(" definition: '10.0'\n")
f.write("- !pf.TPDef\n")
f.write(" id: end\n")
f.write(" definition: '1.0'\n")
f.write("- !pf.TPDef\n")
f.write(" id: anyP\n")
f.write(" definition: P,p\n")
f.write("- !pf.TPDef\n")
f.write(" id: anyS\n")
f.write(" definition: S,s\n")
f.write("- !pf.TPDef\n")
f.write(" id: P\n")
f.write(" definition: P\n")
f.write("- !pf.TPDef\n")
f.write(" id: S\n")
f.write(" definition: S\n")
f.write("- !pf.TPDef\n")
f.write(" id: p\n")
f.write(" definition: p\n")
f.write("- !pf.TPDef\n")
f.write(" id: s\n")
f.write(" definition: s\n")
f.write("ncomponents: 10\n")
f.write("receiver_depth: " + str(receiver_depth) + "\n")
f.write("source_depth_min: " + str(source_depth_min) + "\n")
f.write("source_depth_max: " + str(source_depth_max) + "\n")
f.write("source_depth_delta: " + str(source_depth_delta) + "\n")
f.write("distance_min: " + str(distance_min) + "\n")
f.write("distance_max: " + str(distance_max) + "\n")
f.write("distance_delta: " + str(distance_delta) + "\n")
The created config file should have the following structure:
--- !pf.ConfigTypeA
id: northsea
modelling_code_id: qseis.2006a
earthmodel_1d: |2
0.0 2.06 0.88 1.62 500 250
2.3 2.06 0.88 1.62 500 250
2.3 6.01 3.5 2.72 500 250
11.6 6.01 3.5 2.72 500 250
transition S
11.600 6.51 3.75 2.82 500 250
transition P
21.1 6.51 3.75 2.82 500 250
baselocal
21.1 7.04 3.9 2.98 500 250
30.6 7.04 3.9 2.98 500 250
30.6 8.06 4.48 3.32 1000 500
49.6 8.06 4.48 3.32 1000 500
mantle
49.6 8.25 4.74 3.34 1000 500
94.6 8.25 4.74 3.34 1000 500
94.6 8.5 4.88 3.42 1000 500
294.6 8.5 4.88 3.42 1000 500
sample_rate: 5
component_scheme: elastic10 # number of Green's function components (always use 10 with QSEIS).
tabulated_phases:
- !pf.TPDef
id: begin
definition: '10.0'
- !pf.TPDef
id: end
definition: '1.0'
- !pf.TPDef
id: anyP
definition: P,p
- !pf.TPDef
id: anyS
definition: S,s
- !pf.TPDef
id: P
definition: P
- !pf.TPDef
id: S
definition: S
- !pf.TPDef
id: p
definition: p
- !pf.TPDef
id: s
definition: s
ncomponents: 10
receiver_depth: 0.0
source_depth_min: 0.0
source_depth_max: 45000.0
source_depth_delta: 1000.0
distance_min: 0.0
distance_max: 2000000.0
distance_delta: 5000.0
- To plot the model using Cake application use the following command:
cake plot-model --model=config_northsea_GFs
As a result, this will recreate Figure 3.2.
3.3 How to compute Green’s functions store
- Activate the virtual environment:
conda activate obsrocko
- Initiate Fomosto project:
fomosto init qseis.2006b northsea_gfs
This command will create a folder named northsea_gfs, containing the following folders and files:
northsea_gfs/
|-- config # (1)
`-- extra/
`-- qseis # (2)
Copy config file created in Step 2 to northsea_gfs folder.
To create tabulated phase arrivals run the follwing command:
fomosto ttt
In case of rebuilding the tabulated phase arrivals, a follwoing command should be used instead to overwrite the existing table:
fomosto ttt --force
- To run Green’s function computation use the command:
fomosto build
Parallalisation of the computations between N CPU’s can be done using:
fomosto build --nworkers=N
Without setting the number of CPU’s used for the computations the code will try to use all the available CPU’s.
- Other usefull commands:
For troubleshouting of the computations a command –loglevel=debug can be helpful:
fomosto build --loglevel=debug --nworkers=N
To overwrite the existing Green’s function store use:
fomosto build force
For more information about the structure of the config file and for the clarification on Fomosto package please use the official Fomosto tutorial page.