Pressure, Plate Motion, and Earthquake Size

IN THE SPACE BELOW, WRITE OUT IN FULL AND THEN SIGN THE HONOR PLEDGE:

“I pledge my honor that I have not violated the honor code during this examination.”

PRINT NAME:

If a fellow student has contributed significantly to this work, please acknowledge them here:

Peer(s):

Contribution:

By uploading this assignment through Canvas, I sign off on the document below electronically.

---

Part I: Pressure variations on and within a terrestrial planet

Pressure, force per unit area, plays an important role in making a habitable planet. Atmospheric pressure controls the living condition on the planet’s surface while the pressure deep inside the planet dictate which materials at which states exist in the planet. We study pressure on Earth and other planets to understand the dynamic of the planets and find the condition that life may prosper.

Section A: pressure at the surface

Pressure at the surface mostly comes from the atmosphere layer on top of the rocky planet. The surface pressure dictates whether liquid water can exist on the planet. The variation of the pressure in space and time causes the circulation of air and more importantly, weather.

The air pressure at the Earth’s surface is \(1.013 \times 10^5 \text{ N/m}^2\). This pressure supports the weight of the Earth’s atmosphere. (Weight is the gravitational force: \(\text{weight} = m \times g\).)

TO DO

Question 1 What is the mass of Earth’s atmosphere? The radius of the Earth: \(6371 \text{ km}\). The gravitational acceleration \(g: 9.8 \text{ m/s}^2\). You assume that weight of the Earth’s atmosphere acts equally on Earth’s surface.

Answer

Question 2 The planets closest to us, Mars and Venus, are in many ways similar to Earth, but their atmospheres are quite different.

Planet	Venus	Mars
Radius	6052 km	3397 km
Gravitational acceleration	9.1 m/s\(^2\)	3.8 m/s\(^2\)
Mass of atmosphere	4.5 x 10\(^{20}\) kg	2.2 x 10\(^{16}\) kg

What is the atmospheric pressure on the surface of Mars and Venus?

Answer

Section B: pressure inside the rocky interior

While measuring pressure on the surface is easy, measuring deep inside the interior is challenging if possible, because we may not directly access the environment down below. Instead, we determine the pressure from other physical properties including planet’s density which can be determined from seismological observations and astronomical observations.

Dziewonski and Anderson (1980) made the Preliminary Reference Earth Model (PREM) from normal mode observations (mode of the Earth’s vibrations), travel time observations for seismic waves and astronomic and geodetic data. PREM contains physical properties such as density, seismic wave velocities, and other elastic properties at specific depths including major boundaries such as Mohorovičić discontinuity and core-mantle bounary.

Image: Seismic velocities (\(\alpha\) and \(\beta\)) and density (\(\rho\)) for the Preliminary Reference Earth Model (PREM).

In this problem set, we will use PREM to calculate the mass of the Earth, gravitational acceleration, and pressure at any given depth. We have demonstrated how to calculate the pressure inside a planet \(P(r)\) from the planet’s mass density \(\rho = \rho(r)\) and the pressure at the planet’s surface \(P(R)\).

(1)\[\begin{equation} P(r) = P(R) + \int_r^R \rho(a) g(a) da \end{equation}\]

The gravitational acceration \(g(r)\) inside a planet can be determined using Newton’s law of gravitation with mass inside the sphere of radius \(r\):

(2)\[\begin{equation} g(r) = \frac{GM_{\text{inside}}}{r^2} = \frac{G}{r^2} \int_0^r 4\pi a^2 \rho(a) da \end{equation}\]

Note that the density must satisfy the equation

(3)\[\begin{equation} M_E = \int_0^R 4\pi r^2 \rho(r) dr \end{equation}\]

where \(M_E\) is the planet’s mass and \(R\) is the planet radius. If you wonder about this equation, it is essentially \(\text{mass} = \text{density} \times \text{volume}\) integrated over the entire sphere.

import pint
import pandas as pd
import pint_pandas

def read_mineos_cards(file,header = 3, R = None):
    """
    Read a card deck file of physical properties in mineos format
    
    Input Parameters:
    ----------------
    
    file: mineos card file containing columns with various properties
    
    header: number of lines in the header
    
    R: Mean radius of the planet
    """

    # Get the default unit registry e.g. MKS units
    ureg = pint.get_application_registry()
        
    # set default radius as Earth
    if R is None: R = 6371000.0 * ureg.meter

    names=['radius','rho','vpv','vsv','qkappa','qmu','vph','vsh','eta']
    units =['m','kg/m^3','m/s','m/s','dimensionless','dimensionless','m/s','m/s','dimensionless']
    fields=list(zip(names,units))
    #formats=[np.float for ii in range(len(fields))]
    # modelarr = np.genfromtxt(file,dtype=None,comments='#',skip_header=3,names=fields)
    modelarr = pd.read_csv(file,skiprows=header,comment='#',sep='\s+',names=fields)

    # read the units from last header
    modelarr_ = modelarr.pint.quantify(level=-1)
    
    # Get the depths based on subtracting radius from R
    modelarr_['depth'] = R - modelarr_['radius'].pint.to(R.units)
                            
    return modelarr_

PREM = read_mineos_cards('Files/PREM750_CARDS')

PREM

