Deformation due to a buried point-like pressure source

Visualizing Mogi's model with p5

The following sketch illustrates Mogi's calculation for a sub-surface point-like pressure source. The sphere (large white circle in the illustration) represents the position of the point-like pressure source in the sub-surface. The mouse is used to change the position of the point-like source and the calculation is updated for the position. Adjust the excess pressure, $\Delta P$, using the slider bar. The deformation observed at the surface is shown in terms of vertical displacement ($w$, mm) (white small circles) and horizontal displacement ($u$, yellow small circles). Note that the sign of the horizontal displacement changes across the point-like source. To the left of the source the horizontal displacement is in the negative $x$-direction and to the right of the point-like source the surface displacement is in the positive $x$-direction. In other words, the magnitude of surface displacement is radially symmetric with respect to the point-like source.

For this calculation, The radius of the point-like source is 0.5 km, the shear modulus, $G = 3 \times 10^4$ MPa (30 GPa), and Poisson's ratio is $\sigma = 0.25$.

[Top]

Model calculations

The following image describes the geometry of the problem. The point-like source, centered on point $O$, has radius $a$ and is located at depth $d$, where $a << d$. The displacement is calculated at point $S$, located at the surface a radial distance $x$ from $O$. The point-like pressure source has excess pressure, $\Delta P$, which results in displacement at the Earth's surface (blue shaded line), which is described by the two vectors $w$ and $u$.

The point $S$ is a distance from the center of the pressure source, $O$: $$r = \sqrt{x^2 + d^2}.$$ $x$ is the radial distance and $d$ is the depth (positive downward).

The surface displacements produced by change in excess pressure within a spherical cavity in two dimensions are given by: $$\begin{pmatrix}u\\w\end{pmatrix} = a^3\Delta P\frac{1-v}{G(x^2+d^2)^\frac{3}{2}} \begin{pmatrix}x\\d\end{pmatrix}$$ where:

• $u$ and $w$ are the horizontal $(x)$ and vertical $(z)$ components of surface displacement (m) at point $(x,0,0)$,
• $d$ is the depth (m) to the center of the cavity,
• $a$ is the radius (m) of the spherical cavity $(a << d)$,
• $\Delta P$ is the pressure change (MPa) in the cavity,
• $v$ is Poisson's ratio of the half-space, and
• $G$ is the shear modulus (MPa) of the half-space.

Notice that the displacement is directly proportional to the excess pressure, $\Delta P$. Changing $\Delta P$ changes the amplitude of the displacement. The displacement is inversely proportional to the shear modulus, $G$. The ratio $\frac{\Delta P}{G}$ controls the amplitude of the displacement. The surface displacement is directly proportional to the original volume of the source ($a^3$). The displacement is nonlinearly related to horizontal distance and depth, $x$ and $d$, respectively.

In an elastic medium, the change of volume of the source can be approximated as $\pi a^3 \Delta P/G$. Mogi’s solution is often expressed in terms of volume change to avoid the trade-off between pressure and original volume of the source. Change in volume is related to the pressure acting over the surface area of the sphere. The deformation expected from an equivalent change in volume is: $$\begin{pmatrix}u\\w\end{pmatrix} = \frac{\Delta V (1-v)}{\pi (x^2+d^2)^\frac{3}{2}} \begin{pmatrix}x\\d\end{pmatrix}$$

[Top]

Model assumptions

This model assumes that:

• The Earth's crust is an homogeneous, semi-infinite elastic half-space
• The origin of the pressure anomaly in the subsurface is a small spherical source. Specifically, the model is strictly correct for $a \le \frac{1}{3} d $.

Being a half space the surface does not have any topography, and being homogeneous and isotropic, the model does not account for possible variation of material properties (e.g. lithology or temperature). The model is purely elastic not accounting for permanent deformation.

[Top]

p5.js script

This p5.js script creates the animation shown above. The script calculates the displacement due to the point source shown in the figure above and plots the displacement values along the surface line.

