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)).

Concept of the Greens function stores.

Figure 3.1: Concept of the Greens function stores.

Velocity model used for Green’s function computations (see supplement material in Jerkins et al., 2023).

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

  1. 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.
Velocity model used for Greens function computations (see supplement material in Jerkins et al., 2023.).

Figure 3.3: Velocity model used for Greens function computations (see supplement material in Jerkins et al., 2023.).

  1. 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
  1. 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

  1. Activate the virtual environment:
conda activate obsrocko
  1. 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)
  1. Copy config file created in Step 2 to northsea_gfs folder.

  2. 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
  1. 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.

  1. 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.

References

Heimann, Sebastian, Marius Isken, Daniela Kühn, Henriette Sudhaus, Andreas Steinberg, Simon Daout, Simone Cesca, Hannes Vasyura-Bathke, and Torsten Dahm. 2018. Grond - A Probabilistic Earthquake Source Inversion Framework. GFZ Data Services. https://doi.org/10.5880/GFZ.2.1.2018.003.
Heimann, Sebastian, Hannes Vasyura-Bathke, Henriette Sudhaus, Marius Paul Isken, Marius Kriegerowski, Andreas Steinberg, and Torsten Dahm. 2019. “A Python Framework for Efficient Use of Pre-Computed Green’s Functions in Seismological and Other Physical Forward and Inverse Source Problems.” Solid Earth 10 (6): 1921–35. https://doi.org/10.5194/se-10-1921-2019.
Jerkins, Annie E., Volker Oye, Celso Alvizuri, Felix Halpaap, and Tormod Kværna. 2023. “The 21 March 2022 Mw 5.1 Tampen Spur Earthquake, North Sea: Location, Moment Tensor, and Context.” Bulletin of the Seismological Society of America 114 (2): 741–57. https://doi.org/10.1785/0120230163.
Wang, Rongjiang. 1999. “A Simple Orthonormalization Method for Stable and Efficient Computation of Green’s Functions.” Bulletin of the Seismological Society of America 89 (3): 733–41. https://doi.org/10.1785/bssa0890030733.
Yagi, Yuji, and Yukitoshi Fukahata. 2011. “Introduction of Uncertainty of Green’s Function into Waveform Inversion for Seismic Source Processes.” Geophysical Journal International 186 (2): 711–20. https://doi.org/10.1111/j.1365-246x.2011.05043.x.