/home/globalseismology/.conda/envs/fall2022-sef/lib/python3.9/site-packages/pint_pandas/pint_array.py:648: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  return np.array(qtys, dtype="object", copy=copy)
/home/globalseismology/.conda/envs/fall2022-sef/lib/python3.9/site-packages/pint_pandas/pint_array.py:648: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  return np.array(qtys, dtype="object", copy=copy)

	radius	rho	vpv	vsv	qkappa	qmu	vph	vsh	eta	depth
0	0.0	13088.48	11262.21	3667.8	1327.6	84.6	11262.21	3667.8	1.0	6371000.0
1	6824.0	13088.47	11262.2	3667.79	1327.6	84.6	11262.2	3667.79	1.0	6364176.0
2	13648.0	13088.44	11262.18	3667.78	1327.6	84.6	11262.18	3667.78	1.0	6357352.0
3	20472.0	13088.39	11262.14	3667.75	1327.6	84.6	11262.14	3667.75	1.0	6350528.0
4	27296.0	13088.32	11262.09	3667.72	1327.6	84.6	11262.09	3667.72	1.0	6343704.0
...	...	...	...	...	...	...	...	...	...	...
745	6369800.0	1020.0	1450.0	0.0	57822.5	0.0	1450.0	0.0	1.0	1200.0
746	6370100.0	1020.0	1450.0	0.0	57822.5	0.0	1450.0	0.0	1.0	900.0
747	6370400.0	1020.0	1450.0	0.0	57822.5	0.0	1450.0	0.0	1.0	600.0
748	6370700.0	1020.0	1450.0	0.0	57822.5	0.0	1450.0	0.0	1.0	300.0
749	6371000.0	1020.0	1450.0	0.0	57822.5	0.0	1450.0	0.0	1.0	0.0

750 rows × 10 columns

# grabbing values in 'g/cc' as a pandas series
PREM.rho.pint.to('g/cc')

0      13.088480000000002
1      13.088470000000003
2      13.088440000000004
3      13.088390000000002
4      13.088320000000003
              ...        
745    1.0200000000000002
746    1.0200000000000002
747    1.0200000000000002
748    1.0200000000000002
749    1.0200000000000002
Name: rho, Length: 750, dtype: pint[gram / cubic_centimeter]

# grabbing value in 'g/cc' as a numpy array
PREM.rho.pint.to('g/cc').values.data