var x = 100;
var y = 100;
var radius =  20; //in screen units, 1 km = 40 screen units
var x_zero = 50;
var x_scale = 40; //1 km = 40 screen units
var initial_depth = 2; //initial depth in km
var depth_zero = 300; //where the depth axis starts in screen units
var depth_scale = 40; //keep vert exaggeration 1:1
var deformation_min = 270;  //less than depth_zero 
var deformation_scale = 12;
var deformation_zero;
var del_pressure;
var slider;
     
function setup() {
  var canvas = createCanvas(600,650);
 
  // Move the canvas so it's inside our 
.
  canvas.parent('sketch-holder');
//  canvas.noStroke(); 
  slider = createSlider(10, 40, 30);
  slider.parent('slider');
  slider.position(width-190, 30);
  slider.style('width','180px','background','white');    
  fill(0)
   .strokeWeight(.25);
  del_pressure = 10;
}

function draw() {
  background('#006747');
 
 var mogi = [];
 var x1km, y1km, r;
 var mogi_x, mogi_z;
 var offset;
 var x1,  z1, m, n;
 var del_pressure = slider.value();
 
  // lerp() calculates a number between two numbers at a specific increment. 
  // The amt parameter is the amount to interpolate between the two values 
  // where 0.0 equal to the first point, 0.1 is very near the first point, 0.5 
  // is half-way in between, etc.  
  
  // Here we are moving 5% of the way to the mouse location each frame
  //if (mouseX > 0 && mouseY > 0){
  x = lerp(x,mouseX,0.05);
  y = lerp(y,mouseY,0.05);//}
  
  if (x < x_zero+radius) {x = x_zero+radius;}
  if (x > width - radius) {x = width - radius;} 
  if (y < depth_zero+2*radius) {y = depth_zero+2*radius;} 
  if (y > height - radius) {y = height - radius;}
  ellipse(x, y, 2*radius, 2*radius);
  
  ////calculate and draw deformation
  for (offset=x_zero; offset < width; offset+=10){ //count in screen units
    
     x1 = offset - (x);
     z1 = (y); 
     x1km = x1/x_scale*1000.0;
     y1km = (z1 - depth_zero)/depth_scale*1000.0 +  initial_depth*1000.0;
     r = radius/x_scale*1000.0; 
      
     mogi =gsphere(r, x1km, y1km, del_pressure);
     // convert deformation to integer
     mogi_x = mogi[0]*1000*deformation_scale;
     mogi_z = mogi[1]*1000*deformation_scale
     m = deformation_zero - mogi_z;
     n = deformation_zero - mogi_x;
     // plot point
     ellipse(offset, m, 6);
    fill('yellow');
     ellipse(offset, n, 6);
    fill(255)
  }
  //draw axes and labels
  drawAxes(x_zero, x_scale, depth_zero, depth_scale, deformation_min, deformation_scale, initial_depth);
}


function drawAxes(x_zero, x_scale, depth_zero, depth_scale, deformation_min, deformation_scale, d_i) {
 fill('#CFC493');
 textAlign(RIGHT);
 for (var i = 0; i< 9 ; i++) {
  var  d= i*depth_scale + depth_zero;
  line(x_zero,d,x_zero-10,d);
  text(i + d_i,x_zero-12,d+5);
 }
 
  //draw and label u,w axis
 textAlign(RIGHT);
 for (var i = 0; i< 11 ; i++) {
  var d= deformation_min - 2*i*deformation_scale;
  line(x_zero,d,x_zero-10,d);
  text(i*2 - 6,x_zero-12,d+5);
     if(i*2 - 6 == 0){deformation_zero = d;}
 }
  //draw label the horizontal axis
  textAlign(CENTER);
 for (var i = 0; i<14; i++) {
  var x=i*x_scale + x_zero;
  line(x,depth_zero-35,x,depth_zero-25);
  line(x,depth_zero-5,x,depth_zero+5);
  var x_value =  i;
   text(x_value,x,depth_zero-10);
 }
  // rotate(3.14159/2.0);
   textAlign(CENTER);
   text(slider.value(), width - 200, 40);
   text("Delta P (MPa)", width - 50, 20)
    
   text("z (km)", 100, 350);
   text("x (km)", 500, depth_zero+25);
   
   text("w (mm)", 100, 100);
  fill('yellow')
   text("u (mm)", 100, 120);
   
   stroke(255);
   fill(255);
   //depth arrow
   line(75, 360, 75, 400);
   triangle(70,400,80,400,75,425);
    //horiontal arrow
   line(520, depth_zero+20, 555, depth_zero+20);
   triangle(555,depth_zero+15,555,depth_zero+25,580,depth_zero+20);
    //mgal arrow
   line(75,80, 75, 40);
   triangle(70,40,80,40,75,15);
}


