Learn R Programming

geophys (version 1.2-1)

mogiM: Mogi MOdel

Description

Mogi model deformation returns the deformation from a point source presurized inflation in an elastic medium.

Mogi's model (point source in elastic half-space). computes radial and vertical displacements Ur and Uz, ground tilt Dt, radial and tangential strain Er and Et on surface, at a radial distance R from the top of the source due to a hydrostatic pressure inside a sphere of radius A at depth F, in a homogeneous, semi-infinite elastic body and approximation for A << F (center of dilatation). Formula by Anderson [1936] and Mogi [1958].

Usage

mogiM(R = 1, F = 1, A = 0.1, P = 1e+05, E = 1e+10, nu = 0.25)

Arguments

R
Hoirizontal Distance frm source, m
F
Depth below surface, m, positive down
A
radius of magma chamber
P
hydrostatic pressure change in the sphere
E
elasticity (Young's modulus)
nu
Poisson's ratio

Value

  • list:
  • urradial displacements Ur
  • uzvertical displacements Uz, Uz > 0 = UP
  • dtground tilt Dt
  • erradial strain Er
  • ettangential strain Et on surface

Details

Original paper by Mogi used poisson's ratio equale to 0.25, i.e. lame parameters lambda and nu were equal.

References

Anderson, E.M., Dynamics of the formation of cone-sheets, ring-dikes, and cauldron-subsidences, Proc. R. Soc. Edinburgh, 56, 128-157, 1936.

Mogi, K., Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them, Bull. Earthquake Res. Inst. Univ. Tokyo, 36, 99-134, 1958.

Examples

Run this code
data(PXY)

delV = 2.3E13/(100^3)    ##### (convert to meter^3 from cm^3)
F = 2.8E5/100      ##### (convert to meter from cm    )

EX = seq(from=0, by=100, to= 9000)


Atest = mogiM(R=EX,F=F,A=delV)



 plot(PXY, pch=6, col='purple', xlim=c(0,9), ylim=c(0, 1) )
    ###  model
    lines(EX/1000, Atest$uz/max(Atest$uz))


############ best fit   optimization

library(stats)

    fr<-function(x)
      {

        Atest = mogiM(R=PXY$x*1000 ,F=x[1],A=x[2])


        rms = sum ( (PXY$y - Atest$uz/max(Atest$uz))^2 )

        return(rms)
      }
xin = c(2600, 2.0e+07)

FOUT = optim(xin , fr)

  Btest = mogiM(R=EX,F=FOUT$par[1] ,A=FOUT$par[2])

   plot(PXY, pch=6, col='purple', xlim=c(0,9), ylim=c(0, 1) )
 
 lines(EX/1000, Btest$uz/max(Btest$uz))

Run the code above in your browser using DataLab