array([13.08848, 13.08847, 13.08844, 13.08839, 13.08832, 13.08822,
       13.08811, 13.08798, 13.08783, 13.08766, 13.08746, 13.08725,
       13.08702, 13.08676, 13.08649, 13.0862 , 13.08588, 13.08555,
       13.08519, 13.08482, 13.08442, 13.08401, 13.08357, 13.08311,
       13.08264, 13.08214, 13.08162, 13.08109, 13.08053, 13.07995,
       13.07935, 13.07873, 13.07809, 13.07744, 13.07676, 13.07606,
       13.07534, 13.0746 , 13.07384, 13.07306, 13.07225, 13.07143,
       13.07059, 13.06973, 13.06885, 13.06795, 13.06702, 13.06608,
       13.06512, 13.06413, 13.06313, 13.0621 , 13.06106, 13.06   ,
       13.05891, 13.05781, 13.05668, 13.05553, 13.05437, 13.05318,
       13.05198, 13.05075, 13.0495 , 13.04823, 13.04695, 13.04564,
       13.04431, 13.04296, 13.04159, 13.0402 , 13.03879, 13.03736,
       13.03591, 13.03444, 13.03295, 13.03144, 13.02991, 13.02836,
       13.02679, 13.0252 , 13.02358, 13.02195, 13.0203 , 13.01863,
       13.01693, 13.01522, 13.01349, 13.01173, 13.00996, 13.00816,
       13.00635, 13.00451, 13.00266, 13.00078, 12.99888, 12.99697,
       12.99503, 12.99307, 12.9911 , 12.9891 , 12.98708, 12.98504,
       12.98299, 12.98091, 12.97881, 12.97669, 12.97455, 12.97239,
       12.97021, 12.96801, 12.96579, 12.96355, 12.96129, 12.95901,
       12.9567 , 12.95438, 12.95204, 12.94968, 12.94729, 12.94489,
       12.94247, 12.94002, 12.93756, 12.93508, 12.93257, 12.93005,
       12.9275 , 12.92494, 12.92235, 12.91975, 12.91712, 12.91447,
       12.91181, 12.90912, 12.90641, 12.90368, 12.90094, 12.89817,
       12.89538, 12.89257, 12.88974, 12.88689, 12.88402, 12.88113,
       12.87822, 12.87529, 12.87234, 12.86937, 12.86638, 12.86337,
       12.86034, 12.85729, 12.85421, 12.85112, 12.84801, 12.84488,
       12.84172, 12.83855, 12.83535, 12.83214, 12.82891, 12.82565,
       12.82238, 12.81908, 12.81576, 12.81243, 12.80907, 12.8057 ,
       12.8023 , 12.79888, 12.79545, 12.79199, 12.78851, 12.78501,
       12.78149, 12.77795, 12.7744 , 12.77082, 12.76722, 12.7636 ,
       12.16635, 12.15977, 12.15314, 12.14645, 12.13971, 12.13291,
       12.12605, 12.11914, 12.11218, 12.10515, 12.09807, 12.09093,
       12.08373, 12.07648, 12.06917, 12.06179, 12.05437, 12.04688,
       12.03933, 12.03172, 12.02405, 12.01633, 12.00854, 12.00069,
       11.99278, 11.98481, 11.97678, 11.96868, 11.96053, 11.95231,
       11.94403, 11.93569, 11.92728, 11.91881, 11.91028, 11.90168,
       11.89302, 11.8843 , 11.87551, 11.86666, 11.85774, 11.84875,
       11.8397 , 11.83058, 11.8214 , 11.81215, 11.80284, 11.79346,
       11.784  , 11.77449, 11.7649 , 11.75525, 11.74553, 11.73574,
       11.72588, 11.71595, 11.70595, 11.69589, 11.68575, 11.67554,
       11.66526, 11.65492, 11.6445 , 11.634  , 11.62344, 11.61281,
       11.6021 , 11.59132, 11.58047, 11.56955, 11.55855, 11.54748,
       11.53634, 11.52512, 11.51383, 11.50246, 11.49102, 11.4795 ,
       11.46791, 11.45625, 11.4445 , 11.43269, 11.42079, 11.40882,
       11.39677, 11.38464, 11.37244, 11.36016, 11.3478 , 11.33537,
       11.32285, 11.31026, 11.29758, 11.28483, 11.272  , 11.25909,
       11.2461 , 11.23303, 11.21987, 11.20664, 11.19333, 11.17993,
       11.16645, 11.15289, 11.13925, 11.12553, 11.11172, 11.09783,
       11.08386, 11.0698 , 11.05566, 11.04144, 11.02713, 11.01274,
       10.99826, 10.9837 , 10.96905, 10.95432, 10.9395 , 10.92459,
       10.9096 , 10.89452, 10.87935, 10.8641 , 10.84876, 10.83333,
       10.81781, 10.80221, 10.78651, 10.77073, 10.75486, 10.7389 ,
       10.72285, 10.70671, 10.69048, 10.67416, 10.65775, 10.64124,
       10.62465, 10.60796, 10.59119, 10.57432, 10.55736, 10.5403 ,
       10.52316, 10.50592, 10.48858, 10.47116, 10.45363, 10.43602,
       10.41831, 10.40051, 10.38261, 10.36461, 10.34652, 10.32834,
       10.31005, 10.29168, 10.2732 , 10.25463, 10.23596, 10.2172 ,
       10.19833, 10.17937, 10.16031, 10.14115, 10.1219 , 10.10254,
       10.08309, 10.06353, 10.04388, 10.02412, 10.00427,  9.98432,
        9.96426,  9.9441 ,  9.92384,  9.90348,  5.56645,  5.56175,
        5.55705,  5.55236,  5.54766,  5.54297,  5.53828,  5.53359,
        5.5289 ,  5.52421,  5.51953,  5.51485,  5.51016,  5.50548,
        5.50081,  5.49613,  5.49145,  5.49145,  5.48673,  5.48201,
        5.47729,  5.47257,  5.46785,  5.46313,  5.45842,  5.4537 ,
        5.44899,  5.44427,  5.43956,  5.43485,  5.43013,  5.42542,
        5.42071,  5.416  ,  5.41129,  5.40657,  5.40186,  5.39715,
        5.39244,  5.38773,  5.38302,  5.3783 ,  5.37359,  5.36888,
        5.36417,  5.35945,  5.35474,  5.35002,  5.34531,  5.34059,
        5.33587,  5.33116,  5.32644,  5.32172,  5.317  ,  5.31228,
        5.30755,  5.30283,  5.2981 ,  5.29338,  5.28865,  5.28392,
        5.27919,  5.27446,  5.26972,  5.26498,  5.26025,  5.25551,
        5.25077,  5.24602,  5.24128,  5.23653,  5.23178,  5.22703,
        5.22227,  5.21752,  5.21276,  5.208  ,  5.20323,  5.19847,
        5.1937 ,  5.18893,  5.18415,  5.17938,  5.1746 ,  5.16982,
        5.16503,  5.16024,  5.15545,  5.15065,  5.14586,  5.14106,
        5.13625,  5.13144,  5.12663,  5.12182,  5.117  ,  5.11217,
        5.10735,  5.10252,  5.09769,  5.09285,  5.08801,  5.08316,
        5.07831,  5.07346,  5.0686 ,  5.06374,  5.05887,  5.054  ,
        5.04913,  5.04425,  5.03936,  5.03447,  5.02958,  5.02468,
        5.01978,  5.01487,  5.00996,  5.00504,  5.00012,  4.99519,
        4.99026,  4.98532,  4.98038,  4.97543,  4.97047,  4.96551,
        4.96055,  4.95558,  4.9506 ,  4.94562,  4.94063,  4.93564,
        4.93064,  4.92563,  4.92062,  4.9156 ,  4.91058,  4.90555,
        4.90051,  4.89547,  4.89042,  4.88537,  4.88031,  4.87524,
        4.87016,  4.86508,  4.85999,  4.8549 ,  4.8498 ,  4.84469,
        4.83957,  4.83445,  4.82932,  4.82418,  4.81904,  4.81388,
        4.80873,  4.80356,  4.79839,  4.7932 ,  4.78802,  4.78282,
        4.77761,  4.7724 ,  4.76718,  4.76195,  4.75672,  4.75148,
        4.74622,  4.74096,  4.73569,  4.73042,  4.72513,  4.71984,
        4.71454,  4.70923,  4.70391,  4.69858,  4.69325,  4.6879 ,
        4.68255,  4.67719,  4.67182,  4.66644,  4.66105,  4.65565,
        4.65024,  4.64482,  4.6394 ,  4.63396,  4.62852,  4.62306,
        4.6176 ,  4.61213,  4.60664,  4.60115,  4.59565,  4.59014,
        4.58461,  4.57908,  4.57354,  4.56799,  4.56242,  4.55685,
        4.55127,  4.54568,  4.54007,  4.53446,  4.52883,  4.5232 ,
        4.51755,  4.5119 ,  4.50623,  4.50055,  4.49486,  4.48916,
        4.48345,  4.47773,  4.472  ,  4.46625,  4.4605 ,  4.45473,
        4.44895,  4.44317,  4.44317,  4.4393 ,  4.43543,  4.43156,
        4.42767,  4.42379,  4.4199 ,  4.416  ,  4.4121 ,  4.4082 ,
        4.40428,  4.40037,  4.39644,  4.39252,  4.38858,  4.38465,
        4.3807 ,  3.99214,  3.99097,  3.98981,  3.98864,  3.98748,
        3.98632,  3.98515,  3.98399,  3.98282,  3.98166,  3.98049,
        3.97933,  3.97817,  3.977  ,  3.97584,  3.97584,  3.96744,
        3.95903,  3.95063,  3.94223,  3.93383,  3.92542,  3.91702,
        3.90862,  3.90022,  3.89181,  3.88341,  3.87501,  3.86661,
        3.8582 ,  3.8498 ,  3.8414 ,  3.833  ,  3.82459,  3.81619,
        3.80779,  3.79939,  3.79098,  3.78258,  3.77418,  3.76578,
        3.75737,  3.74897,  3.74057,  3.73217,  3.72376,  3.54326,
        3.53968,  3.53609,  3.53251,  3.52893,  3.52534,  3.52176,
        3.51818,  3.51459,  3.51101,  3.50743,  3.50385,  3.50026,
        3.49668,  3.4931 ,  3.48951,  3.48593,  3.48235,  3.47877,
        3.47518,  3.4716 ,  3.46802,  3.46443,  3.46085,  3.45727,
        3.45368,  3.4501 ,  3.44652,  3.44294,  3.43935,  3.43577,
        3.35949,  3.36016,  3.36082,  3.36148,  3.36214,  3.3628 ,
        3.36346,  3.36413,  3.36479,  3.36545,  3.36611,  3.36677,
        3.36743,  3.36809,  3.36876,  3.36942,  3.37008,  3.37074,
        3.3714 ,  3.37206,  3.37273,  3.37339,  3.37405,  3.37471,
        3.37471,  3.37514,  3.37557,  3.37601,  3.37644,  3.37687,
        3.3773 ,  3.37773,  3.37817,  3.3786 ,  3.37903,  3.37946,
        3.37989,  3.38032,  3.38076,  2.9    ,  2.9    ,  2.9    ,
        2.9    ,  2.9    ,  2.9    ,  2.9    ,  2.9    ,  2.9    ,
        2.9    ,  2.9    ,  2.6    ,  2.6    ,  2.6    ,  2.6    ,
        2.6    ,  2.6    ,  2.6    ,  2.6    ,  2.6    ,  2.6    ,
        2.6    ,  1.02   ,  1.02   ,  1.02   ,  1.02   ,  1.02   ,
        1.02   ,  1.02   ,  1.02   ,  1.02   ,  1.02   ,  1.02   ])