function gsphere(a, x2, z2, del_pressure){

  var G;
  var u;
  var w;
    
/*calculate the surface
deformation due to a point pressure source

u=  3.0/4.0 * del_pressure * r * a**3 / (G * (r**2 + d**2)**(3/2));
w=  3.0/4.0 * del_pressure * d * a**3 / (G * (r**2 + d**2)**(3/2));
distance units in meters
returns meters

a = radius of sphere;
x1 is horizontal offset
z1 is depth; 
*/
//define consts.
G=30000; //shear modulus in MPa //30 GPa
u=  3.0/4.0 * del_pressure * x2 * pow(a, 3) / (G * pow(pow(x2,2) + pow(z2, 2), 3/2));
w=  3.0/4.0 * del_pressure * z2 * pow(a, 3) / (G * pow(pow(x2,2) + pow(z2, 2), 3/2));   

return [u,w]; //return x, z displacement in meters
}

[Top]

python script

The following python code calculates the vertical displacement using the Mogi solution as a function. The calculated displacement is compared with the data collected on deflation of Sakurajima volcano following the 1914 eruption. The survey data were collected using precise line levelling. The displacements are found by calculating the differences in elevations at permanent survey markers between surveys in 1893 and 1914, following the eruption. The data indicate there was dramatic deflation of the volcano edifice following the eruption in June, 1914, likely due to depressurization of a magma reservoir. The survey data were collected by F. Omori and were reported in Mogi (1958).

The script includes python functions for vertical and horizontal displacements. It is necessary to have python modules installed to run the script. The modules required are:

• matplotlib
• numpy

The python script can be cut and paste into a jupyter notebook, or run using any python version 3 environment.

import numpy as np
import matplotlib.pyplot as plt

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 = (1-sigma)*a**3*d*delta_p/(g*(x**2+d**2)**1.5)
    return w  # m displacement

def mogi_horizontal(x,a,d,g,delta_p,sigma):
    # calculate the horizontal displacement (m) at the
    # 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)
    u = (1-sigma)*a**3*x*delta_p/(g*(x**2+d**2)**1.5)
    return u # m displacement

a = 2500.0 #source radius (m)
d = 10000.0 #source depth (m)
sigma = 0.25 #Poisson ratio
delta_p = -400.0 #excess pressure (MPa)
g = 30000.0 #rigidity of elastic half-space (MPa)

# Sakurajima displacements following 1914 eruption
# The data are taken from Mogi (1958)
# The data were collected by line levellin by
# F. Omori and colleagues
# displacement is the difference between 
# measured elevations in 1893 and 1914.
sakura = [(17.0,-291),
(16.2,-291),
(15.75,-302),
(13.5,-407),
(12.45,-446),
(10.4,-527),
(8.7,-658),
(7.85,-770),
(6.7,-894),
(7.75,-776),
(8.25,-703),
(9.15,-608),
(9.85,-526),
(9.7,-536),
(9.55,-560),
(9.3,-577),
(8.5,-684),
(8.55,-690),
(8.95,-681),
(9.15,-685),
(10.05,-616),
(11.1,-481),
(11.9,-427),
(12.15,-369),
(12.1,-348),
(12.5,-301),
(12.3,-306),
(10.95,-433),
(11.05,-426),
(11.0,-453),
(12.0,-356),
(12.4,-344),
(12.5,-363),
(11.6,-401),
(11.2,-468),
(12.6,-379),
(13.95,-314),
(15.8,-242),
(12.15,-369),
(12.5,-338),
(13.75,-260),
(14.15,-244),
(15.4,-206),
(17.25,-155),
(18.8,-143)]

# calculate vertical displacement for this
# range of radial distances (m)
x1 = np.arange(0,20000,100)

