PDF
Propagation of Seismic Waves Through Earth
Introduction
The analysis of the propagation of seismic waves through the internal structure of Earth is a complex topic. The curvature of Earth, the existence of discontinuities and the variation of the material properties with the depth create a complex interaction that makes the analysis quite challenging. This tutorial presents a method to analyze the propagation of seismic waves through Earth’s internal structure. The model uses a 2D axisymmetric geometry to represent the material discontinuities and the variation of properties through the concentric layers of Earth.
Figure 1: Geometry used in the model.
The model uses the Elastic Waves, Time Explicit and the Pressure Acoustics, Time Explicit interfaces to represent the solid and fluid parts of Earth. The model demonstrates the scalability of time explicit interfaces and their applicability to capture wave propagation in large and very large acoustic models (containing many wavelengths).
Model Definition
A simplified excitation consisting of a pulse is used to study the transmission and propagation of the different pressure and shear waves across the inner structure of Earth. This pulse, sometimes called a tone burst pulse, uses the product of a sine wave with frequency f0 multiplied by a sine wave with a frequency of f0/5. Figure 2 shows this pulse, with the highest frequency content of 0.1 Hz.
Figure 2: Tone burst load.
The variation of properties through the internal structure of Earth is taken from the Ref. 1which is also reproduced in Ref. 2. Figure 3 shows the variation of density, speed of pressure waves, and speed of shear waves with the depth. Figure 3 shows several discontinuities between the outermost layers of Earth. The area between 2900 km and 5100 km of depth is named the outer core of Earth. This outer core is made of molten iron and nickel, and therefore shear waves do not propagate across it.
d
Figure 3: Variation of the material properties with depth.
The model is excited with a vertical force following the pulse previously described. This model disregards the presence of continental and oceanic crust and considers Earth as a perfect sphere with concentric layers. This means that the location (latitude and longitude) of the source is not relevant, as the representation is the same for any point on Earth.
Results and Discussion
Time explicit interfaces impose a strict limit on the time step for stability purposes. This time step is proportional to the overall smallest cell wave time scale (the variable elte.wtc or pate.wtc), shown in Figure 4. The cell wave time scale is defined as the ratio between the element size and the fastest wave propagating through the element.
The mesh needs to resolve the wavelength of the slowest wave traveling through the element as well. This requirement imposes different mesh size on different layers of Earth, which in turn translates into different cell wave time scales. As shown in Figure 4, due to the slow speed of the shear waves of the core, the elements with the smallest cell wave time scales are located in the inner core of Earth.
Figure 4: Cell wave time scale.
Figure 5: 3D representation of Earth showing the velocity.1
The model could be revolved to form a 3D image of Earth, as shown in Figure 5. Then in Figure 6, the velocity at different times of the simulation is depicted. This figure reveals a lot of information about the different waves being transmitted and reflected;
a
Note the two different pressure wave pulses, called p-waves in the context of seismology, the first one coming directly from the source while the second pulse coming from its reflection on Earth surface.
b
Two pulses of shear waves, called s-waves in the context of seismology, start to develop at the source. Both p-waves and s-waves are generated at the same time, but as s-waves travel approximately at half the speed of the p-waves, after a given time it is possible to see the two distinct wave fronts.
c
As the p-waves travel through the mantle and reach the outer core, part of their energy is reflected back as p-waves. This is usually called the reflection from the core-mantle boundary.
d
The p-waves reaching the outer core are partially transmitted as p-waves traveling through the liquid. The p-waves that have not been transmitted to the outer core or reflected back continue their travel. This leaves part of Earth without direct p-waves. This is called the p-waves shadow zone and it was quite relevant to identify the internal structure of Earth at the early stages of the development of seismology.
e
As the s-waves reach the core-mantle boundary, part of their energy is reflected back as p-waves and s-waves. As p-waves travel faster, it is easy to identify them further away from the reflection point.
f
The rest of the reflected energy is transmitted as s-waves.
g
As s-waves cannot travel through fluids, the s-waves reaching the core-mantle boundary are partially reflected and partially transformed into p-waves traveling through the liquid.
h
As the p-waves traveling through the outer core reach the inner core, part of their energy is converted to s-waves transmitted through the inner core. Note the short wavelength of these waves, caused by the slow speed of the s-waves in the inner core.
i
The existence of discontinuities in the outer most layers create a series of refracted waves called head waves or von Schmidt waves.
j
The p-waves transmitted through the inner core continue to travel forming a continuous wavefront with the p-waves that have not reached the inner core.
k
The p-waves coming directly from the source finally reach the surface of Earth. This point defines one end of the p-waves shadow zone.
l
The s-waves continue to travel through the mantle without reaching all of Earth. This area where no s-waves are transmitted is called the s-wave shadow zone. As the shear waves cannot travel through fluids, these are the only direct s-waves still existing.
m
The p-waves reach the core-mantle boundary at the other end of Earth.
n
The s-waves finally reach the surface of Earth. This point defines the start of the s-waves shadow zone, marked in blue.
o
The Rayleigh waves, which are slower than the s-waves, can be seen a this point.
p
The p-waves that have traveled through the core finally reach the other side of Earth, defining the other end of the p-waves shadow zone, marked in red.
Figure 6: Velocity plot at 200 s, 500 s, 600 s, 800 s, 1000 s, and 1200 s.
The velocity experienced above the source and at the probe location is shown Figure 7. This shows how the relatively simple response at the source is transformed in a complex time signal showing the arrival of the different types of waves to the probe.
Figure 7: Probe output.
Notes About the COMSOL Implementation
The model uses material properties that vary with depth. A variable called depth is used to define the depth of any point in the model. Interpolation curves with linear interpolation are used to define the material properties using the expression rho3(depth), for example.
The Material Discontinuity feature is used on all the boundaries between solid layers with sudden changes of material properties.
The mesh uses a different size for each of the layers. This size depends on the speed of the slowest waves transmitted through the layer.
Given the large size of the model, some of the variables are removed from the output. This is done to reduce the size of the model file. Stress variables are also quite useful to discern between p-waves and s-waves, as demonstrated in the Ground Motion After Seismic Event: Scattering off a Small Mountain tutorial.
Probes are a very useful way to obtain output that will not be stored in the solution. It is also a good method to obtain signals with more time resolution than the requested output, as by default they will be populated using all of the times used by the solver.
References
1. B.L.N. Kennett, E.R. Engdahl, and R. Buland. 1995. “Constraints on seismic velocities in the earth from travel times,” Geophys. J. Int., vol. 122, pp. 108–124, doi.org/10.1111/j.1365-246X.1995.tb03540.x.
2. Data Services Products: EMC-ak135-f, doi.org/10.17611/DP/9991801.
Application Library path: Acoustics_Module/Elastic_Waves/seismic_waves_earth
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D Axisymmetric.
2
In the Select Physics tree, select Acoustics > Pressure Acoustics > Pressure Acoustics, Time Explicit (pate).
3
Click Add.
4
In the Select Physics tree, select Acoustics > Elastic Waves > Elastic Waves, Time Explicit (elte).
5
Click Add.
6
Click  Study.
7
In the Select Study tree, select General Studies > Time Dependent.
8
Click  Done.
Import the parameters from an external file.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_parameters.txt.
Create a Gaussian spatial distribution that will be used in the load.
Space Distribution
1
In the Home toolbar, click  Functions and choose Global > Analytic.
2
In the Settings window for Analytic, type Space Distribution in the Label text field.
3
In the Function name text field, type G_space.
4
Locate the Definition section. In the Expression text field, type 10e26/sqrt(pi*dS)*exp(-((r - r0)^2 + (z - z0)^2)/dS).
5
In the Arguments text field, type r,z.
6
Locate the Units section. In the table, enter the following settings:
Argument
Unit
r
m
7
In the Function text field, type 1.
8
Locate the Plot Parameters section. In the table, enter the following settings:
Plot
Argument
Lower limit
Upper limit
Fixed value
Unit
√
r
0
100000[m]
0
m
√
z
z0-50000[m]
z0+50000[m]
0
 