Mini-tutorial

• Integration with scipy

SciPy is a collection of mathematical algorithms and convenience functions built on the NumPy extension of Python. SciPy allows us to optimize functions, find roots of functions, interpolate, and integrate. It also provides mathematical and physical constants

Since PREM provided the density at specific depths, and we do not know the analytical form of density, we will use the trapezoidal rule to integrate.

In this mini-tutorial, we will compute a definite integral of \(y = f(x) = 10 - x^2\) from \(x=0\) to \(x=3\). This integral can be evaluated easily the value by hand:

(4)\[\begin{equation} \int_{-3}^3 (10 - x^2) dx = \left(10x - \frac{x^3}{3}\right) \Big|_{x=-3}^{x=3} = 42 \end{equation}\]

We will use scipy.integrate.trapezoid to integrate the function with given (x,y) using Trapezoidal rule.

Focus on the following lines:

from scipy import integrate

integral = integrate.trapezoid(y_trap, x_trap)

Run the next two cells below and then adjust the slider on the figure. Notice the area approaches 42 when the number of trapezoids increases.

import numpy as np
from scipy import integrate

import plotly.graph_objects as go
from plotly.offline import init_notebook_mode, iplot
init_notebook_mode()

# Create figure
fig = go.Figure()

# actual curve
x_true = np.linspace(-4,4,801)
y_true = 10 - x_true ** 2

num_steps = 21

# True_curve
trace_list1 = []

# Outline of Trapezoids + shades
trace_list2 = []

# divides of Trapezoids
trace_list3 = []