# plot vertical displacement calculated with
# Mogi solution in mm vs. km
plt.plot(x1/1000, mogi_vertical(x1,a,d,g,delta_p,sigma)*1000)

# Compare the Mogi solution to observations
# for Salurajima following deflation after
# the 1914 eruption
radial, vert_displace = zip(*sakura)
plt.plot(radial, vert_displace, "bo")
plt.xlabel("radial distance (km)")
plt.ylabel("vertical displacement (mm)")
plt.title("Mogi solution for 1914 Sakurajima deflation")
plt.savefig("sakurajima.png")
plt.show()

[Top]

python script output

The python script outputs this figure, showing the calculated vertical displacement (blue line) as a function of distance from the center of deflation. The displacement data collected by Omori are shown as solid blue circles. Changing the model parameters in the python script will change the fit between the Mogi model and Omori's observations.

[Top]

Volume reduction in the subsurface following the 1914 eruption

The Mogi model provides a means of estimating the volume of material "missing" from the subsurface following the 1914 eruption and comparison with the volume of magma erupted. The Mogi model is fit to the observed displacements (see the python script and graph), then these parameters are used to calculate the change in radial distance from the volcano as a function of vertical displacement. This function defines a solid of revolution, so integration over the range from zero displacement to the maximum displacement (at $x = 0$) gives the volume of the solid, in this case the volume of material missing as reflected in the observed radially symmetric ground deformation.

First, create a lumped parameter to simplify the algebra: $$\zeta = \frac{a^3 |\Delta P| (1-v)}{G}$$ Then, recast Mogi's equation as the radial distance in terms of magnitude of vertical displacement, $|w|$: $$x = \sqrt{\left [ \frac{\zeta d}{|w|}\right ]^{2/3} - d^2}$$ Now, find the volume of the solid of revolution, rotating about the vertical (displacement) axis: $$V = \pi \int_0^{|w_{max}|} \left [\left [ \frac{\zeta d}{|w|}\right ]^{2/3} - d^2 \right ] d|w|$$ where $|w_{max}|$ is the maximum amplitude of vertical displacement at $x = 0$. Solving: $$V = 3 \pi \zeta^{2/3} d^{2/3} |w_{max}|^{1/3} - \pi d^2 |w_{max}|$$ For displacements following the 1914 Sakurajima eruption, the model (python script) uses depth, $d = 10000$ m, excess pressure, $\Delta P = -400$ MPa, rock shear strength, $G = 30000$ MPa, pressure source radius, $a = 2500$ m, and Poisson's ratio, $v = 0.25$. The maximum vertical displacement calculated for these parameters at $x=0$ m is $|w_{max}| = 1.5625$ m. These parameters yield volume, $V = 0.98$ km$^3$. This is the volume missing, represented by the deflation of the ground around the volcano following the eruption. The volume missing due to deflation is not quite the same as the volume of the magma chamber, because the surrounding rock is compressible ($v < 0.5$). Delaney and McTigue (1994) showed that for a point-like pressure source: $$\frac{\Delta V_{surface}}{\Delta V_{chamber}} = 2(1-v) $$ In the above calculations, it is assumed that Poisson's ratio is $v = 0.25$. Given $\Delta V_{surface} = 0.98$ km$^3$, $\Delta V_{chamber} = 0.65$ km$^3$. The volume of lava reported to have erupted in the months to years following the January, 1914 eruption is reported to be 0.5 to 1.5 km$^3$, with a much smaller volume of tephra reported: $0.09 \pm 0.03$ km$^3$ (Ishihara et al., 1981, Todde et al., 2017). It appears that most of the eruption volume was emitted by the time of Omori's June, 1914 survey based on historical accounts, as the isthmus connecting Sakurajima island with the mainland formed by lava flow inundation by April, 1914 (Ishihara et al., 1981). Although there are large uncertainties in the volume estimates, it appears that the volume erupted (in dense rock equivalent) agrees well with the volume needed to account for the observed ground deflation following the 1914 eruption. This correlation between deflation and volume erupted suggests there can also be correlation between magnitude of inflation and available magma to erupt, which is why deformation data and models, like the Mogi model, are widely used to monitor active volcanoes.