9
Click  Plot.
The image should look like this.
Create a tone burst that will be used in the load.
Time Distribution
1
In the Home toolbar, click  Functions and choose Global > Analytic.
2
In the Settings window for Analytic, type Time Distribution in the Label text field.
3
In the Function name text field, type G_time.
4
Locate the Definition section. In the Expression text field, type if(t<2.5*T0,sin(2*pi*f0*t)*sin(2*pi*f0*t/5),0).
5
In the Arguments text field, type t.
6
Locate the Units section. In the table, enter the following settings:
Argument
Unit
t
s
7
In the Function text field, type N.
8
Locate the Plot Parameters section. In the table, enter the following settings:
Plot
Argument
Lower limit
Upper limit
Fixed value
Unit
√
t
0
10*T0
0
s
9
Click  Plot.
The image should look like Figure 2.
Create and import tables that define the density, speed of pressure waves, and speed of shear waves of the different layers of the Earth.
Rho Layer 3
1
In the Home toolbar, click  Functions and choose Global > Interpolation.
2
In the Settings window for Interpolation, type Rho Layer 3 in the Label text field.
3
Locate the Definition section. In the Function name text field, type rho3.
4
Click  Load from File.
5
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_rho3.txt.
6
Locate the Interpolation and Extrapolation section. From the Extrapolation list, choose Linear.
7
Locate the Units section. In the Argument table, enter the following settings:
Argument
Unit
t
km
8
In the Function table, enter the following settings:
Function
Unit
rho3
kg/m^3
9
Click  Plot.
The image should look like this.
Cp Layer 3
1
Right-click Rho Layer 3 and choose Duplicate.
2
In the Settings window for Interpolation, type Cp Layer 3 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cp3.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cp3.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cp3
m/s
8
Click  Plot.
The image should look like this.
Cs Layer 3
1
Right-click Cp Layer 3 and choose Duplicate.
2
In the Settings window for Interpolation, type Cs Layer 3 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cs3.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cs3.txt.
7
Click  Plot.
The image should look like this.
Rho Layer 4
1
Right-click Cs Layer 3 and choose Duplicate.
2
In the Settings window for Interpolation, type Rho Layer 4 in the Label text field.
3
Locate the Definition section. In the Function name text field, type rho4.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_rho4.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
rho4
kg/m^3
8
Click  Plot.
The image should look like this.
Cp Layer 4
1
Right-click Rho Layer 4 and choose Duplicate.
2
In the Settings window for Interpolation, type Cp Layer 4 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cp4.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cp4.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cp4
m/s
8
Click  Plot.
The image should look like this.
Cs Layer 4
1
Right-click Cp Layer 4 and choose Duplicate.
2
In the Settings window for Interpolation, type Cs Layer 4 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cs4.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cs4.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cs4
m/s
8
Click  Plot.
The image should look like this.
Rho Layer 5
1
Right-click Cs Layer 4 and choose Duplicate.
2
In the Settings window for Interpolation, type Rho Layer 5 in the Label text field.
3
Locate the Definition section. In the Function name text field, type rho5.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_rho5.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
rho5
kg/m^3
8
Click  Plot.
The image should look like this.
Cp Layer 5
1
Right-click Rho Layer 5 and choose Duplicate.
2
In the Settings window for Interpolation, type Cp Layer 5 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cp5.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cp5.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cp5
m/s
8
Click  Plot.
The image should look like this.
Cs Layer 5
1
Right-click Cp Layer 5 and choose Duplicate.
2
In the Settings window for Interpolation, type Cs Layer 5 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cs5.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cs5.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cs5
m/s
8
Click  Plot.
The image should look like this.
Rho Layer 6
1
Right-click Cs Layer 5 and choose Duplicate.
2
In the Settings window for Interpolation, type Rho Layer 6 in the Label text field.
3
Locate the Definition section. In the Function name text field, type rho6.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_rho6.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
rho6
kg/m^3
8
Click  Plot.
The image should look like this.
Cp Layer 6
1
Right-click Rho Layer 6 and choose Duplicate.
2
In the Settings window for Interpolation, type Cp Layer 6 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cp6.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cp6.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cp6
m/s
8
Click  Plot.
The image should look like this.
Rho Layer 7
1
Right-click Cp Layer 6 and choose Duplicate.
2
In the Settings window for Interpolation, type Rho Layer 7 in the Label text field.
3
Locate the Definition section. In the Function name text field, type rho7.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_rho7.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
rho7
kg/m^3
8
Click  Plot.
The image should look like this.
Cp Layer 7
1
Right-click Rho Layer 7 and choose Duplicate.
2
In the Settings window for Interpolation, type Cp Layer 7 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cp7.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cp7.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cp7
m/s
8
Click  Plot.
The image should look like this.
Cs Layer 7
1
Right-click Cp Layer 7 and choose Duplicate.
2
In the Settings window for Interpolation, type Cs Layer 7 in the Label text field.
3
Locate the Definition section. In the Function name text field, type cs7.
4
Click  Clear Table.
5
Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file seismic_waves_earth_cs7.txt.
7
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
cs7
m/s
8
Click  Plot.
The image should look like this.
Group the interpolation curves to facilitate the navigation through the model.
Cp Layer 3 (cp3), Cs Layer 3 (cs3), Rho Layer 3 (rho3)
1
In the Model Builder window, under Global Definitions, Ctrl-click to select Rho Layer 3 (rho3), Cp Layer 3 (cp3), and Cs Layer 3 (cs3).
2
Right-click and choose Group.
Layer 3 Properties
In the Settings window for Group, type Layer 3 Properties in the Label text field.
Cp Layer 4 (cp4), Cs Layer 4 (cs4), Rho Layer 4 (rho4)
1
In the Model Builder window, under Global Definitions, Ctrl-click to select Rho Layer 4 (rho4), Cp Layer 4 (cp4), and Cs Layer 4 (cs4).
2
Right-click and choose Group.
Layer 4 Properties
In the Settings window for Group, type Layer 4 Properties in the Label text field.
Cp Layer 5 (cp5), Cs Layer 5 (cs5), Rho Layer 5 (rho5)
1
In the Model Builder window, under Global Definitions, Ctrl-click to select Rho Layer 5 (rho5), Cp Layer 5 (cp5), and Cs Layer 5 (cs5).
2
Right-click and choose Group.
Layer 5 Properties
In the Settings window for Group, type Layer 5 Properties in the Label text field.
Cp Layer 6 (cp6), Rho Layer 6 (rho6)
1
In the Model Builder window, under Global Definitions, Ctrl-click to select Rho Layer 6 (rho6) and Cp Layer 6 (cp6).
2
Right-click and choose Group.
Layer 6 Properties
In the Settings window for Group, type Layer 6 Properties in the Label text field.
Cp Layer 7 (cp7), Cs Layer 7 (cs7), Rho Layer 7 (rho7)
1
In the Model Builder window, under Global Definitions, Ctrl-click to select Rho Layer 7 (rho7), Cp Layer 7 (cp7), and Cs Layer 7 (cs7).
2
Right-click and choose Group.
Layer 7 Properties
In the Settings window for Group, type Layer 7 Properties in the Label text field.
Now that the material properties have been added to the model, proceed to generate the geometry.
Geometry 1
1
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2
In the Settings window for Geometry, locate the Units section.
3
From the Length unit list, choose km.
Circle 1 (c1)
1
In the Geometry toolbar, click  Circle.
2
In the Settings window for Circle, locate the Size and Shape section.
3
In the Radius text field, type r_earth.
4
In the Sector angle text field, type 180.
5
Locate the Rotation Angle section. In the Rotation text field, type -90.
6
Click to expand the Layers section. In the table, enter the following settings:
Layer name
Thickness (km)
Layer 1
th1
Layer 2
th2
Layer 3
th3
Layer 4
th4
Layer 5
th5
Layer 6
th6
Layer 7
th7
7
Click  Build Selected.
8
Click the  Zoom Extents button in the Graphics toolbar.
Add a point at the source location to make sure that the mesh captures this point.
Point 1 (pt1)
1
In the Geometry toolbar, click  Point.
2
In the Settings window for Point, locate the Point section.
3
In the r text field, type r0.
4
In the z text field, type z0.
5
Click  Build Selected.
Add a point away from the source that will be used as a probe to output the resulting pressure and velocity.
Point 2 (pt2)
1
In the Geometry toolbar, click  Point.
2
In the Settings window for Point, locate the Point section.
3
In the r text field, type r_earth*sin(phi0).
4
In the z text field, type r_earth*cos(phi0).
5
Click  Build All Objects.
The geometry should look like Figure 1.
Definitions
Variables 1
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
In the table, enter the following settings:
Name
Expression
Unit
Description
depth
r_earth-sqrt(r^2+z^2)
m
Depth
This variable defines the depth. The material models require this depth to obtain the material properties at a given point.
Given the size of the model, some of the variables will not be stored in the output. Create a few probes where all variables can be obtained, even those not stored in the output.
VV - Source
1
In the Definitions toolbar, click  Probes and choose Point Probe.
2
In the Settings window for Point Probe, type VV - Source in the Label text field.
3
Locate the Source Selection section. Click  Clear Selection.
4
Select Point 18 only.
5
Locate the Expression section. In the Expression text field, type v2z.
Note that the probes at the source will experience a much severe movement, so it makes sense to adapt the units of the probes to account for this.
6
From the Table and plot unit list, choose mm/s.
7
Select the Description checkbox. In the associated text field, type Vertical Velocity.
VV- Probe
1
Right-click VV - Source and choose Duplicate.
2
In the Settings window for Point Probe, type VV- Probe in the Label text field.
3
Locate the Source Selection section. Click  Clear Selection.
4
Select Point 21 only.
5
Locate the Expression section. In the Expression text field, type v2r*sin(phi0)+v2z*cos(phi0).
Change the unit of this probe as the expected pressure is a thousand times lower than the one in the source.
6
From the Table and plot unit list, choose µm/s.
HV - Probe
1
Right-click VV- Probe and choose Duplicate.
2
In the Settings window for Point Probe, type HV - Probe in the Label text field.
3
Locate the Expression section. In the Expression text field, type v2r*cos(phi0)-v2z*sin(phi0).
4
Select the Description checkbox. In the associated text field, type Horizontal Velocity.
The outer core of the Earth is liquid, so it is captured in the model through the Pressure Acoustics, Time Explicit physics.
Pressure Acoustics, Time Explicit (pate)
1
In the Model Builder window, under Component 1 (comp1) click Pressure Acoustics, Time Explicit (pate).
2
In the Settings window for Pressure Acoustics, Time Explicit, locate the Domain Selection section.
3
Click  Clear Selection.
4
Select Domains 6 and 10 only.
The rest of the layers of the Earth are solid, so they are captured in the model through the Elastic Waves, Time Explicit physics.
Elastic Waves, Time Explicit (elte)
1
In the Model Builder window, under Component 1 (comp1) click Elastic Waves, Time Explicit (elte).
2
Select Domains 1–5, 7–9, and 11–15 only.
Elastic Waves, Time Explicit Model 1
1
In the Model Builder window, under Component 1 (comp1) > Elastic Waves, Time Explicit (elte) click Elastic Waves, Time Explicit Model 1.
2
In the Settings window for Elastic Waves, Time Explicit Model, locate the Linear Elastic Material section.
3
From the Specify list, choose Pressure-wave and shear-wave speeds.
Create a body load that will be used as an earthquake in the model.
Body Load 1
1
In the Physics toolbar, click  Domains and choose Body Load.
2
Select Domain 13 only.
3
In the Settings window for Body Load, locate the Force section.
4
From the Load type list, choose Total force.
5
Specify the Ftot vector as
0
r
G_space(r,z)*G_time(t)
z
Use the Material Discontinuity feature on each of the interior boundaries where there is a sudden change in material properties.
Material Discontinuity 1
1
In the Physics toolbar, click  Boundaries and choose Material Discontinuity.
2
Select Boundaries 26–29 and 36–39 only.
Proceed to create the different materials existing in the model. The outermost and innermost layers have constant properties, while the rest of the materials will have properties that depend on depth.
Materials
Layer 1
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, type Layer 1 in the Label text field.
3
Locate the Geometric Entity Selection section. Click  Clear Selection.
4
Select Domains 1 and 15 only.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Pressure-wave speed
cp
cp1
m/s
Pressure-wave and shear-wave speeds
Shear-wave speed
cs
cs1
m/s
Pressure-wave and shear-wave speeds
Density
rho
rho1
kg/m³
Basic
Layer 2
1
Right-click Layer 1 and choose Duplicate.
2
In the Settings window for Material, type Layer 2 in the Label text field.
3
Locate the Geometric Entity Selection section. Click  Clear Selection.
4
Select Domains 2 and 14 only.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
rho2
kg/m³
Basic
Pressure-wave speed
cp
cp2
m/s
Pressure-wave and shear-wave speeds
Shear-wave speed
cs
cs2
m/s
Pressure-wave and shear-wave speeds
Layer 3
1
Right-click Layer 2 and choose Duplicate.
2
In the Settings window for Material, type Layer 3 in the Label text field.
3
Locate the Geometric Entity Selection section. Click  Clear Selection.
4
Select Domains 3 and 13 only.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
rho3(depth)
kg/m³
Basic
Pressure-wave speed
cp
cp3(depth)
m/s
Pressure-wave and shear-wave speeds
Shear-wave speed
cs
cs3(depth)
m/s
Pressure-wave and shear-wave speeds
Layer 4
1
Right-click Layer 3 and choose Duplicate.
2
In the Settings window for Material, type Layer 4 in the Label text field.
3
Locate the Geometric Entity Selection section. Click  Clear Selection.
4
Select Domains 4 and 12 only.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
rho4(depth)
kg/m³
Basic
Pressure-wave speed
cp
cp4(depth)
m/s
Pressure-wave and shear-wave speeds
Shear-wave speed
cs
cs4(depth)
m/s
Pressure-wave and shear-wave speeds
Layer 5
1
Right-click Layer 4 and choose Duplicate.
2
In the Settings window for Material, type Layer 5 in the Label text field.
3
Locate the Geometric Entity Selection section. Click  Clear Selection.
4
Select Domains 5 and 11 only.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
rho5(depth)
kg/m³
Basic
Pressure-wave speed
cp
cp5(depth)
m/s
Pressure-wave and shear-wave speeds
Shear-wave speed
cs
cs5(depth)
m/s
Pressure-wave and shear-wave speeds
Layer 6
1
In the Model Builder window, right-click Materials and choose Blank Material.
2
In the Settings window for Material, type Layer 6 in the Label text field.
3
Select Domains 6 and 10 only.
4
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
rho6(depth)
kg/m³
Basic
Speed of sound
c
cp6(depth)
m/s
Basic
Layer 7
1
In the Model Builder window, under Component 1 (comp1) > Materials right-click Layer 5 (mat5) and choose Duplicate.
2
In the Settings window for Material, type Layer 7 in the Label text field.
3
Locate the Geometric Entity Selection section. Click  Clear Selection.
4
Select Domains 7 and 9 only.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
rho7(depth)
kg/m³
Basic
Pressure-wave speed
cp
cp7(depth)
m/s
Pressure-wave and shear-wave speeds
Shear-wave speed
cs
cs7(depth)
m/s
Pressure-wave and shear-wave speeds
Layer 8
1
In the Model Builder window, right-click Materials and choose Blank Material.
2
In the Settings window for Material, type Layer 8 in the Label text field.
3
Locate the Geometric Entity Selection section. Click  Clear Selection.
4
Select Domain 8 only.
5
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Pressure-wave speed
cp
cp8
m/s
Pressure-wave and shear-wave speeds
Shear-wave speed
cs
cs8
m/s
Pressure-wave and shear-wave speeds
Density
rho
rho8
kg/m³
Basic
Add the Acoustic-Structure Boundary, Time Explicit multiphysics coupling to connect both physics.
Multiphysics
Acoustic–Structure Boundary, Time Explicit 1 (asbte1)
1
In the Physics toolbar, click  Multiphysics Couplings and choose Boundary > Acoustic–Structure Boundary, Time Explicit.
2
Select Boundaries 30, 31, 34, and 35 only.
Proceed to generate the mesh. To limit the size of the model, each of the layers will use a different size that is driven by the wavelength of the slowest wave traveling through that layer divided by 1.5.
Mesh 1
Mapped 1
1
In the Mesh toolbar, click  Mapped.
2
In the Settings window for Mapped, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 1, 2, 14, and 15 only.
Size 1
1
Right-click Mapped 1 and choose Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
Use a slightly finer mesh in the outermost layer, as the Rayleigh waves are slightly slower than the shear waves.
4
Locate the Element Size Parameters section.
5
Select the Maximum element size checkbox. In the associated text field, type cs1/f0/2.0.
6
Select the Minimum element size checkbox. In the associated text field, type cs1/f0/2.0.
7
Click  Build Selected.
Size
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1 click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Coarser.
Convert 1
1
In the Mesh toolbar, click  Modify and choose Convert.
2
In the Settings window for Convert, click  Build Selected.
Free Triangular 1
In the Mesh toolbar, click  Free Triangular.
Size 1
1
Right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 3 and 13 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type cs3(th1+th2)/f0/1.5.
Size 2
1
In the Model Builder window, right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 4 and 12 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type cs4(th1+th2+th3)/f0/1.5.
Size 3
1
Right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 5 and 11 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type cs5(th1+th2+th3+th4)/f0/1.5.
Size 4
1
Right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 6 and 10 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type cp6(th1+th2+th3+th4+th5)/f0/1.5.
Size 5
1
Right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 7 and 9 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type cs7(th1+th2+th3+th4+th5+th6)/f0/1.5.
Size 6
1
Right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 8 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type cs8/f0/1.5.
8
Click  Build All.
The mesh should look like this.
Study 1
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,100,1400).
4
In the Model Builder window, click Study 1.
5
In the Settings window for Study, locate the Study Settings section.
6
Clear the Generate default plots checkbox.
7
Clear the Generate convergence plots checkbox.
8
In the Study toolbar, click  Get Initial Value.
Results
1
In the Model Builder window, expand the Results > Datasets node, then click Results.
2
In the Settings window for Results, locate the Save Data in the Model section.
3
From the Save plot data list, choose On.
Study 1/Solution 1 (3) (sol1)
1
In the Model Builder window, under Results > Datasets right-click Study 1/Solution 1 (sol1) and choose Duplicate.
This dataset is created only for postprocessing purposes.
Selection
1
In the Results toolbar, click  Attributes and choose Selection.
2
In the Settings window for Selection, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundaries 25–31 and 34–41 only.
Mirror 2D 1
In the Results toolbar, click  More Datasets and choose Mirror 2D.
Mirror 2D 2
1
In the Results toolbar, click  More Datasets and choose Mirror 2D.
2
In the Settings window for Mirror 2D, locate the Data section.
3
From the Dataset list, choose Study 1/Solution 1 (3) (sol1).
Revolution 2D 1
1
In the Results toolbar, click  More Datasets and choose Revolution 2D.
2
In the Settings window for Revolution 2D, click to expand the Revolution Layers section.
3
In the Start angle text field, type -90.
4
In the Revolution angle text field, type 225.
Probe Plot
1
In the Model Builder window, under Results click Probe Plot Group 1.
2
In the Settings window for 1D Plot Group, type Probe Plot in the Label text field.
3
Locate the Legend section. From the Position list, choose Manual.
4
In the x-position text field, type 0.15.
5
In the y-position text field, type 0.1.
Result Templates
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Study 1/Solution 1 (1) (sol1) > Acoustic–Structure Boundary, Time Explicit 1 > Cell Wave Time Scale (asbte1).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
This plot, from the Result Templates, for the multiphysics coupling shows the cell wave time scale for both physics in a single plot.
Results
Cell Wave Time Scale (asbte1)
1
In the Settings window for 2D Plot Group, locate the Data section.
2
From the Dataset list, choose Mirror 2D 1.
3
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
4
Locate the Color Legend section. Select the Show maximum and minimum values checkbox.
5
Select the Show units checkbox.
6
Click to expand the Number Format section. Select the Manual color legend settings checkbox.
7
In the Precision text field, type 6.
Surface 1
1
In the Model Builder window, expand the Cell Wave Time Scale (asbte1) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Rainbow.
4
In the Cell Wave Time Scale (asbte1) toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
The image should look like Figure 4.
Result Templates
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Study 1/Solution 1 (1) (sol1) > Elastic Waves, Time Explicit > Velocity Magnitude (elte).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Velocity Magnitude (elte)
1
In the Settings window for 2D Plot Group, locate the Data section.
2
From the Dataset list, choose Mirror 2D 2.
Surface 1
1
In the Model Builder window, expand the Velocity Magnitude (elte) node, then click Surface 1.
2
In the Settings window for Surface, locate the Data section.
3
From the Dataset list, choose Mirror 2D 1.
4
Locate the Expression section. In the Expression text field, type if(isnan(pate.v_inst),elte.vel,pate.v_inst).
This expressions shows the velocity for both physics in a single plot.
5
From the Unit list, choose µm/s.
6
Select the Description checkbox. In the associated text field, type Velocity Magnitude.
7
Click to expand the Range section. Select the Manual color range checkbox.
8
In the Minimum text field, type 0.
9
In the Maximum text field, type 50.
10
Locate the Coloring and Style section. From the Color table list, choose GrayScale.
11
From the Color table transformation list, choose Reverse.
12
In the Velocity Magnitude (elte) toolbar, click  Plot.
13
Click the  Zoom Extents button in the Graphics toolbar.
Velocity (3D)
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Velocity (3D) in the Label text field.
3
Locate the Color Legend section. Select the Show units checkbox.
Surface 1
1
Right-click Velocity (3D) and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type 0.
4
Click to expand the Title section. From the Title type list, choose None.
5
Click to expand the Quality section. From the Evaluation settings list, choose Manual.
6
From the Resolution list, choose Custom.
7
In the Element refinement text field, type 6.
Image 1
1
Right-click Surface 1 and choose Image.
2
In the Settings window for Image, locate the File section.
3
In the Filename text field, type data:///physics/images/earth.jpg.
4
Locate the Mapping section. From the Mapping list, choose Spherical.
5
Find the Axis subsection. From the Axis type list, choose Cartesian.
6
In the x text field, type 0.18.
7
In the y text field, type 0.45.
8
In the z text field, type 1.
9
Find the Angle subsection. In the Rotation text field, type 270.
10
In the Velocity (3D) toolbar, click  Plot.
11
Click the  Show Grid button in the Graphics toolbar.
Selection 1
1
In the Model Builder window, right-click Surface 1 and choose Selection.
2
Select Domains 1 and 15 only.
3
In the Settings window for Selection, locate the Revolution Selection section.
4
Clear the Evaluate the start cap checkbox.
5
Clear the Evaluate the end cap checkbox.
Surface 2
1
In the Model Builder window, right-click Velocity (3D) and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type elte.vel.
4
Select the Description checkbox. In the associated text field, type Velocity.
5
From the Unit list, choose µm/s.
6
Click to expand the Range section. Select the Manual color range checkbox.
7
In the Minimum text field, type 0.
8
In the Maximum text field, type 50.
9
Locate the Coloring and Style section. From the Color table list, choose Thermal.
10
Locate the Quality section. From the Evaluation settings list, choose Manual.
11
From the Resolution list, choose Custom.
12
In the Element refinement text field, type 6.
Selection 1
Right-click Surface 2 and choose Selection.
Material Appearance 1
1
In the Model Builder window, right-click Surface 2 and choose Material Appearance.
2
In the Settings window for Material Appearance, locate the Appearance section.
3
From the Appearance list, choose Custom.
4
From the Material type list, choose Soil.
5
Locate the Color section. Select the Use the plot’s color checkbox.
Selection 1
1
In the Model Builder window, click Selection 1.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose All domains.
4
Locate the Revolution Selection section. Clear the Evaluate the mantle checkbox.
Surface 3
1
In the Model Builder window, right-click Velocity (3D) and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type pate.v_inst.
4
From the Unit list, choose µm/s.
5
Locate the Title section. From the Title type list, choose None.
6
Locate the Range section. Select the Manual color range checkbox.
7
In the Maximum text field, type 50.
8
Locate the Coloring and Style section. From the Coloring list, choose Gradient.
9
From the Top color list, choose Black.
10
From the Bottom color list, choose Custom.
11
Click Define custom colors.
12
Set the RGB values to 253, 185, and 19, respectively.
13
Click Add to custom colors.
14
Click Show color palette only or OK on the cross-platform desktop.
15
Clear the Color legend checkbox.
Selection 1
Right-click Surface 3 and choose Selection.
Material Appearance 1
1
In the Model Builder window, right-click Surface 3 and choose Material Appearance.
2
In the Settings window for Material Appearance, locate the Appearance section.
3
From the Appearance list, choose Custom.
4
From the Material type list, choose Soil.
5
Locate the Color section. Select the Use the plot’s color checkbox.
Selection 1
1
In the Model Builder window, click Selection 1.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose All domains.
4
Locate the Revolution Selection section. Clear the Evaluate the mantle checkbox.
Increase the table size as there will be more than 10000 cells in the probe table.
Probe Table 1
1
In the Model Builder window, expand the Results > Tables node, then click Probe Table 1.
2
In the Settings window for Table, locate the Display section.
3
In the Maximum number of rows displayed text field, type 20000.
Turn off storage of all the variables except for the velocity in both physics.
Study 1
Solver Configurations
In the Model Builder window, expand the Study 1 > Solver Configurations node.
Solution 1 (sol1)
1
In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) node, then click Dependent Variables 1.
2
In the Settings window for Dependent Variables, locate the General section.
3
From the Defined by study step list, choose User defined.
4
In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 node, then click Eigenvalues, Structural (comp1.asbte1.eig).
5
In the Settings window for Field, locate the General section.
6
From the Store in output list, choose None.
7
In the Model Builder window, click Eigenvectors, Structural (comp1.asbte1.veig).
8
In the Settings window for Field, locate the General section.
9
From the Store in output list, choose None.
10
In the Model Builder window, click Strain Tensor, Voigt Notation (comp1.e).
11
In the Settings window for Field, locate the General section.
12
From the Store in output list, choose None.
13
In the Model Builder window, click Eigenvalues, Downside (comp1.elte.mde1.eigd).
14
In the Settings window for Field, locate the General section.
15
From the Store in output list, choose None.
16
In the Model Builder window, click Eigenvalues, Upside (comp1.elte.mde1.eigu).
17
In the Settings window for Field, locate the General section.
18
From the Store in output list, choose None.
19
In the Model Builder window, click Eigenvectors, Downside (comp1.elte.mde1.veigd).
20
In the Settings window for Field, locate the General section.
21
From the Store in output list, choose None.
22
In the Model Builder window, click Eigenvectors, Upside (comp1.elte.mde1.veigu).
23
In the Settings window for Field, locate the General section.
24
From the Store in output list, choose None.
25
In the Model Builder window, click Acoustic Pressure (comp1.p).
26
In the Settings window for Field, locate the General section.
27
From the Store in output list, choose None.
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, click to expand the Results While Solving section.
3
Select the Plot checkbox.
Plot the velocity while the analysis is computing, as this will not slow down the calculation and will allow you to check how the simulation progresses.
4
In the table, enter the following settings:
Plot group
Plot window
Velocity (3D)
Graphics
The analysis takes around 10 hours in a workstation. Time Explicit physics are quite well suited to parallel run, so it is highly recommended to run this model in a cluster using several nodes.
5
In the Study toolbar, click  Compute.
Results
Material Properties
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Material Properties in the Label text field.
3
Locate the Data section. From the Time selection list, choose Last.
4
Click to expand the Title section. From the Title type list, choose Label.
5
Locate the Plot Settings section. Select the Two y-axes checkbox.
6
Select the x-axis label checkbox. In the associated text field, type Depth (km).
7
Select the y-axis label checkbox. In the associated text field, type Speed of waves (km/s).
8
Select the Secondary y-axis label checkbox. In the associated text field, type Density (g/cm<sup>3</sup>).
9
Locate the Legend section. From the Position list, choose Middle right.
Line Graph 1
1
Right-click Material Properties and choose Line Graph.
2
Select Boundaries 1–8 only.
3
In the Settings window for Line Graph, locate the y-Axis Data section.
4
In the Expression text field, type if(isnan(pate.c),elte.cp,pate.c).
5
From the Unit list, choose km/s.
6
Click to expand the Title section. From the Title type list, choose None.
7
Locate the x-Axis Data section. From the Parameter list, choose Expression.
8
In the Expression text field, type depth.
9
Click to expand the Coloring and Style section. From the Width list, choose 2.
10
Click to expand the Legends section. Select the Show legends checkbox.
11
From the Legends list, choose Manual.
12
In the table, enter the following settings:
Legends
Pressure wave
13
Click to expand the Quality section. From the Evaluation settings list, choose Manual.
14
From the Resolution list, choose Extra fine.
15
From the Smoothing list, choose Everywhere.
16
In the Material Properties toolbar, click  Plot.
Line Graph 2
1
Right-click Line Graph 1 and choose Duplicate.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type if(isnan(pate.c),elte.cs,0).
4
Locate the Legends section. In the table, enter the following settings:
Legends
Shear wave
Line Graph 3
1
Right-click Line Graph 2 and choose Duplicate.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type if(isnan(pate.rho),elte.rho,pate.rho).
4
From the Unit list, choose g/cm^3.
5
Locate the y-Axis section. Select the Plot on secondary y-axis checkbox.
6
Locate the Legends section. In the table, enter the following settings:
Legends
Density
Annotation 1
1
In the Model Builder window, right-click Material Properties and choose Annotation.
2
In the Settings window for Annotation, locate the y-Axis section.
3
Select the Plot on secondary y-axis checkbox.
4
Locate the Coloring and Style section. Clear the Show point checkbox.
5
In the Material Properties toolbar, click  Plot.
Annotation 2
1
Right-click Annotation 1 and choose Duplicate.
2
In the Settings window for Annotation, locate the Position section.
3
In the z text field, type 13.65.
4
In the Material Properties toolbar, click  Plot.
The image should look like Figure 3.
Cycle through the different plot groups to reproduce the figures in the results section.
Probe Plot
1
In the Model Builder window, under Results click Probe Plot.
2
In the Probe Plot toolbar, click  Plot.
The results should look like Figure 7.
Surface 1
1
In the Model Builder window, under Results > Velocity Magnitude (elte) click Surface 1.
2
In the Settings window for Surface, locate the Data section.
3
From the Solution parameters list, choose From parent.
4
In the Velocity Magnitude (elte) toolbar, click  Plot.
Cycle through the different times to reproduce the Figure 6.
Velocity Magnitude (elte)
1
In the Model Builder window, click Velocity Magnitude (elte).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Time (s) list, choose 1400.
4
In the Velocity Magnitude (elte) toolbar, click  Plot.
Velocity (3D)
1
In the Model Builder window, click Velocity (3D).
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Time (s) list, choose 1000.
4
In the Velocity (3D) toolbar, click  Plot.
The results should look like Figure 5.
1
NASA Goddard Space Flight Center Image by Reto Stöckli (land surface, shallow water, clouds). Enhancements by Robert Simmon (ocean color, compositing, 3D globes, animation). Data and technical support: MODIS Land Group; MODIS Science Data Support Team; MODIS Atmosphere Group; MODIS Ocean Group Additional data: USGS EROS Data Center (topography); USGS Terrestrial Remote Sensing Flagstaff Field Center (Antarctica); Defense Meteorological Satellite Program (city lights).