# Add traces, one for each slider step
for number_of_trapezoids in range(1,1+num_steps):
    trace_list1 += [go.Scatter(x=x_true, y=y_true, 
                               visible=True,
                               line=dict(color="#000000", width=2),
                               mode="lines",
                               name=r"$y = 10 - x^2$")]
    
    # integrate y = f(x) from x=0 to x=3
    x_trap = np.linspace(-3,3,number_of_trapezoids+1)
    y_trap = 10 - x_trap ** 2

    # integration takes place here
    integral = integrate.trapezoid(y_trap, x_trap)\

    trace_list2 += [go.Scatter(
            visible=False,
            line=dict(color="#E00000", width=2),
            name=str(number_of_trapezoids) + " trapezoids",
            fill='tozeroy',
            x=x_trap,
            y=y_trap)]
    
    x_plot = x_trap
    y_plot = 0 * y_trap
    
    for jj in range(number_of_trapezoids+1):
        x_plot = np.concatenate((x_plot, [np.nan], [x_trap[jj], x_trap[jj]]))
        y_plot = np.concatenate((y_plot, [np.nan], [0, y_trap[jj]]))
    
    trace_list3 += [go.Scatter(
            visible=False,
            line=dict(color="#E00000", width=2),
            name="area = %.3f"%integral,
            x=x_plot,
            y=y_plot)]
    

fig = go.Figure(data=trace_list1+trace_list2+trace_list3)

# Make 10th trace visible
fig.data[10-1].visible = True
fig.data[10-1+num_steps].visible = True
fig.data[10-1+2*num_steps].visible = True

steps = []
for i in range(num_steps):
    # Hide all traces
    step = dict(
        method = 'restyle',
        args = ['visible', [False] * len(fig.data)],
        label = str(i + 1)
    )
    # Enable the two traces we want to see
    step['args'][1][i] = True
    step['args'][1][i+num_steps] = True
    step['args'][1][i+2*num_steps] = True
    
    # Add step to step list
    steps.append(step)

sliders = [dict(
    active=9,
    steps = steps,
    currentvalue={"prefix": "Number of trapezoids: "}
)]

fig.layout.sliders = sliders

iplot(fig, show_link=False)

TO DO

Question 3 Plot the density profile. Include the location of the inner core boundary (radius = 1221.5 km) and the core-mantle boundary (3480 km). Comment on relative distribution of lighter versus heavier material inside the Earth between the regions-mantle, inner and outer core. Calculate the mass of the Earth by integration. Does the calculated value of the entire Earth’s mass match the actual Earth’s mass?

Answer

Question 4 Compute the gravitational acceleration at any radii shown in PREM. Where inside the Earth is the gravitational acceleration maximum? Please explain the factors that might lead to this observation.

Answer

Question 5 Compute the pressure based on the density profile of PREM. What are the pressures at the core-mantle boundary and the center of the Earth? How do these values compare to the atmospheric pressure on the Earth’s surface?

Answer

---

Part II: Plate Motion

Below is a sketch that shows the (simplified) geometry of the lithospheric plates near the coast of Oregon and Washington. Three plates are present: the North America plate, the Pacific plate, and the Juan de Fuca plate.

Figure: Schematic of the plate configuration offshore Washington State and the Paciic Northwest.

Assume (1) that the spreading of the Juan de Fuca plate and Pacific plate is perpendicular to the ridge as drawn, and is 40 mm/yr, and (2) that the strike-slip motion between the North America plate and the Pacific plate is parallel to the North America—Pacific plate boundary and is 56 mm/yr. The orientation of the ridge is 10\(^{\circ}\) East of North. The orientation of the North America—Pacific plate boundary is 35\(^{\circ}\) West of North.

import pandas as pd
import json
import numpy as np
import plotly.express as px
import plotly.graph_objects as go

def cutoff_track(lon, lat, cutoff_lon):
    """
    Insert nan between two points that cross the cutoff lontitude. 
    Avoid horizontal line when plot the track
    
    Input Parameters:
    ----------------
    
    lon: longitude array
    lat: latitude array
    cutoff_lon: divide the curve at this cutoff longitude
    
    Return:
    ------
    
    lon, lat: updated arrays
    """
    xtrk = np.mod(lon - cutoff_lon, 360)
    is_cross = np.where(np.abs(xtrk[1:-1]-xtrk[0:-2]) > 90)[0]
    nans = np.empty(np.size(is_cross))
    nans[:] = np.nan
    lon = np.insert(lon, is_cross+1, nans)
    lat = np.insert(lat, is_cross+1, nans)
    return lon, lat

def load_plate_boundaries(filename):
    """
    Input Parameters:
    ----------------
    
    filename: input JSON file name containing the curve locations
    """    

    # Opening JSON file
    f = open(filename)

    # returns JSON object as a dictionary
    boundary = json.load(f)
    
    f.close()

    boundary = np.array(boundary[1])
    lons = boundary[:,0]
    lats = boundary[:,1]

    lons, lats = cutoff_track(lons, lats, 180)
    
    return lons, lats

