Calculate the vertical displacement using the Mogi model#
In this notebook, we will define the Mogi model to measure vertical displacements and apply it on data collected at Sakurajima volcano following the 1914 eruption. The following python modules are required in order to run the script:
matplotlib
numpy
math
Source: https://gscommunitycodes.usf.edu/geoscicommunitycodes/public/numeracy/numeracy_mogi/mogi.html#python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import math
Define Mogi Model for vertical displacement#
Refresher for Python arithmetic operators: https://www.w3schools.com/python/gloss_python_arithmetic_operators.asp
def mogi_vertical(x,a,d,g,delta_p,sigma):
# Calculate the vertical displacement (m)
# at the surface due to a subsurface
# pressure source using Mogi's solution
# x is horizontal distance (m), radial from source
# a is source radius (m) must be << d
# d is depth (m)
# g is shear modulus (MPa), typically 3e5 or less
# delta_p is the excess pressure of the source (MPa)
# sigma is Poisson's ratio (typically 0.25)
# w = [ what is the equation? ]
w = ??????
return w # displacement (m)
check_solution = mogi_vertical(5000,500,2500,2.5e5,100,0.25)
print(check_solution)
# does your value equal to 0.0005366563145999495 ?
Survey data from Sakurajima Volcano#
The survey data were collected using precise line levelling on Sakurajima volcano by F. Omori and were reported in Mogi (1958). The displacements are found by calculating the differences in elevations at permanent survey markers on Sakurajima Volcano between surveys in 1893 and 1914, following the 1914 eruption.
# The data is stored as a list in the following format:
# 1st column: radial distance of permanent survey markers from the center of summit in km
# 2nd column: differences in elevations in mm (i.e., vertical displacements)
# between the two surveys at the each corresponding marker
sakura = [[17.0, -291.0],
[16.2, -291.0],
[15.75, -302.0],
[13.5, -407.0],
[12.45, -446.0],
[10.4, -527.0],
[8.7, -658.0],
[7.85, -770.0],
[6.7, -894.0],
[7.75, -776.0],
[8.25, -703.0],
[9.15, -608.0],
[9.85, -526.0],
[9.7, -536.0],
[9.55, -560.0],
[9.3, -577.0],
[8.5, -684.0],
[8.55, -690.0],
[8.95, -681.0],
[9.15, -685.0],
[10.05, -616.0],
[11.1, -481.0],
[11.9, -427.0],
[12.15, -369.0],
[12.1, -348.0],
[12.5, -301.0],
[12.3, -306.0],
[10.95, -433.0],
[11.05, -426.0],
[11.0, -453.0],
[12.0, -356.0],
[12.4, -344.0],
[12.5, -363.0],
[11.6, -401.0],
[11.2, -468.0],
[12.6, -379.0],
[13.95, -314.0],
[15.8, -242.0],
[12.15, -369.0],
[12.5, -338.0],
[13.75, -260.0],
[14.15, -244.0],
[15.4, -206.0],
[17.25, -155.0],
[18.8, -143.0]]
# Rearrange data to 2 separate arrays
marker_radial_dist = np.array([item[0] for item in sakura])
obs_disp = np.array([item[1] for item in sakura])
# Plot observation
plt.plot(marker_radial_dist, obs_disp, "ko")
plt.xlabel("Radial distance from summit (km)")
plt.ylabel("Vertical displacement (mm)")
plt.title("Vertical displacements at Sakurajima after the 1914 eruption")
plt.show()
Forward modelling#
Take a guess on the source depth and excess pressure and see how the Mogi model fits on Sakurajima volcano.
## Input for Mogi Model
# variable
d = ??????? #source depth (m)
delta_p = ??????? #excess pressure (MPa)
# fixed
a = 2500.0 #source radius (m)
sigma = 0.25 #Poisson ratio
g = 30000.0 #shear modulus (MPa)
# ===================================== #
# calculate vertical displacement for this range of radial distances (m)
# !! note the units (e.g., m, mm, km, ...)
x1 = np.arange(0,20000,100) # in meters
# Calculate vertical displacement using Mogi model
# !! multiply by 1000 as Mogi model uses m and observed measurements are in mm.
syn_disp_in_meter = mogi_vertical(x1,a,d,g,delta_p,sigma)
syn_disp = syn_disp_in_meter * 1000 # convert meter to millimeter
# Compare the Mogi solution to observations for Salurajima after the 1914 eruption
# !! note the units when plotting
x1_in_km = x1/1000
# Mogi solution in red line
plt.plot(x1_in_km, syn_disp,'r')
# Observed displacements in black dots
plt.plot(marker_radial_dist, obs_disp, "ko")
plt.xlabel("radial distance from summit (km)")
plt.ylabel("vertical displacement (mm)")
plt.title("Vertical displacements at Sakurajima after the 1914 eruption")
plt.show()
# uncomment the next line if you want to save the figure:
#plt.savefig("sakurajima.png")
Finding the ‘best’ solution#
How do we know if our model, and if our chosen model parameters is any good?
One way to judge is by measuring the root mean square error (RMSE) between the predicted displacement from our model and the actual observed displacements.
The smaller RMSE, the smaller the difference between the predicted and observed values, an indicator of a good-fitting model.
# calculate the predicted displacements at the survey markers from the previously chosen parameters
y_predicted = mogi_vertical(marker_radial_dist*1000,a,d,g,delta_p,sigma)*1000
y_actual = obs_disp
# Calcuate the Mean Square Error and take the square root to obtain Root Mean Square Error
MSE = np.square(np.subtract(y_actual,y_predicted)).mean()
RMSE = math.sqrt(MSE)
print("The RMSE for Mogi Model with source depth", d, "[m] and excess pressure of", delta_p,"[MPa] is ", RMSE)
Grid-search for best-fitting Mogi model#
To find a best-fitting Mogi model, we will test a range of source depths and excess pressures and estimate their corresponding RMSE. Instead of changing the value manually one by one, we will use “nested loops” to iterate through the values.
# set a range of depth and excess pressure
depth_range = np.arange(2500,30000,500) # 2500 m to 30000 m in 500 m increment. cannot be shallower than source radius
delta_p_range = np.arange(-1000,0,50) # -1000 MPa to 0 MPa in 50 MPa interval
# fix other paramters
a = 2500.0 #source radius (m)
sigma = 0.25 #Poisson ratio
g = 30000.0 #shear modulus (MPa)
# ===================================== #
# pre-allocate RMSE array
RMSE_array = np.zeros((len(depth_range),len(delta_p_range)))
# loop through depth and excess pressure
for index0,d in enumerate(depth_range):
for index1,delta_p in enumerate(delta_p_range):
y_predicted_tmp = mogi_vertical(marker_radial_dist*1000,a,d,g,delta_p,sigma)*1000
MSE_tmp = np.square(np.subtract(y_actual,y_predicted_tmp)).mean()
RMSE_array[index0,index1] = math.sqrt(MSE_tmp)
# plot results
extent = np.min(delta_p_range), np.max(delta_p_range),np.min(depth_range), np.max(depth_range)
plt.contourf(RMSE_array,50,cmap=plt.cm.plasma,origin='lower',extent=extent)
plt.xlabel("Excess pressure [MPa]")
plt.ylabel("Source depth [m]")
cbar = plt.colorbar()
cbar.ax.set_ylabel('RMSE', rotation=90)
# find best-fitting model with the smallest RMSE
index = np.where(RMSE_array == np.amin(RMSE_array))
plt.plot(delta_p_range[index[1]], depth_range[index[0]], 'ro')
print("The best-fitting Mogi Model has a source depth of", float(depth_range[index[0]]),"[m] and excess pressure of", float(delta_p_range[index[1]]),"[MPa]")