def add_arrow(fig, p0, p1, label_txt="", c="black"):
    """
    Draw an arrow starting from p0 pointing at p1
    
    Input Parameters:
    -----------------
    
    fig:       target figure object
    p0:        [lon, lat] where the arrow starts
    p1:        [lon, lat] where the arrow heads
    label_txt: label text     [default: ""]
    c:         arrow color    [default: "black"]
    """
    # vector of the arrow
    v = p1-p0
    # normalized vector
    w = v/np.linalg.norm(v)
    # orthogonal on v
    u  =np.array([-v[1], v[0]])

    widh = 0.15
    l = 0.75

    P = p1-l*w
    S = P - widh*u
    T = P + widh*u
    fig.add_trace(
        go.Scattergeo(
            lon=[p0[0], p1[0], S[0], np.nan, p1[0], T[0]], 
            lat=[p0[1], p1[1], S[1], np.nan, p1[1], T[1]],
            mode='lines',
            line_color=c,
            line_width=2,
            showlegend=False,
            text=label_txt
        )
    )
    
def add_oval(fig, center, a, b, rot_angle, label_txt="", c="black"):
    """
    Draw an oval
    
    Input Parameters:
    -----------------
    
    fig:       target figure object
    center:    [lon, lat] of the oval's center
    a:         semi-major axis
    b:         semi-minor axis
    rot_angle: rotation angle counter clockwise from horizontal position (in degrees)
    label_txt: label text     [default: ""]
    c:         arrow color    [default: "black"]
    """
    rot_angle = rot_angle * np.pi / 180
    
    theta = np.linspace(0, 2*np.pi, 101)

    x_oval_horizontal = a * np.cos(theta)
    y_oval_horizontal = b * np.sin(theta)

    # rotation
    lon_oval = center[0] + x_oval_horizontal * np.cos(rot_angle) - y_oval_horizontal * np.sin(rot_angle)
    lat_oval = center[1] + x_oval_horizontal * np.sin(rot_angle) + y_oval_horizontal * np.cos(rot_angle)
    
    fig.add_trace(
        go.Scattergeo(
            lon=lon_oval, 
            lat=lat_oval,
            mode='lines',
            line_color=c,
            line_width=4,
            showlegend=False,
            text=label_txt
        )
    )


def plot_earthquake_map(catalog):
    """
    Input Parameters:
    ----------------
    
    catalog: dataframe containing earthquakes
    """   
    
    rlons, rlats = load_plate_boundaries('../PS_Plate_Tectonics/Files/ridge.json')
    tlons, tlats = load_plate_boundaries('../PS_Plate_Tectonics/Files/transform.json')
    hlons, hlats = load_plate_boundaries('../PS_Plate_Tectonics/Files/trench.json')

    fig1a = px.line_geo(lon=rlons, lat=rlats, width=1200, height=600)
    fig1a.update_traces(line={'color':'skyblue', 'width':3})
    
    fig1b = px.line_geo(lon=tlons, lat=tlats, width=1200, height=600)
    fig1b.update_traces(line={'color':'skyblue', 'width':3})
    
    fig1c = px.line_geo(lon=hlons, lat=hlats, width=1200, height=600)
    fig1c.update_traces(line={'color':'skyblue', 'width':3})
    
    fig2 = px.scatter_geo(catalog, lon="ep.lon", lat="ep.lat", size="centroid.MW", color="ep.depth", width=600, height=600, 
                     labels={"ep.lon":"longitude (degrees)", "ep.lat":"latitude (degrees)", "centroid.MW":"Moment Magnitude",
                             "ep.depth":"depth (km)", "ep.timestamp":"origin time"}, 
                     hover_data=['ep.timestamp'], 
                     title=" ")
    
    fig = go.Figure(data=fig2.data + fig1a.data + fig1b.data + fig1c.data, layout = fig2.layout)
    fig.update_geos(lataxis_showgrid=True, lonaxis_showgrid=True)
    fig.update_layout(
        geo = dict(
            lonaxis = dict(
                showgrid = True,
                gridwidth = 0.5,
                range= [ -138.0, -118.0 ],
                dtick = 5
            ),
            lataxis = dict (
                showgrid = True,
                gridwidth = 0.5,
                range= [ 37.0, 57.0 ],
                dtick = 5
            )
        )
    )
    
    # Adds arrows
    length_scale = 25
    
    # Juan de Fuca ridge: Juan de Fucca arrow
    p0 = np.array([-129.5, 46.4])
    length = 40 / length_scale
    direction = 100 * np.pi/180
    p1 = p0 + length * np.array([np.sin(direction), np.cos(direction)])
    add_arrow(fig, p0, p1, "40 mm/yr")
    
    # Juan de Fuca ridge: Pacific plate arrow
    p0 = np.array([-130.5, 46.6])
    length = 40 / length_scale
    direction = 280 * np.pi/180
    p1 = p0 + length * np.array([np.sin(direction), np.cos(direction)])
    add_arrow(fig, p0, p1, "40 mm/yr")
    
    # North American plate arrow
    p0 = np.array([-134, 57])
    length = 56 / length_scale
    direction = 145 * np.pi/180
    p1 = p0 + length * np.array([np.sin(direction), np.cos(direction)])
    add_arrow(fig, p0, p1, "56 mm/yr")
    
    # North American plate arrow
    p0 = np.array([-135, 53])
    length = 56 / length_scale
    direction = 325 * np.pi/180
    p1 = p0 + length * np.array([np.sin(direction), np.cos(direction)])
    add_arrow(fig, p0, p1, "56 mm/yr")
    
    # Juan de Fuca ridge: Juan de Fuca plate arrow
    p0 = np.array([-127, 41.4])
    length = 40 / length_scale
    direction = 100 * np.pi/180
    p1 = p0 + length * np.array([np.sin(direction), np.cos(direction)])
    add_arrow(fig, p0, p1, "40 mm/yr", "white")
    
    # Juan de Fuca ridge: Pacific plate arrow
    p0 = np.array([-128, 41.6])
    length = 40 / length_scale
    direction = 280 * np.pi/180
    p1 = p0 + length * np.array([np.sin(direction), np.cos(direction)])
    add_arrow(fig, p0, p1, "40 mm/yr")
    
    # Add oval
    add_oval(fig, [-128.7, 43.7], 2.5, 1.2, -20, "", "limegreen")
    fig.add_annotation(x=.25, y=.40,
            text="???",
            font=dict(
                color="darkgreen",
                size=20
            ),
            showarrow=False)
    
    fig.show()

# Read earthquake catalog
cmt_cat = pd.read_table("../PS_Plate_Tectonics/Files/Global_CMT_1972_2018.csv", delimiter=',')
cmt_cat = cmt_cat.sort_values(by=['ep.depth', 'centroid.MW'])

/tmp/ipykernel_129260/2244039337.py:2: DtypeWarning:

Columns (15) have mixed types. Specify dtype option on import or set low_memory=False.

Review

Feel free to adjust the filter parameters in the cell below to have more or less earthquakes on the map. You may hover the dots to get the information of the earthquakes. The colors of the lines indicate the type of plate boundaries: Red is a Trench, Green is a Transform, and Blue is a Ridge boundary.

# Filter to improve readablility
Mw_min    = 5                                 # Moment magnitude
Mw_max    = 10                                # Moment magnitude
depth_min = 0                                 # minimum depth in km
depth_max = 6400                              # maximum depth in km
datetime_min = '1000-01-01 00:00:00.000000'   # 'yyyy-mm-dd hh:mm:ss.SSSSSS'
datetime_max = '2050-01-01 00:00:00.000000'   # 'yyyy-mm-dd hh:mm:ss.SSSSSS'
lon_min = -138
lon_max = -118
lat_min = 37
lat_max = 57

# Don't change these two statements
mask = (cmt_cat['centroid.MW'] >= Mw_min) & (cmt_cat['centroid.MW'] <= Mw_max) & (cmt_cat['ep.depth'] >= depth_min) & (cmt_cat['ep.depth'] <= depth_max) & (cmt_cat['ep.timestamp'] >= datetime_min) & (cmt_cat['ep.timestamp'] <= datetime_max)
mask = mask & (cmt_cat['ep.lon'] >= lon_min) & (cmt_cat['ep.lon'] <= lon_max) & (cmt_cat['ep.lat'] >= lat_min) & (cmt_cat['ep.lat'] <= lat_max)
plot_earthquake_map(cmt_cat[mask])

/tmp/ipykernel_129260/2372063229.py:18: RuntimeWarning:

invalid value encountered in remainder

Mini-tutorial

• Scalars and Vectors in Python

scalars are quantities having only a magnitude (e.g. a temperature of 300 K). It is nonsense to ask for direction of a scalar. On the other hand, vectors are quantities having a magnitude as well as a direction (e.g. a velocity with a magnitude of 10 m/s and a direction of 30\(^{\circ}\) East of North). Vectors can be represent as an array of numbers specifying the quantity for each direction. For example, a velocity with a magnitude of 10 m/s and a direction of 30\(^{\circ}\) East of North can be represented as \(v = \begin{pmatrix} 5 \\ 5\sqrt{3} \end{pmatrix}\) m/s where the first number is the velocity in the east direction and the second number is the velocity in the north direction. This is similar to “vector” object in programming language. In Python, this means you can use numpy.array to represent vectors. The conversion of the magnitude and the (azimuthal) direction to the values in east and north direction can be done as follow:

import numpy as np

magnitude = 10
direction = 30 # in degrees

direction_radian = direction * np.pi / 180   # convert to radians

# vector
v_east  = magnitude * np.sin(direction_radian)
v_north = magnitude * np.cos(direction_radian)
v = np.array([v_east, v_north])

# display a vector
v

array([5.        , 8.66025404])

Figure: Horizontal and Vertical components of a vector.

You may add two vectors by adding up element-wise from the two vectors. For numpy.array, you may simply add two vectors together like w = u + v where u and v are two vectors.

# vector addition
u = np.array([3, -2])
v = np.array([2,  1])

# the result of vector addition
w = u + v

# display the result vector
w

array([ 5, -1])

Figure: Example of an addition of two vectors.

TO DO

Question 6 What is the motion (the direction with respect to North in degrees, and the magnitude in mm/yr) of the Juan de Fuca plate with respect to the North America plate? Draw an arrow showing this direction on the map.

Hint: Because the Juan de Fuca plate is small, you can consider this a ‘flat Earth’ problem. Also, recall that if three plates interact, and the velocity vector of plate \(A\) with respect to plate \(B\) is \(_B\textbf{V}_A\), and that of plate \(B\) with respect to plate \(C\) is \(_C\textbf{V}_B\), then \(_C\textbf{V}_A= _C\textbf{V}_B + _B\textbf{V}_A\)

Answer

Question 7.1 Based on the relative plate motion, what type of plate boundary is this? (Reminder: there are three types of plate boundaries: convergent or collision boundaries, divergent or spreading boundaries, and transform boundaries.)

Answer

Question 7.2 Based on plate motions, figure out what kind of boundary is denoted by the oval in the interactive figure above. What type of fault (strike-slip, normal, or reverse) is giving rise to these earthquakes?

Answer

---

Part III: Quake Magnitude

Seismologists were for a long time content to use the intensity scale and the magnitude scale to measure and compare the sizes of earthquakes. Both of these scales are empirical, even though determining a magnitude involves making a measurement on a seismogram. Neither measure depends directly on the size of the earthquake fault or the amount of slip that occurred on it. These scales can be contrasted against the seismic moment (or scalar moment), usually denoted by \(M_0\). This quantity can be routinely and accurately determined from seismograms. To do so is more complicated than measuring a magnitude, but the advantage of the seismic moment is that it is directly related to physical quantities. The seismic moment is

\[ M_0 =\mu \times A \times u, \]

where \(A\) is the area of the fault that slipped, \(u\) is the average amount of slip on the fault, and \(\mu\) is the shear modulus, which is a property of the rocks surrounding the fault. The shear modulus describes how much force is necessary to elastically bend the rock – a typical value for rocks in the Earth’s crust is \(\mu= 3.3 \times 10^{10} \text{N}/\text{m}^2\).

Figure: Comparison of moment, magnitudes, fault area, and fault slip for four earthquakes. Sizes of rectangles indicate the areas of the faults that slip. \(M_s\), the magnitude obtained from amplitudes of surface-wave vibrations, saturates for events with \(M_w > 8\) and so is no longer a useful measurement of earthquake size (Modified from Figure 4.6-3 of Stein and Wysession, Introduction to Seismology, Earthquakes and Earth Structure).

We can calculate the seismic moment of the greatest known earthquake, the 1960 Chile earthquake, by using estimates of the fault dimensions and slip: the length of the fault \(L∼1000\) km, the width \(W∼200\) km and the total slip \(u∼20\) meters: \(M_0 = 1000 \times 10^3 \times 200 \times 10^3 \times 20 \times 3 \times 10^{10} = 1.2 \times 10^{23} \text{Nm}\) (Newton \(\times\) meter). Physically this is not an energy but a moment, i.e. a force times a length of lever arm.

Figures: (top) Histogram of earthquakes of various moment magnitudes, (middle) fraction of the cumulative moment released by various earthquakes, (bottom) cumulative moment with time since modern recordings of earthquakes. Data Source: Global CMT Project

Take away message: Most earthquakes are small, so the total moment (or energy) released is dominated by a small number of large earthquakes.

TO DO

Question 8 The seismic moment of the 1994 Northridge (L.A.) earthquake was \(2 \times 10^{19}\) Nm. The slip appears to have been 1.3 meters. How big was the fault area that slipped?

Answer

Question 9 If the 20 meters of slip that occurred in the 1960 Chile earthquake corresponded to the accumulated plate motions between the Nazca and South America plates since the previous great earthquake in this region — when would that earthquake have occurred? [Useful fact: Approximately 8 km of the Nazca plate has subducted beneath South America in the last 100,000 years.]

Answer

Question 10 When an earthquake happens, the energy stored in the elastically deformed plates is converted into other energy forms, primarily frictional heat and elastic wave energy (most of which turns into heat, too). It turns out that it is very difficult to measure the amount of energy \(E\) released in an earthquake accurately but, to within a factor of 5, it appears to be proportional to the seismic moment. A relationship which seems to be on average correct is \(E= M_0 \times 5 \times 10^{-5}\). The energy \(E\) is measured in newton-meters (Nm) or joules (J)

Estimate the energy released in the 1994 Northridge (L.A.) earthquake. Compare the value you get with U.S. total daily energy consumption (\(~5 \times 10^{17}\) J).

Answer

Question 11 Instead of getting rid of the magnitude scale, seismologists decided to create a new magnitude based on the seismic moment. This was done so as not to confuse the public (and themselves), to avoid dealing with very large numbers like \(10^{20}\), and to be able to compare old magnitude estimates with modern seismic moments, and so on. From the seismic moment (\(M_0\) measured in Nm) one can calculate a magnitude \(M_W\) called the moment magnitude using the equation

\[ M_W = \frac{2}{3} (\log_{10}M_0 - 9.05) \]

The constants in this equation were chosen so that the new scale coincides with Richter’s traditional magnitude scale. That is, an earthquake that is measured with traditional means to have a Richter magnitude of 6.5 will typically also have a seismic moment that corresponds to a moment magnitude of 6.5. However, for an individual earthquake, it is not uncommon for the the traditional magnitude to be wrong by half a magnitude unit. It is particularly difficult to measure the magnitudes of very large earthquakes using the ‘old-fashioned’ approach.

What is the moment magnitude of the 1960 Chile earthquake?

Answer

Question 12 From these relationships, calculate how much more energy is released in a \(M_W=7\) than in a \(M_W=6\) earthquake (Note: Express this as a fraction or a ratio).

Answer

Question 13 Why do we tend to have large earthquakes (\(M_W > 7.5\)) in the subduction zones? Think about the factors that contibute to the scalar moment based on the equations above.

Answer