Modeling Crustal Deformation near Active Faults and

Volcanic Centers—A Catalog of Deformation Models 

U.S. Department of the Interior

U.S. Geological Survey

Techniques and Methods 13–B1

Chapter 1 of

Section B, Modeling of Volcanic Processes

Book 13, Volcanic Monitoring

Modeling Crustal Deformation near Active 

Faults and Volcanic Centers—A Catalog of 

Deformation Models

By Maurizio Battaglia, Peter F. Cervelli, and Jessica R. Murray

Chapter 1 of

Section B, Modeling of Volcanic Processes

Book 13, Volcanic Monitoring

Techniques and Methods 13-B1

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director 

U.S. Geological Survey, Reston, Virginia: 2013

For more information on the USGS—the Federal source for science about the Earth, its natural and living 

resources, natural hazards, and the environment, visit http://www.usgs.gov or call 1–888–ASK–USGS.

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials 

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:

Battaglia, Maurizio, Cervelli, P.F., and Murray, J.R., 2013, Modeling crustal deformation near active faults and 

volcanic centers—A catalog of deformation models: U.S. Geological Survey Techniques and Methods, book 13, 

chap. B1, 96 p., http://pubs.usgs.gov/tm/13/b1.
iii

Contents

Abstract...........................................................................................................................................................1

Introduction.....................................................................................................................................................1

Geodetic Transformations ............................................................................................................................2

ITRF05 to ITRF00....................................................................................................................................4

XYZ to LLH ..............................................................................................................................................5

LL to UTM ...............................................................................................................................................6

Local Coordinates (ENU) .....................................................................................................................8

XYZ to ENU....................................................................................................................................8

ENU to XYZ....................................................................................................................................8

Spherical Source (Magma Chamber).........................................................................................................9

Volume Change ...................................................................................................................................11

Ground Tilt............................................................................................................................................11

Internal Deformation and Strain.......................................................................................................12

Leading-Order Solution.............................................................................................................12

First Free-Surface Correction..................................................................................................13

Higher Order Cavity Correction ...............................................................................................14

Sixth-Order Free-Surface Correction.....................................................................................15

Internal Deformation .................................................................................................................16

Strain .........................................................................................................................................16

Examples of Surface Deformation, Tilt, Internal Deformation, and Strain for a 

Spherical Source...................................................................................................................19

Prolate Spheroid Source (Magma Conduit) ............................................................................................20

Coordinates and Displacement (YANGDISP) ........................................................................21

Primitive (YANGINT) ..................................................................................................................21

General Parameters (YANGPAR).............................................................................................24

Verification..................................................................................................................................25

Volume Change ...................................................................................................................................28

Ground Tilt............................................................................................................................................29

Strain ..................................................................................................................................................30

Examples of Parameters Calculated for a Prolate Spheroid.......................................................30

Sill-Like Source ............................................................................................................................................31

Volume Change ...................................................................................................................................35

Ground Tilt............................................................................................................................................35

Strain ..................................................................................................................................................36

Examples of Parameters Calculated for a Sill-Like Structure.....................................................36

Surface Deformation and Ground Tilt for Rectangular Dikes and Faults ...........................................37

Displacement.......................................................................................................................................38

Parameters for Displacements and Displacement Gradients Functions fi .............................40

Elementary Function f

1 for a Strike-Slip Fault..............................................................................42

Elementary Displacement ........................................................................................................42

Displacement Gradient .............................................................................................................42
iv

Contents—Continued

Surface Deformation and Ground Tilt for Rectangular Dikes and Fault—Continued

Elementary Function f2 for a Dip-Slip Fault...................................................................................42

Elementary Displacement ........................................................................................................42

Displacement Gradient .............................................................................................................43

Elementary Function f3 for a Tensile Crack .................................................................................43

Elementary Displacement ........................................................................................................43

Displacement Gradient .............................................................................................................43

Verification...........................................................................................................................................44

Strike-Slip Fault..........................................................................................................................44

Dip-Slip Fault...............................................................................................................................45

Tensile Crack ..............................................................................................................................46

Examples of Surface Deformation and Ground Tilt for Three Types of Faults..........................48

Internal Deformation and Strain for Rectangular Dikes and Faults ...................................................49

Source Geometry.......................................................................................................................49

Displacement..............................................................................................................................49

Step 1—Elementary Displacement and Gradient .......................................................50

Step 2—Combination of Values......................................................................................51

Step 3—Computation of Deformation in Okada’s Coordinate System.....................51

Step 4—Rotate Deformation Back to Cartesian Coordinate System.......................51

Displacement-Gradient Matrix................................................................................................52

Combination of Values for X-Derivatives .....................................................................52

Combination of Values for Y-Derivatives.......................................................................52

Computation of X and Y Derivatives ..............................................................................53

Rotation of the Displacement-Gradient Matrix Back to Cartesian 

Coordinate System..............................................................................................53

Strain .........................................................................................................................................54

General Parameters ...........................................................................................................................54

Okada (1992, table 6)..................................................................................................................54

Okada (1992, equation 14).........................................................................................................54

Okada (1992, footnote, table 6) ................................................................................................55

Singularity Condition iii .............................................................................................................55

Singularity Condition iv .............................................................................................................55

I Parameters ...............................................................................................................................56

Singularity Condition ii ....................................................................................................56

J and K Parameters...................................................................................................................56

Parameters E–Q .........................................................................................................................57

Function fA (Infinite Medium)...........................................................................................................58

Function fB (Free Surface Deformation) .......................................................................................60

Function fC (Depth Multiplied Term).............................................................................................62

Verification and Examples.................................................................................................................65

Strike-Slip Fault..........................................................................................................................65

Dip-Slip Fault...............................................................................................................................67

Tensile Crack ..............................................................................................................................69
v

Contents—Continued

Discretized Faults.........................................................................................................................................71

Fault Segments and Subsegments ..................................................................................................73

Tiles ..................................................................................................................................................74

Verification of the Discretized Fault Model ....................................................................................75

Strike-Slip Fault..........................................................................................................................76

Dip-Slip Fault...............................................................................................................................77

Tensile Crack ..............................................................................................................................78

Smoothing Operator D.................................................................................................................................79

Building the Smoothing Operator.....................................................................................................79

Computing the Location of the Center Point of Each Tile....................................................80

Computing Distances Between the Center Points of Tiles .................................................81

Computing the Smoothing Operator D ...................................................................................82

Case 1—Fault does not Break the Surface..................................................................82

Case 2—Fault Breaks the Surface ................................................................................83

Case 3—Dike Opening at the Base of the Dislocation...............................................83

D Matrix.......................................................................................................................................84

Verification...........................................................................................................................................84

Example for a Right-Lateral Strike-Slip Fault .................................................................................86

Summary........................................................................................................................................................88

References Cited..........................................................................................................................................88

Appendix 1. Mathematical Methods for Computing Displacement ...................................................90

Figures

1. Diagram showing coordinate systems......................................................................................3

2. Diagram showing definition of boundary conditions and geometry for a spherical

source ..........................................................................................................................................10

3. Graphs showing free-surface deformation caused by a pressurized spherical

magma chamber ........................................................................................................................10

4. Diagram showing comparison between the dimensionless volume change of a

pressurized sphere calculated by use of equation 19 and a finite element method 

numerical model .......................................................................................................................11

5. Diagram showing comparison between calculations of ground tilt in the east and

north direction by a two-dimensional analytical model, and a finite element

method numerical model ..........................................................................................................12

6. Diagram showing definition of dimensionless boundary conditions and 

coordinates for a spherical source ........................................................................................13

7. Graph showing first free-surface correction ........................................................................14

8. Graph showing higher order cavity correction at the free surface ..................................15

9. Graph showing sixth-order surface correction at ζ = 0 .....................................................16

10. Graphs showing internal deformation caused by a spherical source ..............................17

11. Graphs showing internal strain caused by a spherical source..........................................18
vi

Figures

12. Diagrams showing definition of coordinates system and boundary-value

problem for a spheroidal source..............................................................................................20

13. Flow chart for coding Yang and others’ (1988) expressions for displacement, tilt,

and strain ....................................................................................................................................20

14. Graphs showing free-surface deformation caused by the inflation of a vertical

prolate spheroid ( φ = ° 0 ) for different values of the geometric aspect ratio A ...............25

15. Graphs showing free-surface deformation caused by the inflation of a tilted

( θ unequal to 0) and rotated (φ unequal to 0) prolate spheroid .....................................26

16. Graphs showing internal deformation caused by the inflation of a tilted and 

rotated prolate spheroid............................................................................................................27

17. Graph showing verification of the empirical correction for the volume change

of a pressurized prolate spheroid ............................................................................................28

18. Graphs showing ground tilt caused by the inflation of a tilted and rotated prolate

spheroid........................................................................................................................................29

19. Graph showing internal strain caused by the inflation of a tilted and rotated 

prolate spheroid..........................................................................................................................30

20. Graph showing definition of boundary conditions and geometry for a sill-like

source ..........................................................................................................................................33

21. Graph showing matlab function fialko01.m flow chart for calculating and coding

the parameters for a sill-like source .......................................................................................33

22. Graphs showing surface deformation due to the inflation of a sill-like source;

internal deformation...................................................................................................................34

23. Graph showing volume change for a sill-like source ..........................................................35

24. Graphs showing ground tilt caused by the inflation of a sill-like source..........................35

25. Graphs showing internal strain caused by the inflation of a sill-like source ...................36

26. Program flow chart for rectangular dislocations .................................................................37

27. Diagrams showing Geometry of the dislocation model in a Cartesian coordinate

system; Geometry of the source model .................................................................................37

28. Vector and profile plots of the displacement field of a strike-slip 

fault ..............................................................................................................................................44

29. Vector and profile plots of the displacement field of a dip-slip 

fault ..............................................................................................................................................45

30. Vector and profile plots of the displacement field of a 

tensile crack ................................................................................................................................46

31. Profile plots showing tilt caused by a tensile crack.............................................................47

32. Program flow chart for rectangular dislocations..................................................................49

33. Vector and profile plots showing displacement fieldand straincaused by a

strike-slip fault ............................................................................................................................65

34. Graphs showing vector and profile plots of the displacement field and strain 

caused by a dip-slip fault .........................................................................................................67

35. Graphs showing vector and profile plots of the displacement field and strain

caused by a tensile crack .........................................................................................................69

36. Graph showing representation of the fault geometry by segments, subsegments,

and tiles ........................................................................................................................................71

37. The Matlab script for coding a discretized fault includes four major steps.....................72

38. Diagrams showing segment, subsegment, and tile geometry............................................74
vii

Figures—Continued

39. Graphs showing fault geometry and vector plot and profile plots of the

surface displacement caused by a right-lateral strike-slip fault........................................76

40. Graphs showing fault geometry and vector plot and profile plots of the surface

displacement caused by a reverse dip-slip fault ..................................................................77

41. Graphs showing fault geometry and vector plot and profile plots of the surface

displacement caused by the opening of a tensile crack .....................................................78

42. Image showing grid and tile geometry for smoothing matrix D .........................................80

43. Grid used to compute the second-order finite-difference operator D; fault divied into 

P=2 rows, N=9 columns, and 18 tiles of different width; smoothing operator D for the 

fault in B. ......................................................................................................................................81

44. Graph showing fault geometry used to validate the algorithm for the smoothing

operator D....................................................................................................................................84

45. Color maps for the smoothing operator, slip distribution, and finite difference 

second derivative 2 ∇ sfor the fault geometry shown in figure 44.......................................85

46. Graph showing fault geometry employed in the example ...................................................86

Tables

1. List of earthquake and volcano deformation sources described in this report.................2

2. The 14 transformation parameters between ITRF05 and ITRF00 .........................................4

3. Comparison of IRTF00 and IRTF05 parameters from itrf052itrf00.m with results from

the National Geodetic Survey code HTDP for epochs 2000 and 2010 ................................4

4. Comparison of Matlab function xyz2llh.m with NGS code HTDP ........................................5

5. Comparison of coordinates for the International Global Survey site FLIN computed

in the UTM system in ArcGIS and ll2utm.m (datum WGS 84) ...............................................7

6. Conversion of coordinates of GPS sites near the Augustine volcano from IRTF 

to ENU.............................................................................................................................................8

7. Surface deformation, tilt, and internal deformation .............................................................19

8. Corrections for internal deformation ......................................................................................19

9. Internal deformation and strain ...............................................................................................19

10. General parameters ...................................................................................................................24

11. Surface deformation and ground tilt .......................................................................................30

12. Internal deformation and strain................................................................................................30

13. Surface deformation and tilt.....................................................................................................36

14. Internal deformation and strain ...............................................................................................36

15. Strike-slip fault ...........................................................................................................................48

16. Dip-slip fault.................................................................................................................................48

17. Tensile crack ..............................................................................................................................48

18. Strike-slip fault (internal deformation and strain; function fA) ..........................................58

19. Dip-slip fault (internal deformation and strain; function fA)...............................................58

20. Tensile crack (internal deformation and strain; function fA)..............................................59

21. Scaling used to compute the dimensional displacement (function fA)............................59

22. Strike-slip fault (internal deformation and strain; function fB) ..........................................60

23. Dip-slip fault (internal deformation and strain; function fB) ..............................................60

24. Tensile crack (internal deformation and strain; function fB) .............................................61

25. Scaling (internal deformation and strain; function fB)........................................................61
viii

Acronyms

EDM Electronic Distance Meter

ENU local east-north-up Cartesian coordinate system

FEM Finite Element Method

GPS Global Positioning System

GTSM Gladwin Tensor Strainmeters

InSAR Interferometric synthetic aperture radar

ITRF International Terrestrial Reference Frame

LLH Latitude, longitude, height

LL Latitude, longitude

PBO Plate Boundary Observatory

UTM Universal Transverse Mercator

WGS 84 World Geodetic System 1984

Tables—Continued

26. Strike-slip fault (internal deformation and strain; function fC) .........................................62

27. Dip-slip fault (internal deformation and strain; function fC) ..............................................63

28. Tensile crack (internal deformation and strain; function fC)..............................................64

29. Scaling (internal deformation and strain; function fC)........................................................64

30. Internal deformation and strain around a strike-slip fault...................................................66

31. Internal deformation and strain around a dip-slip fault ......................................................68

32. Internal deformation and strain around a tensile crack. ....................................................70

33. Parameters for strike-slip fault segments .............................................................................76

34. Parameters for dip-slip fault segments...................................................................................77

35. Parameters for tensile-crack fault segments........................................................................78

36. Values of fault parameters .......................................................................................................85

37. Values of fault parameters........................................................................................................86

38. Values of the smoothing operator D .......................................................................................87

39. Values of the slip vector U ........................................................................................................87

40. Values of the product DU .........................................................................................................87

41. Values of the finite-difference second derivative 2 ∇ s ........................................................87
Introduction 1

Abstract

This manual provides the physical and mathematical concepts for selected models used to 

interpret deformation measurements near active faults and volcanic centers. The emphasis is on 

analytical models of deformation that can be compared with data from the Global Positioning 

System (GPS) receivers, Interferometric synthetic aperture radar (InSAR), leveling surveys, 

tiltmeters and strainmeters. Source models include pressurized spherical, ellipsoidal, and 

horizontal penny-shaped geometries in an elastic, homogeneous, flat half-space. Vertical dikes 

and faults are described following the mathematical notation for rectangular dislocations in 

an elastic, homogeneous, flat half-space. All the analytical expressions were verified against 

numerical models developed by use of COMSOL Multyphics, a Finite Element Analysis 

software (http://www.comsol.com). In this way, typographical errors present were identified 

and corrected. Matlab scripts are also provided to facilitate the application of these models. 

Introduction

Deformation of the Earth’s crust can arise from tectonic and volcanic forces and from 

human activities such as aquifer withdrawal or geothermal exploitation. Mathematical 

models of how the crust deforms in response to different physical processes are required 

to distinguish among possible sources of deformation. These models can also help to 

characterize deformation-driving processes and to constrain source location, size, orientation, 

and strength. This information is valuable for hazards forecasting and mitigation, assessment 

of anthropogenic environmental impact, land-use planning, and other applications. A key 

assumption behind geodetic monitoring of faults and volcanoes is that ground deformation of 

the Earth’s surface reflects tectonic and volcanic processes at depth (for example, fault slip and 

mass transport) that are transmitted to the surface through the mechanical properties of the crust 

(Dzurisin, 2003). 

This paper discusses the application of analytical models to fit and interpret ground 

deformation. Analytical models offer a closed-form description of the source of crustal 

deformation (table 1). This means that, in principle, it is possible to readily infer the relative 

importance of any of the source parameters. Although analytical models are based on numerous 

simplifications (for example, the assumption that the crust is a homogenous, isotropic, elastic, 

flat half-space) that make the set of differential equations describing the problem tractable, 

they can take into account a vast array of source geometries. The careful use of analytical 

models, together with high-quality data sets, can yield valuable insights into the nature of the 

deformation source.

This report is intended as a practical reference for readers interested in applying 

mathematical models to investigate volcano and earthquake physics. Under each of the models 

described, the report first explains the assumptions, applications, and limitations of the model, 

and then introduces the corresponding analytical formulation. Finally, the report provides tables 

with numerical examples and Matlab scripts to aid in the use and coding of the formulas. 

Modeling Crustal Deformation near Active 

Faults and Volcanic Centers—A Catalog of 

Deformation Models

By Maurizio Battaglia1,2, Peter F. Cervelli1, and Jessica R. Murray1

1

 U.S. Geological Survey, 345 Middlefield Road, Menlo Park, California, 94025.

2

 Department of Earth Sciences, Sapienza University of Rome, Italy.
2 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

A zipped file with MATLAB scripts can be downloaded from pubs.usgs.gov/tm/13/

b1. The scripts are organized in folders, with each folder corresponding to a chapter in this 

publication (for example, the folder named “Geodetic transformations” contains MATLAB 

scripts for the geodetic transformation discuss in the section “Geodetic transformations” of 

this report). Folders may include subfolders with results from COULOMB 3.3 (“coulomb”) 

and Finite Element Method models (“FEM”), scripts verifying the algorithms (“verification”) 

and MATLAB functions (“functions”). The zipped file includes MATLAB functions for all 

the sources discussed in this publication and examples of the inversion of Global Positioning 

System (GPS) data for a spherical, spheroidal and sill-like source (the MATLAB Optimization 

Toolbox is required to run these examples). Online help can be obtained for each MATLAB 

function by typing “nameofthefunction help” in the MATLAB command line.

This publication is also designed to allow the user to turn directly to the desired model 

without reading any other sections (although a general understanding of the concepts and 

calculus behind the models from other literature will could enhance user experience, see for 

example Segall, 2010). As a result, some repetitions will be apparent if the manual is read 

sequentially. Many of the models described in this publication were taken from other credited 

or standard sources, but the authors themselves derived or rearranged many equations for this 

publication. All formulas have been checked for possible errors and verified against numerical 

models (COULOMB 3.3; http://earthquake.usgs.gov/research/modeling/coulomb/; FEM 

models). Since our emphasis is on describing the characteristics of the models and how they are 

used, the reader interested in the derivation and calculus of original equations is referred to the 

original literature. To limit the number of typographical errors, the formulas have been edited 

directly from the corresponding Matlab functions by using Mathtype. 

Geodetic Transformations

Coordinates can be expressed in systems defined by various datums, all of which are 

related to each other through geometrical transformations. Geodetic datums describe the size 

and shape of the Earth, and the origin, orientation, and time derivatives of the coordinate 

system (fig. 1). The list below includes four common coordinate systems:

1. global Cartesian (XYZ) system—International Terrestrial Reference Frame 2005 (ITRF05; 

http://itrf.ensg.ign.fr); units are in meters.

2. global geographic system, either latitude, longitude, height (LLH) or latitude, longitude 

(LL)—International Terrestrial Reference Frame 2000 (ITRF00–WGS 84 ellipsoid; http://

itrf.ensg.ign.fr); units are in degrees for latitude and longitude, and in meters for height. 

3. local east-north-up (ENU) Cartesian coordinate system; units are in meters.

4. local Cartesian (XY) system—Universal Transverse Mercator coordinate system (UTM–

WGS 84 ellipsoid); units are in meters.

Table 1. List of earthquake and volcano deformation sources described in this report.

Model Source Reference

Volcano source

Sphere Magma chamber McTigue (1987)

Spheroid Conduit Yang and others (1988)

Pennyshaped crack Sill Fialko and others (2001)

Tensile dislocation Dike Okada (1985), Okada (1992)

Earthquake source

Rectangular dislocation Dip- and strike-slip fault 

(single rectangular segment)

Okada (1985), Okada (1992)

Superposition of rectangular 

segments

Fault composed of several 

rectangular segments

Jonsson and others (2002)
Geodetic Transformations 3

λ

φ

Z

N E U

Y

X

Figure 1. Coordinate systems. XYZ, global Cartesian coordinate 

system; ENU, local Cartesian coordinate system; ( ) λ φ, global 

geographic coordinate system (latitude, longitude).

WGS 84 (World Geodetic System), released in 1984 and last revised in 2004 (National Imagery 

and Mapping Agency (NIMA), 2004) is the standard datum used by the Global Positioning 

System. 

GPS solutions from processing software like GIPSY (https://gipsy-oasis.jpl.nasa.gov) or 

GAMIT/GLOBK (https://www-gpsg.mit.edu/~simon/gtgk) are usually defined in the ITRF 

(XYZ) Cartesian coordinate system. GPS cooordinates can also be referenced to a local 

topocentric datum. The origin of this datum is an arbitrary point on the surface of the Earth 

(fig. 1). It has three right-handed orthogonal axes: E (for “east”) is in the local horizontal plane 

and points to the geographic east, N (for “north”) is in the local horizontal plane and points 

to the geographic north, U (for “up”) is vertical (in other words, perpendicular to the local 

equipotential surface) and points upwards; units are meters. UTM is a local Cartesian (XY) 

coordinate system based on the Universal Transverse Mercator projection and defined for the 

whole Earth (U.S. Army, 1987); units are meters. Modeling deformation velocities from GPS 

requires the transformation of the original ITRF coordinates into local Cartesian coordinates. 

Coordinates in ITRF05 can be transformed into ITRF00 by means of a 14-parameter 

Helmert transformation. ITRF00 and WGS 84 (G1150), the latest realization of the WGS 84 

coordinate system, are the same within one cm (NIMA, 2000 and 2004). The OpenGIS Web 

Map Service (www.opengeospatial.org/standards/wms) offers georeferenced maps in UTM 

(WGS84). The report provides formulas and Matlab functions to transform coordinates from 

1. ITRF05 to ITRF00,

2. ITRF (XYZ) to LLH,

3. LL to UTM (XY), and

4. ITRF (XYZ) to the east-north-up coordinate system (ENU).
4 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

ITRF05 to ITRF00

The 14-parameter Helmert transformation is (Dawson and Woods, 2010)

X

Y

Z

d dr t t

d dr t t

d dr t t

x x

y y

z z

00

00

00

0

0

0

















=

+ − ( )

+ − ( )

+ − ( )





















+ +  + − ( )  

















 1 0

05

05

05

sc src t t

X

Y

Z

R , (1)

where d d d x y z 

  is the origin shift in m, xyz   dr dr dr   is the shift rate in m/yr, t is the time (epoch) in years,

0t = 2000.0, sc is the dimensionless scaling in ppb and src the scaling rate in ppb/yr (table 2).

The rotation matrix Ris the sum of steady (Rs) and time dependent (Rt) components (see table 3 for a comparison 

between itr052itrf00.m and standard NGS software)

R R s t =

−

−

−





















=

( ) − − − ( )

−

1

1

1

0 R R 0 0

R R

R R

Rr t t Rr t t

z y

z x

y x

z y

and Rr t t Rr t t

Rr t t Rr t t

z x

y x

( ) − ( ) −

( ) − − − ( )





















0 0

0 0

0

0

. (2)

Table 3. Comparison of IRTF00 and IRTF05 parameters from itrf052itrf00.m with results from the 

National Geodetic Survey code HTDP for epochs 2000 and 2010. 

[The difference ∆ between the coordinates is less than 1 millimeter. See http://www.ngs.noaa.gov/TOOLS/Htdp/Htdp.

shtml; ITRF, International Terrestrial Reference Frame; HTDP, horizontal time-dependent positioning; coordinates for 

both epochs are for site FLIN.; Δ, difference]

Epoch 2000.01.01

X 

(meters)

Y 

(meters)

Z 

(meters)

ITRF05 -766174.473 -3611375.309 5184056.243

ITRF00 (HDTP) -766174.473 -3611375.311 5184056.239

ITRF00 -766174.473 -3611375.311 5184056.239

∆ -0.0002 -0.0002 0.0003

Epoch 2010.05.14

X 

(meters)

Y 

(meters)

Z 

(meters)

ITRF05 -766174.661 -3611375.333 5184056.209

ITRF00 (HDTP) -766174.664 -3611375.337 5184056.191

ITRF00 -766174.664 -3611375.337 5184056.191

∆ 0.0001 -0.0002 -0.0001

Table 2. The 14 transformation parameters between ITRF05 and ITRF00. 

[See http://itrf.ensg.ign.fr/ITRF_solutions/2005/tp_05-00.php; ITRF, International Terrestrial Reference Frame; dx, dy, and dz, Cartesian parameters measured in 

meters; Rx, Ry, and Rz, angular coordinates; mas, milliarcseconds; sc, scaling; ppb, part per billion; drx, dry, and drz, Cartesian shift rates; m/yr, meters per year; 

Rrx, Rry, and Rrz, angular shift rates; mas/yr, milliarcseconds per year; src, scaling rate; ppb/yr; part per billion per year].

dx (meters) dy (meters) dz (meters) Rx (mas) Ry (mas) Rz (mas) Sc (ppb)

0.0001 -0.0008 -0.0058 0.000 0.000 0.000 0.40

drx (m/yr) dry (m/yr) drz

 (m/yr) Rrx (mas/yr) Rry (mas/yr) Rrz (mas/yr) Src (ppb/yr)

-0.0002 0.0001 -0.0018 0.000 0.000 0.000 0.08
Geodetic Transformations 5

XYZ to LLH

The formulas implemented in xyz2llh.m compute longitude, latitude (both in decimal 

degrees), and height (in meters) from the ITRF00 cartestian coordinates XYZ. Geodetic 

positioning LLH is in the WGS 84 (G1150) reference frame. Horizontal datum is referred to 

WGS 84, and the vertical datum to the WGS 84 ellipsoid. Sign conventions are positive for E 

and longitude; positive for N and latitude. The constants for the WGS 84 ellipsoid are (NIMA, 

2000, paragraph 3.2)

 

a

f

=

=

6378137

1 298 257223563

semimajor axis, in meters, and

/ . flattening.

(3)

The geometrical parameters of the ellipsoidal coordinate system (Hoffmann-Wellenhof 

and others, 1997, equations 10.11–10.13) are

 

b f a

e f f

= − ( )

= −

1

2 2 2

semi-minor axis, in meters

eccentricity squared

second numerical eccentricity

radi

e a b b

p Y X

2

2 2 2

2 2

φ = −

= +

( )

us of a parallel

(4)

The longitude (φ), latitude (λ), and height (H) are given by

φ = − 



 



 − tan 1 X

Y , (5)

λ

θ

θ

θ φ = 













 = 



 





+

−

− − tan , tan

sin

cos

1 1 2

3

2 3

Z e b

p e a

Za

pb , (6)

2 2

, and

1 in cos s

p a H N N

e

= − λ =

− λ

, (7)

where 

 N is the radius of curvature. 

Finally, longitude and latitude are converted to degrees and the sign is changed so that E 

longitude is positive (table 4).

Table 4. Comparison of Matlab function xyz2llh.m with NGS code HTDP. 

[See http://www.ngs.noaa.gov/TOOLS/Htdp/Htdp.shtml. Values are for the International Global Navigation Satellite 

System Survey site FLIN. ITRF, International Terrestrial Reference Frame; WGS, World Geodetic System; HDTP, 

horizontal time-dependent positioning; DD, decimal degrees; X, Y, and Z, Cartesian coordinates; Long, longitude; Lat, 

latitude; Δ, difference]

X (meters) Y (meters) Z (meters)

ITRF05 -766174.473 -3611375.309 5184056.243

Long (DD) Lat (DD) Height (meters)

WGS 84 (HDTP) -101.978033 54.725584 311.5070

WGS 84 -101.978033 54.725584 311.5068

∆ -0.000000 -0.000000 -0.0002
6 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

LL to UTM

The Matlab function ll2utm.m computes the E (east, meters) and N (north, meters) 

coordinates on the UTM grid from the WGS 84 (G1150) longitude (φ) and latitude (λ). 

According to the Defense Mapping Agency (1989) report, “The computations (* * *) are 

accurate to the nearest 0.001 arc second for geographic coordinates and to the nearest 0.01 

meter for grid coordinates”, where 0.001 arc second = 7 3 10− × decimal degrees = 9 4.848 10− ×

rad. The UTM horizontal datum is the WGS 84 ellipsoid, with sign convention positive for E 

longitude and positive for N latitude.

The first step is to set the constants for the WGS 84 ellipsoid (NIMA, 2000, 

paragraph 3.2):

a

f

=

=

6378137

1 298 257223563

semimajor axis, in meters, and

/ . flattening.

(8)

and the WGS 84 ellipsoid parameters (Defense Mapping Agency, 1989, paragraph 2–2.1): 

b f a

e f f

= − ( )

= −

1

2 2 2

semiminor axis, in meters

eccentricity squared

e a b b second numerical eccentricity

v a

e

2

2 2 2

2 2 1

φ

φ

= −

=

−

( )

sin

radius of curvature in the prime vertical

(9)

( ) ( )

( ) ( )

( ) ( )

( )

( ) ( )

( ) ( ) ( ) ( )

23 45

2 34 5

23 45

34 5

4 5

2

 1 5 4 81 64

 3 / 2 7 8 55 64

 15 /16 0.75

(35 / 48) 11 /16

315 / 512 

sin 2 sin 4 sin 6 sin 8

 

f n f

Aa n nn nn

B an n n n n

C an n n n

D an n n

E an n

SA B C D E

φ

φ

φ

φ

φ

φφ φ φ φ

= −

= −+ − + −    

= −+ − +    

= −+ −    

= −+

= −

= φ− φ + φ − φ + φ

, (10)

where S is the meridional arc. 

The UTM projection parameters (Defense Mapping Agency, 1989, paragraph 2–2.2) are 

the longitude measured with respect to the central meridian (λ0 ), the central scale factor k0 , 

false easting (FE), and false northing (FN):

0 0 , 0.9996

0 or 10,000,000 500,000

k

FN FE

∆λ = λ − λ =

= =

, (11)

where FN is 0 for the northern hemisphere and 10,000,000 for the southern hemisphere.
Geodetic Transformations 7

The following terms are used to calculate the general equations which follow in this 

section (Defense Mapping Agency, 1989, paragraph 2–2.3):

T Sk

T v k

T T e e p

1 0

2 0

3 2

2 2

2

2

2

2

2

12 5 9 4

=

=

= ( ) − + +

sin cos

cos tan cos

φ φ

φ φ φ φ cos

4 ( ) φ

(12)

T T e e 4 2 p

4 2 4

2

2 2

2 = ( ) cos t φ φ 360 61 58 − + an tan c φ φ + − 270 330 os tan c φ os φ

2

2

2 4

2

3 6 2

2

2 4

2

4

445 324 680

88

φ

φ φ φ φ φ φ φ

φ

(

+ + − +

+

e e e

e

cos cos tan cos

cos tan cos tan cos

cos

8 2

2

3 6 2

2

4 8

5 2

6

600 192

20160

φ φ φ φ φ

φ

− − φ φ ) = ( )e e

T T 1385 3111 543 2 4 6

6 0

( ) − + −

=

tan tan tan

cos

φ φ φ

T v φk

T T e

T T

7 6

2 2

2

2

8 6

4 2

6 1

120 5 18

= ( )( ) − +

= ( ) − +

cos tan cos

cos tan ta

φ φ φ

φ φ

φ

n cos tan cos

cos cos

4

2

2 2

2

2

2

2 4

2

3 6

14 58

13 4 6

φ φ φ φ

φ φ

φ φ

φ φ

( + − +

+ + −

e e

e e 4 24

5040 61 479

2

2

2 4 2

2

3 6

9 6

6

tan cos tan cos

cos ta

φ φ φ φ

φ

φ φ e e

T T

− ) = ( ) − n tan tan 2 4 6 ( ) φ φ + − 179 φ

The general formulas for the conversion of geographic coordinates to north (N) and east 

(E) grid coordinates (Defense Mapping Agency, 1989, paragraph 2–5; table 5) are

N FN T T T T T

E FE T T T T

= + + + + +

= + + + +

1 2

2

3

4

4

6

5

8

6 7

3

8

5

9

7

∆ ∆ ∆ ∆

∆ ∆ ∆ ∆

λ λ λ λ

λ λ λ λ

(13)

Table 5. Comparison of coordinates for the International Global Survey site FLIN computed in 

the UTM system in ArcGIS and ll2utm.m (datum WGS 84). 

[WGS, World Geodetic System; UTM, Universal Transverse Mercator; GIS, geographic information system; DD, 

decimal degrees; X, Y, and Z, global Cartesian coordinates; Long, longitude; Lat, latitude; Δ, difference]

Long (DD) Lat (DD)

WGS 84 (G1150) -101.9780327 54.7255844

E (meters) N (meters)

UTM (ArcGIS) 308230.414 6068325.668

UTM 308230.413 6068325.666

∆ -0.001 -0.002
8 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Local Coordinates (ENU)

XYZ to ENU

The Matlab function xyz2enu.m transforms global Cartestian coordinates (XYZ) into local 

east-north-up (ENU) Cartesian coordinates. The origin of the local coordinate system is set to 

the minimum value of the (XYZ) coordinates. The transformation comprises two steps:

1. Transform the origin from XYZ coordinates ( ) 00 0 X YZ , , to LL coordinates ( ) 0 0 φ λ, , and 

2. Convert XYZ coordinates to ENU (Hoffmann-Wellenhof and others, 1997, p. 149) through 

the equation

E

N

U

















=

−

− −

sin cos

sin cos sin sin cos

cos co

λ λ

φ λ φ λ φ

φ

0 0

0 0 0 0 0

0

0

s c λ φ os sin s λ φ in 0 0 0 0

0

0

0

















−

−

−

















X X

Y Y

Z Z

. (14)

ENU to XYZ

This function (enu2xyz.m) implements the inversion of the XYZ to ENU transformation. 

A local reference point (for example, the origin) in XYZ is used to tie the coordinates to ITRF 

(see example in table 6 for the Augustine volcano in Alaska):

X

Y

Z

















=

− −

−

sin sin cos cos cos

cos sin sin cos

λ φ λ φ λ

λ φ λ

0 0 0 0 0

0 0 0 φ λ

φ φ

0 0

0 0

0

0

0 0

sin

cos sin

































+

















E

N

U

X

Y

Z

. (15)

Table 6. Conversion of coordinates of GPS sites near the Augustine volcano from IRTF to ENU.

[See Cervelli and others (2006) for original coordinates; ITRF, International Terrestrial Reference Frame; ENU, east north up (local 

coordinates); (λ0

, φ 0) are the latitude and longitude of the local origin (Augustine volcano); (X0,Y0,Z0) are the ITRF Cartesian 

coordinates of the local origin (Augustine volcano). Footnotes 1 and 2 are the numerical values of (λ0

, φ0) and (X0,Y0,Z0)]

ITRF coordinates

X 

(meters)

Y 

(meters)

Z 

(meters)

AC59 -2900773.239 -1440890.357 5476476.756

AUGL -2911728.045 -1461015.567 5465156.878

AV01 -2915632.507 -1456172.186 5464818.582

AV02 -2916881.537 -1458858.111 5463145.488

AV03 -2913037.130 -1456338.750 5466001.512

AV04 -2915070.452 -1456916.823 5465417.529

AV05 -2914535.420 -1458049.359 5465541.274

ENU coordinates

E 

(meters)

N 

(meters)

U 

(meters)

AC59 2606.223 9556.053 8918.654

AUGL 15009.630 -5668.214 -595.141

AV01 -7450.800 240.078 2191.218

AV02 -3437.328 6941.531 -1299.731

AV03 2404.740 1320.257 -752.401

AV04 679.336 -816.519 4293.685

AV05 -16943.480 11095.571 238.721

1

(λ0, φ0): (-153.394502 59.325519).

2

(X0,Y0,Z0): (-2916881.537 -1461015.567 5463145.488).
Spherical Source (Magma Chamber) 9

Spherical Source (Magma Chamber)

The deformation due to an expanding or contracting magma chamber has frequently been 

modeled by a dilatation source in an elastic half space. The most commonly used is the Mogi, 

or point, source (Masterlark, 2007). The model simulates a small spherical source embedded in 

a homogeneous, isotropic, elastic, flat half-space (fig. 2). 

The appeal of Mogi’s model lies in its combination of computational simplicity and 

remarkable ability to predict radially symmetric deformation caused by magma intrusion. The 

accuracy of an interpretation based on Mogi’s model, however, is subject to the suitability of 

the initial assumptions—an often overlooked consideration (Masterlark, 2007). For example, 

Mogi’s point-source model can explain stresses and displacements far away from the magma 

chamber, but the stresses are infinite at the source. In contrast, McTigue’s (1987) formulation 

provides an analytical solution with higher order terms that account for the finite shape of a 

spherical body. Thus, the local stresses at, and far from, the boundary of a chamber remain 

finite and can be calculated. McTigue’s results (1987) for vertical ( uz ) and radial ( ur) surface 

deformation caused by a pressurized ( ∆P ) spherical magma chamber of radius a and at 

position ( , x y , ) z 0 0 0 (fig. 2) can be written in the form

u

Pa z

r z

a

z z = − ( )

( ) +

−





 



 +

( ) − − ( ) −

− 1 1 1

2 7 5

15 2

4 7

3

0

2

0

2 3 2

0

3

ν

μ

ν

ν

∆ ν

5

0

2

2

0

2 ( ) ν +





































z

r z

, (16)

and

u

Pa r

r z

a

z r = − ( )

( ) +

−





 



 +

( ) − − ( ) −

− 1 1 1

2 7 5

15 2

4 7 5

3

2

0

2 3 2

0

3

ν

μ

ν

ν

∆ ν

( ) ν +





































z

r z

0

2

2

0

2 , (17)

where z0 is the depth of the source (a positive number), ν the Poisson’s ratio, and μ the 

shear modulus. A direct consequence of the assumption of a point source of dilatation is that 

the magma-chamber radius a and pressure change ∆P are inseparable because 3 ∆Pa is the 

strength of the point singularity. This is why point source models typically calculate volume 

rather than pressure changes.

When the radius a is small relative to the depth z0 , so that ( )3

0 a z is much less than 

1, the Mogi assumptions are reasonable, and the McTigue (1987) correction is not needed. 

For example, when a z = 0 2 , the correction is only 12 percent, and it decreases to 3 percent 

for a z = 0 3 . This means that, if the overall accuracy of the geodetic survey is within 

3–12 percent, the Mogi’s model is sufficient to fit the deformation. If the site effects induced by 

shallow heterogeneities and topography account for as little as 3–12 percent of the deformation 

signal, the Mogi’s model is again a good approximation. 

We can use equations 16 and 17 and to compute the three components of surface 

deformation (for example, the components measured by GPS receivers) aligned to geographic 

east (x-axis), north (y-axis) and up (z-axis), as shown in figure 3 (see table 7 at the end of 

chapter for numerical examples). 

u u u u

x

r u yr u x y z r r z 

  = 



 



. (18)
10 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

−5000 0 5000

−0.04

−0.03

−0.02

−0.01

0.00

0.01

0.02

0.03

0.04

X−displacement, in meters

−5000 0 5000

x, in meters

−5000 0 5000

EXPLANATION

FEM

3D

2D

Figure 3. Free-surface deformation caused by a pressurized spherical magma chamber. Comparison between 

the three-dimensional (3D) semianalytical solution (equation 22), the two-dimensional (2D) analytical model 

(equations 16 and 17), and a numerical finite element method (FEM) model of the same source. The source is 

at the origin of the coordinate system, and its parameters are a = 500 meters, z a 0 = = 2 1,000 meters, and 

∆P μ = 0 001273. .

Figure 2. Definition of boundary conditions and geometry for a spherical source (modified from 

McTigue, 1987). A, Vertical cross section showing he boundary conditions: the normal stresses 

σ σ rz zz and are equal to zero at the free surface, whereas the normal stress σn is equal to the 

pressure change at the surface of the spherical magma chamber, and the tangential stress σt is 

zero. r and z are local cylindrical coordinates, a is the radius of the spherical source, (x0

,y0,z0) is the 

source location, uz

 is the uplift. B, The local coordinate system at the earth’s surface. x (east) and y 

(north) are local Cartesian coordinates, ur

 is the radial displacement.

men12-7629_fig02

z0

a

z

σ r zz = σrz = 0

σn = -∆P

σt

 = 0

uz

ur

y (N)

x (E)

(x,y)

r2 = (x-x0)2 + (y-y0)2

A B

§

 

r

(x0,y0,z0)

§

(x0,y0,z0)
Spherical Source (Magma Chamber) 11

Volume Change

The change in the volume of the sphere is

∆ ∆ V a P a

z = +





 



















 π

μ

3

0

4

1 . (19)

When a z0

4 ( ) is much smaller than 1, equation 19 will return the well known result for the volume change caused by a 

pressurized point source (Amoruso and Crescentini, 2009; fig. 4).

Figure 4. Comparison between the dimensionless volume change ( ∆ ∆ V a ( ) π μ P 3 ) of 

a pressurized sphere calculated by use of equation 19 (analytical, solid line) and a finite 

element method numerical model (FEM, open circles). a/z0

 is the ratio of the radius of the 

magma chamber to the depth of the centroid of the chamber below the free surface. Source 

parameters for the FEM model are a = 1 000, meters and ∆P μ = 0 001273. .

a/z0

∆V/(πa3∆P/μ)

0.95

1.00

1.05

1.10

1.15

1.20

0.0 0.2 0.4 0.6 0.8

EXPLANATION

FEM

Analytical

Ground Tilt

Measuring tiny changes in the slope angle or tilt of the ground near a volcano is one of the oldest methods for monitoring 

deformation caused by moving magma. When magma forces the ground upward, the slope of adjacent areas will usually tilt 

away from the center of uplift by a small amount. Conversely, if the ground subsides as a consequence of depressurization of 

a magma body, the slope of adjacent areas will tilt toward the center of subsidence. The components of the tilt vector are the 

horizontal derivatives of the vertical deformation uz (fig. 5; see table 7 at the end of chapter for numerical example).

du

dx

Pa xz

r z

a

z

z = −( ) −

( ) +

− −





 



1 3  ( ) + − 1

7 5

1 5 1 18

3

0

2

0

2 5 2

0

3

ν

μ ν

ν

∆ . .75 2

1

0

2

2

0

2

3

0

2

( ) − +





































= −( ) −

ν

ν

μ

z

r z

du

dy

Pa yz

r

z ∆

( ) +

− −





 



 ( ) + − − ( )

+







 z

a

z

z

r z 0

2 5 2

0

3

0

2

2

0

2 3 1

7 5

1 5 1 18 75 2

ν . . ν ν





























. (20)
12 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Internal Deformation and Strain

The internal deformation and strain is better described using the dimensionless notation 

shown in fig. 6.The spatial coordinates are scaled by the depth of the spherical source z0 . The 

radius of the spherical cavity is defined by ε = a z0 . The deformation is scaled by z0 , the 

pressure change ∆P, and the shear modulus μ

ξ ψ ρ ξ ψ

μ μ

ζ ρ

= − = − = +

= =

x x

z

y y

z

u

z P

u

u

z P

u z

D

D r

D

D

0

0

0

0

2 2

3

0

3

3

0

3

( )

( )

( )

(

∆ ∆

) ( ) = ρ ζ u

3D

, (21)

The internal deformation caused by a spherical source is given by the superposition of a 

leading-order solution u ( ) 0 , a first free-surface correction u ( )1 , a higher order cavity correction 

u ( ) 3 , and a sixth-order free-surface correction u ( ) 6 (see tables 7 and 9 at the end of chapter for 

numerical examples):

u u u u u

u u u

D

D

ρ ρ ρ ρ ρ

ζ ζ ζ

ε ε ε

ε ε

( ) ( ) ( ) ( ) ( )

( ) ( ) ( )

3 0 3 1 3 3 6 6

3 0 3 1

= + + +

= + + 3 3 6 6

u u ζ ζ ε ( ) ( ) +

. (22)

At the free surface ( ζ = 0 ), equations return equations 16 and 17.

Leading-Order Solution

The internal displacements in cylindrical coordinates system are 

u

u

ρ

ζ

ε ρ

ρ ζ

ε ζ

ρ ζ

( )

( )

0 3

2 2 3 2

0 3

2 2 3 2

1

4 1

1

4

1

1

=

+ − ( ) 

 

= −

+ − ( ) 

 

. (23)

At the free surface ( ζ = 0 ), equations return McTigue’s (1987) equations 14 and 15.

Figure 5. Comparison between calculations of ground tilt in the east and north 

direction by a two-dimensional (2D) analytical model (equation 20), and a finite 

element method (FEM) numerical model. As in figure 3, source parameters are 

a = 500 meters, z a 0 = = 2 1,000 meters, and ∆P μ = 0 001273. .

men12-7629_fig05

−5000 0 5000

−1.0

−0.5

0.0

0.5

1.0

−1.0

−0.5

0.0

0.5

1.0

Tilt (east), in radians x 10-4

−5000 0 5000

X, in meters Tilt (north), in radians x 10-4

FEM

2D

EXPLANATION
Spherical Source (Magma Chamber) 13

First Free-Surface Correction

The displacements for the first-order free-surface correction are given by the appropriate 

combination of the following Hankel transforms (McTigue, 1987, p. 12939) 

u t e J t dt

u

t

ρ

ζ

ρ

ρ ζ σ ν ζ ρ

ρ ζ τ ν

( , )

( , )

= − ( ) −  

 ( )

= − ( ) −

−

∞

∫

1

2

1 2

1

2

2 1

1

0

1 t e J t dt

u t e J tR dt

t

t

ζ ρ

ρ ζ σ ν ζ

ζ

ζ

ζ



  ( )

= −  ( ) −  

 ( )

−

∞

∫ 1

0

0

0

1

2 ( , ) 2 1

∞

∞

∫

∫ = − ( ) −  

 u t e J ( ) tR dt t

ζ

ζ ( , ρ ζ) τ ν ζ 1

2

1 2 1 0

0

, (24)

where σ τ = = −

1

1

2

te t

 and R = + ρ ζ 2 2 . The integrals in equation 24 can be solved

numerically. At the free surface ( ζ = 0 ), equations return McTigue’s (1987) equations 24 

and 25 (fig. 7; see table 8 at the end of chapter for numerical examples). 

men12-7629_fig06

ζ0

a

ζ

σζζ = σρζ = 0 ρ

σn = -1

στ = 0

uζ

uρ

ψ (N)

ξ (E)

(ξ,ψ)

A B

§

 

ρ

(ξ0,ψ0,ζ0)

§

(ξ0,ψ0,ζ0)

Figure 6. Definition of dimensionless boundary conditions and coordinates for a spherical 

source. A, Vertical cross section showing the boundary conditions: the normal stresses 

σ σ ρζ ζζ and are equal to zero at the free surface, whereas the normal stress σn is equal -1 

at the surface of the spherical magma chamber, and the tangential stress στ is zero. ρ and ζ

are local cylindrical coordinates, a is the radius of the spherical source, (x0,y0,z0) is the source 

location, uζ is the uplift. B, The local coordinate system at the earth’s surface. ξ (east) and ψ

(north) are local Cartesian coordinates, uρ is the radial displacement.
14 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Higher Order Cavity Correction

The displacements for the higher order cavity correction are given by

u u u

u u u

R

R

ρ θ

ζ θ

θ θ

θ θ

( )

( )

sin cos

cos sin

3

3

= +

= −

, (25)

where 

sin ; θ cos ρ θ ζ

ρ ζ

= = −

= + ( ) −

R R

R

1

1 2 2

. (26)

McTigue’s (1987) equations 38–41 provide these expressions for uR and uθ :

u

D

R

P C

D

R

P

R

u C

R = − + − ( ) − 















= −

1

2

5 4 3

2

0

3

2 0 2

3 2

3

2

2

2

2

3

( )

( )

( )

( )

ν

θ 1 2 1

2

2 1 3

2 2

2 ( ) − +

















ν

θ

D

R R

dP

d

( )

, (27)

where the C and D coefficients are

C D 2 D 3 3

0

3 3

2

3 5 5 2

25 7 5

1

12

2

4 7 5

( ) ( ) ( ) = , −

( ) − = − + = −( ) − ε

ν

ν

ε

ν

ε

ν

ν

 and , (28)

Figure 7. First free-surface correction. Blue points, numerical solution 

of equation 24 at ζ = 0 ; red/green line, analytical representation of 

free-surface deformation, equations 16 and 17. uρ , dimensionless radial 

deformation; uζ , dimensionless vertical deformation. 

men12-7629_fig07

0 1 2 3 4 5

0

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50 First Free Surface Correction

Dimensionless distance ρ

Dimensionless deformation

Type of solution

uζ

uρ

Analytical

Numerical

EXPLANATION
Spherical Source (Magma Chamber) 15

and the Legendre’s polynomial and derivatives are

P P dP

d 0 2

2 2 1 1

2 = = , c ( ) 3 1 os θ − = 3cosθ

 and θ . (29)

Equation 25 return McTigue’s (1987) equations 46 and 47 (fig. 8) at the the free surface

(ζ = 0 ).

0 1 2 3 4 5 −0.005

0.000

0.005

0.010

0.015

0.020

0.025

0.030

0.035

0.040 Higher Order Cavity Correction

Dimensionless distance ρ

Dimensionless deformation

Type of solution uζ

uρ

Analytical

Numerical

EXPLANATION

Figure 8. Higher order cavity correction at the free surface, ζ = 0. Blue 

points, numerical solution of equation 25 at ζ = 0 ; red/green line, analytical 

representation of surface deformation, equations 16 and 17. uρ , dimensionless 

radial deformation; uζ , dimensionless vertical deformation. 

Sixth-Order Free-Surface Correction

The procedure to determine the sixth-order free-surface correction is identical to that for 

obtaining the first-order surface effect (equation 24). The following expressions describe the 

normal and shear stresses at the free surface:

u t t

e

u t

e

t

t

ρ

ζ

σ

ν

σ

ν

→ + ( ) −

→ −

−

−





3

2 7 5

7 5

2

2

, (30)

and

τ

ν

1

2

7 5  t

e−t

− . (31)

For ζ = 0 and the stresses given by equations 30 and 31, equation 24 returns McTigue’s 

(1987) equations 50 and 51 (fig. 9; see table 8 at the end of chapter for numerical examples).
16 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Internal Deformation

Figure 10 on the next page shows the verification for the expression for internal 

deformation (equation 22).

Strain

Strainmeters can resolve changes in strain of less than one part per billion (1 mm in 

1,000 km) over short periods, which makes them ideal for capturing transient deformation 

over time intervals ranging from seconds to months. Gladwin Tensor Strainmeters (GTSMs), 

the instruments used by the Plate Boundary Observatory (PBO), are designed to measure three 

components that describe the horizontal strain tensor: the areal strain εa, and shear strains γ1

and γ2 . Traditional engineering analysis has led to these quantitative definitions of quantities 

εa , γ1 , and γ2 : 

ε ε ε γ ε ε γ ε

ε ε

a xx yy xx yy xy

xx x yy

u y

x

u D D

=+ =− =

= ∂

∂ = ∂

, , ;

,

( ) (

1 2 2

3 3

and 

) ( ) ( )

, . ∂ = = ∂

∂

+

∂

∂















y 

u

y

u

x xy yx x y

D D

 and ε ε

1

2

3 3

. (32)

Because of the complexity of the equation 22 for the internal deformation, we compute the 

strain components shown in equation 32 by use of numerical derivatives (see appendix 1).

Figure 9. Sixth-order surface correction at ζ = 0. Blue points, numerical 

solution (equations 24, 30, and 31) at ζ = 0 ; red/green line, analytical 

representation of surface deformation, equations 16 and 17. uρ , 

dimensionless radial deformation; uζ , dimensionless vertical deformation.

men12-7629_fig09

0 1 2 3 4 5

−0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7 Sixth-order surface correction

Dimensionless distance ρ

Dimensionless deformation

Analytical

Numerical

Type of solution uζ

uρ

EXPLANATION
Spherical Source (Magma Chamber) 17

men12-7629_fig10

−0.03

−0.02

−0.01

0.00 0.00

0.01

0.02

0.03

X, in meters

FEM

3D

−0.03

−0.02

−0.01

0.01

0.02

0.03

0

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

−5000 0 5000

−0.05

0.00 0.00

0.05

X−displacement, in meters

−5000 0 5000

−0.05

0.05

Y−displacement, in meters 

−5000 0 5000

0

0.05

0.10

0.15

0.02

0.25

0.03

0.35

Z−displacement, in meters

EXPLANATION

Figure 10. Internal deformation caused by a spherical source. The first row shows internal deformation at 0z a = 0.5 , and the 

second row internal deformation at 0z a = . Blue points (3D), analytical solution, equation 22; dashed green line (FEM), finite 

element method numerical solution. Source parameters for the FEM model are 0z = 1,000 meters and ∆ μ= P 0.001273 . It 

was necessary to change some of the signs of ζ in the original integrals A7– A18 in McTigue (1987) to improve the fit between 

the analytical and the numerical solution at depth ζ greater than zero. 
18 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Figure 11. Internal strain caused by a spherical source. The first row shows internal strain at 0z a = 0.5 , and the second 

row internal strain at z a 0 = . Blue points (3D), analytical solution, equation 32; dashed green line (FEM), finite element 

method numerical solution. Source parameters for the FEM model are 0z = 1,000 meters and ∆ μ= P 0.001273 . 

men12-7629_fig11

−0.5

0.0

0.5

1.0

1.5

2.0

2.5

3.0

FEM

3D

−7

−6

−5

−4

−3

−2

−1

0

1

−7

−6

−5

−4

−3

−2

−1

0

1

X, in meters

Areal strain εa, × 10-4

Shear strain γ1 , × 10-4

Shear strain γ2, × 10-4

−5000 0 5000

0

2

4

6

8

−5000 0 5000

−2.5

−2.0

−1.5

−1.0

−0.5

0.0

0.5

−5000 0 5000

−2.0

−1.5

−1.0

−0.5

0.0

0.5

1.0

EXPLANATION
Spherical Source (Magma Chamber) 19

Table 7. Surface deformation, tilt, and internal deformation.

[Surface deformation, tilt and internal deformation at z = 0 from equations 18, 20, and 22. See figures 3 and 5. Parameters are a = 500 meters, 0 z a = = 2 1, 000

meters, and ∆ μ= P 0.001273. x and y, local Cartesian coordinates; , x y u u , and z u displacements at the free surface z ∂ ∂ u x

and z ∂ ∂ u y , tilt components; u u x

D

y

() () D , 3 3 , and, uz

( ) 3D , internal deformation]

x

(103

 meters)

y

(103 meters)

ux

(10-3 meters)

uy

(10-3 meters)

uz

(10-3 meters)

duz

/dx

(× 10-6)

duz

/dy

(× 10-6)

ux

(3D)

(10-3 meters)

uy

(3D)

(10-3 meters)

uz

(3D)

(10-3 meters)

-4.0 -4.0 -2.4949 -2.4949 0.6237 0.2275 0.2275 -2.4973 -2.4973 0.6305

0.0 0.0 0.0000 0.0000 134.7482 0.0000 0.0000 0.0000 0.0000 138.3455

2.0 2.0 8.8603 8.8603 4.4301 2.9846 2.9846 8.8631 8.8631 4.4876

Table 8. Corrections for internal deformation. 

[Corrections at ζ=0 from equations 24, 25, 30, and 31. See figures 7 through 9. Parameters are a = 500 meters, 0 z a = = 2 1, 000 meters and

∆P μ = 0.001273. ρ, dimensionless distance, equation 21; ur

( ) 1 and uz( ) 1 , first-order correction; ur

( ) 3

 and uz

( ) 3 , higher order correction;

ur

( ) 6

 and , uz

( ) 6

 sixth-order correction]

ρ ur

(1) (X 10-3) uz(1) (X 10-3) ur(3) (X 10-3) uz(3) (X 10-3) ur(6) (X 10-3) uz(6) (X 10-3)

0.0 0.0000 500.0000 0.0000 35.1751 0.0000 673.9130

2.0 89.4426 44.7213 -0.3950 0.9202 18.6663 8.9443

4.0 28.5339 7.1329 -0.4097 0.0697 1.7514 -0.0219

 Table 9. Internal deformation and strain. 

[Internal deformation and strain (z = a=500 meters) from equations 22 and 32. See figures 10 and 11. Parameters are a = 500 meters, 0

z a = = 2 1, 000

meters and ∆ = P μ 0.001273. x and y, local Cartesian coordinates; u u x

D

y

( ) D , 3 3( )and uz

( ) 3D , internal deformation; ε γ a 1 , and γ 2 , strain. The strain 

component γ1 ≡ 0 for x y = ] 

x

(X 103

 meters)

y

(X 103 meters)

ux

(3D)

(10-3 meters)

uy

(3D)

(10-3 meters)

uz

(3D)

(10-3 meters)

εa

(X 10-6)

γ2

(X 10-6)

x

(X 103 meters)

y

(X103 meters)

γ1

(X 10-6)

 -4.0 -4.0 -2.3320 -2.3320 0.5283 -0.4610 -1.6271 -4.0 0.0 -3.9171

 0.0 0.0 0.0000 0.0000 251.5401 639.9259 -0.0000 0.0 0.0 0.0000

 2.0 2.0 7.1229 7.1229 3.9175 -1.3665 -8.4895 2.0 0.0 -16.7005

Examples of Surface Deformation, Tilt, Internal Deformation, and Strain for a Spherical Source
20 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Prolate Spheroid Source (Magma Conduit)

A simple model of a volcanic system includes two principal elements: a magma reservoir and a conduit through which 

magma may reach the surface. When the volcano is quiescent, the conduit will close, allowing pressure to build up in the 

reservoir. Yang and others (1988)1 formulated an approximate but accurate half-space double-force and center-of-dilatation 

solutions for a dipping prolate ellipsoid in an elastic half-space. The general solution for a prolate spheroid depends on seven 

parameters (fig. 12): the dimensionless pressure change ∆P μ, the geometric aspect ratio A ba = between the semimajor axis a

and the semiminor axis b, the source location ( , x y , ) z 0 0 0 , the dip angle θ (measured from the free surface), and the strike (or 

azimuth) angle φ (measured clockwise from the positive North direction).

The expressions for the three components of deformation in the near field are quite complex. Coding the expressions for 

near-field displacement, tilt, and strain required writing the four nested functions described below (fig. 13). The formulas are 

introduced in the same order as their use in the corresponding Matlab functions.

x (East)

y (North)

θ

φ

z

b

a

(x0, y0)

σn = στ = 0

σn = -∆P

στ = 0

z0

yp

x θ p

z

b

a

z0

(0,0)

(side view)

A B

xp

yp

(0,0)

a

b

(top view)

C

Figure 12. Diagrams showing A, definition of coordinates system and boundary-value problem for a spheroidal source. x, y and z, 

local Cartesian coordinates; x0

, y0 and z0, source location; θ, dip angle; φ, strike (or azimuth) angle; a, semimajor axis; b, semiminor 

axis; the normal stress σn and shear stress στ are equal to zero at the free surface, whereas the normal stress σn is equal to the 

pressure change at the surface of the spherical magma chamber, and the shear stress στ is zero. B and C, coordinate system used 

by Yang and others (1988). xp

, yp and z, local Cartesian coordinate system; the source centroid is at the origin of the coordinate 

system, and the strike φ is equal to zero.

YANG

(main function)

calls YANGDISP to

compute 

3D deformation,

strain, and

tilt

YANGDISP

calls YANGPAR and 

YANGINT to compute the 

displacement caused by ó

a pressurized ellipsoid

YANGPAR

computes the input parameters of 

the spheroidal model

YANGINT

computes the óelementary deformation (primitive) 

ò

Figure 13. Flow chart for coding Yang and others’ (1988) expressions for displacement, tilt, 

and strain.

1

Corrections to Yang and others’ (1988) equations are reported in Newman and others (2006), and in this publication.
Prolate Spheroid Source (Magma Conduit) 21

Coordinates and Displacement (YANGDISP)

To compute the displacement, coordinates must be transformed from the Cartesian 

coordinate system (UTM) onto the coordinate system used by Yang and others (1988). The 

displacement vectors then can be back-transformed onto the original system (fig. 12). 

We first translate and rotate the coordinates 

x x x y y

y x x y y

p

p

= − ( ) − − ( )

= − ( ) + − ( )

cos sin

sin cos

φ φ

φ φ

0 0

0 0

, (33)

and then we compute the general expressions 

 

Up U c U c

Up U c U c

Up U c U

x

y

z

= − = − = −

= − = − = −

= = +

1 1

2 2

3 3

( ) ( )

( ) ( )

( )

ξ ξ

ξ ξ

 ξ ( ) ξ = −c

(34)

for the displacement in Yang and others’ (1988) coordinate system, where

 * †

†

0 11

1, 2,3

( , ,, ,, , ,,,,,, )

iii

i ipp

UUU i

U U x y zz a b ab P

=+ =

= θ xμν

(35)

are the primitive of the displacement for a prolate ellipsoid. The formulas for Ui

*

 and † Ui are 

given by equations 49 and 50. Finally, the displacement vector Up is back-transformed onto the 

original coordinate system (see figs. 14 and 15 in the Verification section, p. 25):

u Up Up

u Up Up

u Up

x x y

y x y

z z

= ⋅ + ⋅

= − ⋅ + ⋅

=

cos sin

sin cos

φ φ

φ φ . (36)

Primitive (YANGINT)

To compute the primitive * Ui and † Ui , we must first introduce new coordinates and 

parameters (see Yang and others, 1988, p. 4251):

 

 2 3

1 2 3 03 0

112 223 3

222 222

1 1 23 2

33 33

2 2 3 22 3

32 3 3 2 3

3 3

1 23

3 3

cos sin

 

 

 sin cos sin cos

cos sin cos sin

 

p p x x x y x zz x zz

yxy x y x y x

rx x q

R y

x x

rx x q x x

r r

yyR y

q

y y

q

x =x θ x =x θ

= = = − =+

= = −x = −x = +x

= θ− θ = θ+ θ

= θ+ θ =− θ+ θ

= −x

= ++ = ++

= +x

(37)

To make the formulas valid for dip angles θ unequal to 90 (in other words, a spheroid 

that is not vertical), we must correct the parameter C0 (see Yang and others, 1988, p. 4251, for 

the original expression):

0 0 C z = θ sin . (38)
22 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

To prevent a “divide by zero” warning at the origin, we must add 15

0 1 10− ε=⋅ at the 

denominator in the formula for β

β

θ θ

θ ε = + + ( )( ) +

+

q R q

y

2 2 3

1 0

cos s 1 in

cos

. (39)

Additional parameters used to simplify the expressions for the deformation are

3 13 2 3 2 3 3 3

3 1 3 3 3 2 3 2 3 ln( ) ln( ) ln( )

dr R r dq R q dy R y

lr R r lq R q ly R y

=+ = + = +

=+ = + = +

, (40)

 

A a

R dr b lr r

dr

A a

R dq

b lq

q

dq

1

1

1 3

1 3

3

3

1

1

2 3

1 3

3

3

*

*

= + +

 +



 





= − − +

 −

ξ

ξ















, (41)

1 1 3 3

1 2

2 13 2 2 3 3 3

3 3

3 13 2

1 2

A lr A lq R R

A R r lr A R q lq

r q A RA R

R R

x x =+ = −

=− = −

x x =+ = −

, and (42)

B a

R

b A a

R

b A

B

C

R

* 

 

= +





 



 + − ( ) +





 





= ( ) + −

1

1

1 2

1

2

1 2

0

2

2 3 4 2 ν

ξ ξ A C lq 2 0 3 −

. (43)

The following four parameters F are called in * Ui and † Ui (see equations 49 and 50). They 

were modified from the original formulation of Yang and others (1988) to improve the fit with 

the numerical FEM solution for internal deformation (fig. 16):

F

F

F

F

1

2

1

2

0

0

0

0

*

*

=

=

=

=

. (44)
Prolate Spheroid Source (Magma Conduit) 23

The parameters f also are called in * Ui and † Ui (see equations 49 and 50). While the 

formulas for f

1 and 2f are from Yang and others (1988), the equation for f3 was corrected 

by Newman and others (2006):

( )

1

1 1

1 2 1 12 13 3 3 3

3

3 tan sin 2 .* tan 2 4 y f y ly y lq q y lq x dy cos cos

− x − β = + θ − + β+ − θ θ

. (45)

f y

dy q lq q ly y R y 2

2

3

2 2 3 2 3 1

1

2 3

3

= + ( − + 2 + − ( ) ξ −

θ

θ θ β θ

cos

sin sin tan cos )

 − + 2 ( ) − 2

2 3 3 3 3 cos

cos

θ

θ

A x ly q lq

, and (46)

f q lq q ly y A q ly 3 2 3 2 3 1

1

2 3 3

1 = 2 2

cos

sin tan sin

θ ( ) − + θ β + + θ ξ − − . (47)

The following coefficients are used to compute the displacement

( )

2 3 * † *† and 2

16 1

ab

c c cP x = = μ −ν

. (48)

The displacement components are given by equations 2 to 7 of Yang and others (1988). 

They are the primitive of equation 1 of Yang and others (1988). Compared to the original 

formulation, we deleted the z-contribution in U2

* and U3* to improve the fit to internal 

deformation computed by use of the FEM model (fig. 16)

U c A y A y F y

U c A r A

1 1 1 1 1 1 1

2 1 2 1

3 4

3 4

* * * * *

* * * sin

= +  ( ) − +  



= + ( ) −

ν

θ ν * * * *

* * * *

cos

-cos

q F q B F

U c A r A q

2 1 2 2

3 1 2 1 2 3 4

( ) + + − ( ) 

 

= + ( ) −

θ

θ ν ( ) − + + ( ) 

  F q 1 2 B F2

* * * sin θ

, (49)

( ( ) ) ( )( )

( ( ) ) ( )( )

( ) ( ) ( ( ) )

( ( ) ) ( )( )

† †

1 1 1 11 11 1

† †

2 1 2 12 12 2

22 3 32

† †

3 1 2 12 12

 3 4 41 1 2

 sin 3 4 4 1 1 2

 +4 1 cos cos 3 4

 cos - 3 4 4 1 1 2

U c Ay Ay Fy f

U c Ar Aq Fq f

AA A AF

U c Ar Aq Fq

= + − ν + − −ν − ν    

= θ + − ν + − −ν − ν 



−ν θ + + θ − − ν − 



= θ + − ν + + −ν − ν



( ) ( ) ( ( ) ( ) )

3

22 3 32 4 1 sin sin 3 4 2 3 4

f

AA A AF B





+ −ν θ + + θ + − ν + − − ν 





. (50)
24 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Table 10. General parameters.

 [after Yang and others, 1988, p. 4250]

Expression Description Source

c a = − b 2 2 Prolate ellipsoid focus

acratio = + ( - a c) (a c)

coef ab 1

2

= π2

den1 = π8 1( - ν)

Q d = 3 en1

R d = ( - 1 2 ) en ν 1

I coef ac acratio c a = − 1 ( ) + 2 3 2 ( ) ln( )

I coef a c ac acratio c aa = − 1 ( ) + + 3 2 4 5 2 3( ) 2 ( ) ln( )

a R I 11 a = − 2 4 ( ) π Newman and others (2006)

a R I 12 a = −2 4 ( ) + π Newman and others (2006)

a Qa I RI 21 aa a

2 = + −1 Newman and others (2006)

a Q 22 a Iaa aI R Q 2 = − − − ( ) 2 Newman and others (2006)

den2 = + 3 2 λ μ

num a a 2 22 12 = − 3

den a a a a 3 1 = − 1 22 12 21

num a a 3 11 21 = − 3

( )( ) †

2 2 33 P P den num num den =∆ μ − 2

Pressure Newman and others (2006)

P P den num num den * = ∆ ( ) 1 2 2 2 ( ) λ λ + + ( ) μ 3 3

Pressure Newman and others (2006)

2 †

1a bP = −2

Force

( ) ( ) 22 † *

1b bcP P 3 = + −ν 21 2 

Pressure

General Parameters (YANGPAR)

The general parameters necessary to compute the displacement for a prolate spheroid are listed in table 10 in the same order 

as their use in the code. 
Prolate Spheroid Source (Magma Conduit) 25

−10 −5 0 5 10

−0.20

−0.15

−0.10

−0.05

0.00

0.05

0.10

0.15

−0.20

−0.15

−0.10

−0.05

0.00

0.05

0.10

0.15

U, in meters

−10 −5 0 5 10

z=0, θ=90°, φ=0°, A=0.25

X, in kilometers

V, in meters

−10 −5 0 5 10

-0.10

-0.05

0.00

0.15

0.20

0.25

W, in meters

−10 −5 0 5 10

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

−10 −5 0 5 10

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

z=0, θ=90°, φ=0°, A=0.75

X, in kilometers

FEM

Analytical

−10 −5 0 5 10

0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

EXPLANATION

Figure 14. Free-surface deformation caused by the inflation of a vertical prolate spheroid ( φ = ° 0 ) for different values of the 

geometric aspect ratio A. Comparison between the three-dimensional semianalytical solution (analytical; equation 36), and a 

numerical finite element method model (FEM) of the same source. Source parameters are x y 0 0 = = 0 , z0 = 3 000 , meters, 

length of semimajor axis a = 1000 meters, ∆P μ = 0 1. , μ = 9 6. Gigapascals and ν = 0 25. . Because of a mathematical 

singularity at θ = 90 , we used θ = 89 99.  . U: horizontal deformation (x or East component); V: horizontal deformation (y or 

North component); W: vertical deformation (z or Up component).

Verification

We tested the complete derivation by comparing analytical results to three-dimensional FEM solutions for vertical, tilted, 

and rotated prolate spheroids (figs. 14 and 15).

Results for the internal deformation are shown in figure 16. Although there is no significant difference between our 

formulation and that of Yang and others (1988) for shallow depths, our formulation provides a better fit to the FEM solution for 

greater depths (fig. 16, third row).
26 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

men12-7629_fig15

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

1.0

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

1.0

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

-0.6

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.5

1.0

1.5

2.0

2.5

−10 −5 0 5 10

U, in meters

−10 −5 0 5 10

X, in kilometers

V, in meters

FEM

Analytical

−10 −5 0 5 10

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

W, in meters

z=0, θ=30°, φ=60°, A=0.5 z=0, θ=30°, φ=60°, A=0.5 z=0, θ=30°, φ=60°, A=0.5

EXPLANATION

Figure 15. Free-surface deformation caused by the inflation of a tilted ( θ unequal to 0) and rotated (φunequal to 0) prolate 

spheroid. Comparison between the three-dimensional semianalytical solution (analytical; equation 36), and a numerical 

finite element method model (FEM) of the same source. Source parameters are x y 0 0 = = 0 , z0 = 3 000 , meters, length 

of semimajor axis a = 1 000, meters, ∆P μ = 0 1. , μ = 9 6. Gigapascals and ν = 0 25. . U: horizontal deformation (x or East 

component); V: horizontal deformation (y or North component); W: vertical deformation (z or Up component). See table 11 for 

numerical examples.
Prolate Spheroid Source (Magma Conduit) 27

men12-7629_fig16

Yang (1988)

FEM

Analytical

U, in meters

V, in meters

W, in meters

−10 −5 0 5 10 −10 −5 0 5 10

X, in kilometers

−10 −5 0 5 10

z=0, θ=30°, φ=60°, A=0.5

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

1.0

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

1.0

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

1.0

1.2

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

-0.6

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

-0.6

−0.4

−0.2

1.2

−0.6

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.5

1.0

1.5

2.0

2.5

0.0

0.5

1.0

1.5

2.0

2.5

0.0

0.5

1.5

2.5

3.5

1.0

2.0

3.0

4.0

4.5

EXPLANATION

Figure 16. Internal deformation caused by the inflation of a tilted and rotated prolate spheroid. Comparison between the 

three-dimensional semianalytical solution described in this publication (analytical; equation 36), the original formulation by 

Yang and others (1988), and a numerical finite element method model (FEM) of the same source. Source parameters are 

x y 0 0 = = 0, z0 = 3 000 , meters, length of semimajor axis a = 1 000, meters, ∆P μ = 0 1. , μ = 9 6. Gigapascals and 

ν = 0 25. . First row—displacements at the free surface; second row—displacements at shallow depths (z equals 150 meters); 

third row—displacements at greater depths (z equals 1,500 meters). U: horizontal deformation (x or East component); V: 

horizontal deformation (y or North component); W: vertical deformation (z or Up component). See table 11 for numerical 

examples.
28 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Volume Change

The volume change of the ellipsoidal cavity is given by

∆ ∆

∆

V V P p

P

T

= −

( ) + − 















1 2

2 1

3 ν

ν μ , (51)

where V is the volume of the ellipsoidal cavity, and p

T

 is the trace of the stress inside the 

ellipsoidal cavity (Amoruso and Crescentini, 2009). According to Amoruso and Crescentini 

(2009), the formula obtained by Tiampo and others (2000) for the volume change of an 

ellipsoid, ∆ ∆ V P = 3 4 V ( ) μ , was derived by using a formula for strain which is valid for 

spheres only. As a consequence, Tiampo and others’ (2000) expression is correct for spheres but 

incorrect for ellipsoids and always underestimates the volume change for an ellipsoid cavity.

The rigorous expression 51 can be approximated by the following empirical formulation 

for the volume change (fig. 17):

∆ ∆ V V P A = − A+















 μ

2

3

07 13 . . 7 . (52)

The constants in equation 52 are determined by the polynomial best fitted to the numerical 

values from Amoruso and Crescentini (2009) and several FEM models (fig. 17).

Figure 17. Verification of the empirical correction 

(equation 52) for the volume change of a pressurized 

prolate spheroid. ∆ ∆ V V( ) P μ is the dimensionless 

volume change; 1/A, the reciprocal of the aspect ratio.

1.0

1.1

1.2

1.3

1.4

1 2 3 4 5 6 7 8 9 10 11

A-1

Empirical correction

Amoruso and Crescentini (2009)

Finite Element Method solution

∆V/(V∆P/μ)

EXPLANATION
Prolate Spheroid Source (Magma Conduit) 29

Ground Tilt

Because of the complexity of Yang and others’ (1988) expressions for the displacement, the ground tilt is computed by 

using finite-difference derivatives (see appendix 1). Figure 18 below shows the verification of ground-tilt calculations.

−10 −5 0 5 10

−4

−3

−2

−1

0

1

2

3

4

X, in kilometers

Tilt (east), in radians X 10-4

−10 −5 0 5 10

−4

−3

−2

−1

0

1

2

3

4

5

FEM

Analytical

Tilt (north), in radians X 10-4

z=0, θ=60°, φ=30°, A=0.5

EXPLANATION

Figure 18. Ground tilt caused by the inflation of a tilted and rotated prolate spheroid. 

Comparison between the three-dimensional semianalytical solution described in this publication 

(analytical; equation 36 and appendix 1), and a numerical finite element method model (FEM) of 

the same source. Source parameters are x y 0 0 = = 0, z0 = 3 000 , meters, length of semimajor 

axis a = 1 000, meters, ∆P μ = 0 1. , μ = 9 6. Gigapascals and ν = 0 25. . See table 11 for numerical 

examples.
30 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Strain

Finite-difference derivatives (see appendix 1) are used to compute the components of the displacement-gradient tensor 

(equation 32), verified in figure 19.

Figure 19. Internal strain caused by the inflation of a tilted and rotated prolate spheroid. Comparison between the threedimensional semianalytical solution described in this publication (analytical; equation 36 and appendix 1), and a numerical finite 

element method model (FEM) of the same source. εa is the areal strain, γ1 and γ2 the shear strains. Source parameters are 

x y 0 0 = = 0, z0 = 3 000 , meters, length of semimajor axis a = 1 000, meters, ∆P μ = 0 1. , μ = 9 6. Gigapascals and ν = 0 25. . See 

table 12 for numerical examples.

−10 −5 0 5 10

0

5

10

15

−10 −5 0 5 10

−14

−12

−10

−8

−6

−4

−2

0

2

X, in kilometers

FEM

Analytical

−10 −5 0 5 10

−5

−4

−3

−2

−1

0

1

2

z=150, θ=30°, φ=60°, A=0.5

EXPLANATION

εa, X 10-4

γ1, X 10-5

γ2, X 10-4

Examples of Parameters Calculated for a Prolate Spheroid

Table 11. Surface deformation and ground tilt.

[See equation 36 and appendix 1. Parameters are x0 = 0 meters, y0 = 0 meters, z0 = 3, 000 meters, a = 1 000, meters, A = 0 5. , ∆P μ = 0 1. , μ = 9 6.

GigapascalPa, ν = 0 25. , θ = 30, φ = 60; x and y, local Cartesian coordinates; u u x y , and uz, displacements at the free surface ∂ ∂ u x z and ∂ ∂ u y z , tilt 

components]

x 

(103

 meters)

y

(103 meters) 

ux

(meters)

uy

(meters)

uz

(meters)

duz

/dx

(X 10-3)

duz

/dy

(X 10-3)

-4.0 -4.0 -0.1591 -0.1736 0.1427 0.0368 0.0429

0.0 0.0 0.0625 0.0361 2.3388 0.2229 0.3862

2.0 2.0 0.6528 0.6737 1.0415 -0.3015 -0.2824

Table 12. Internal deformation and strain.

[See equation 36 and appendix 1. Parameters are x0 = 0 meters, y0 = 0 meters, z0 = 3 000 , meters, a = 1 000, meters, A = 0 5. , ∆P μ = 0 1. , μ = 9 6.

GigapascalPa, ν = 0 25. , θ = 30, φ = 60; x and y, local Cartesian coordinates; u u x y , and uz, internal deformation; ε γ a , 1, and γ 2 , strain]

x 

(103

 meters)

y

(103 meters)

z

(103 meters)

ux

(meters)

uy

(meters)

uz

(meters)

ea

(X 10-3)

γ1

(X 10-3)

γ2

(X 10-3)

-4.0 -4.0 1.5 -0.1464 -0.1597 0.1652 -0.0041 -0.0163 -0.0922

0.0 0.0 1.5 0.5138 0.2966 3.7764 4.1151 -0.4081 -0.7067

2.0 2.0 1.5 0.5923 0.6309 0.8141 0.0412 -0.0226 -0.6365
Sill-Like Source 31

A simple 3–D model of a sill intrusion is a horizontal penny-shaped crack in a semiinfinite elastic body (Fialko and others, 2001). The analytical expressions of Fialko and others 

(2001) are appropriate for a sill-like source with a radius as much as five times larger than its 

depth.

We consider a horizontal penny-shaped crack with radius a and depth z0 in an elastic halfspace (fig. 20A). The vertical axis in figure 20A is the axis of symmetry (positive downwards), 

with an origin at the crack center. The surface of the half-space is assumed to be stress-free; r

and z are spatial coordinates, and the index in parentheses corresponds to the upper and lower 

domains of the crack (i=1 and 2 in fig. 20B). The boundary conditions for the problem are 

shown in figure 20A. Figure 21 shows the implementation of the model in the Matlab function 

fialko01.m

We translate and normalize all variables having dimensions of length with respect to 

the crack radius a and all variables having dimensions of pressure with respect to the shear 

modulus μ of the elastic half-space:

x

x x

a y

y y

a r x y

z z z

a h za

= − = − = +

= − =

0 0 2 2

0 0

, , ,

and

. (53)

The general solution for the vertical and horizontal deformation of a half-space is given by 

equations 12 and 13 of Fialko and others (2001):

U B z h A z h z

( ) sinh 1

0

1 2

2 1

= − 

( ) ( ) − + ( ) ( ) ( ) + 

+ −

∞

∫ ν ξ ξ

 ( ) ( ) ν − + ( ) ( ) ( ) + 

 ( )

= − ( ) + + ( )

A z h B z h J r

U A z h r

ξ ξ ξ ξ

ν ξ

cosh

( )

0

1 1 2

d 

B z h

B z h A



( ) ( ) ( ) + 

+ − ( ) ( ) + + ( )

∞

∫ sinh

cosh

ξ

ν ξ ξ

0

 2 1 ( ) ( ) z h +  J r  1 ( ) ξ ξd 

. (54)

All of the variables in equation 54 above are nondimensional, unless noted otherwise. The 

original equations 12 and 13 of Fialko and others (2001) have been changed in equation 54 as 

follows:

3. 2 1 Uz

( ) and 2 1 Ur( ) must be replaced by Uz( )1 and Ur( )1 , respectively, and 

4. d h ξ ξ sinh ( ) must be replaced by dx.

The functions A and B in equation 54 can be defined through the functions Φ and Ψ as 

follows:

A h h

B h h

h

h

= + ( ) ( ) +

= − ( ) ( ) −

−

−

e

e

ξ

ξ

ξ ξ

ξ ξ

Ψ Φ

Ψ Φ

1

1

. (55)

Sill-Like Source
32 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

The expression ( ) 1 e − −2ξh in the original equation 24 of Fialko and others (2001) should be 

replaced by e−ξh . Finally, the functions Φ and Ψ can be expressed in terms of finite Fourier 

transforms:

Φ

Ψ

ξ ξ φ

ξ ξ

ξ ξ ψ

( ) = ( ) ( )

( ) = ( ) − ( ) 



 



 ( )

∫

∫

sin t t dt

t dt

0

1

0

1 sin t

t

cos t

. (56)

The integrals in equation 56 can be solved numerically (for example, by using the GaussLegendre quadrature method; Press and others, 1992, section 4.5). For a hydrostatic pressure 

distribution inside a penny-shaped crack with the boundary conditions given in figure 20A, the 

functions φ and y are the solutions of the following Fredholm equations of the second kind:

φ

π π

φ ψ

ψ

π

t

t T t r r T t r r dr

t T t r

( ) = − +  ( ) ( ) + ( ) ( )  



( ) = ( )

∫ 2 2

2

1 3

0

1

4

, ,

 , φ ψ ( ) r T+ ( ) t r ( ) r dr   ∫ 2 

0

1

,

. (57)

where

T t r h

t r

P t r P t r

T h

tr

P t r P t r

1

3

1 1

2 2 2

4,

,

( ) = −

( )

 ( ) − ( ) +  



= −  ( ) − + ( )  

 + −  ( ) + + ( )  

( ) = − ( ) − ( ) + − − ( )

h P t r P t r

T t r t r t r

h

r P P r t r

3 3

3

2

4 4 , 2 P t r t r P t r

T T t r T

1 1

4 3

 ( ) − + + ( ) ( ) +   { }

( ) , =

. (58)

and

P s h s

h s

P s h s h s h s

h s

1

2 2

2 2 3

2

2 2

4 2 2 4

2 2

12

4

4 8 2

4

( ) = −

( ) +

( ) = + ( ) +

+ −

+

log

( )

( ) = − +

( ) +

( ) = −

( ) +

2

3

4 2 2 4

2 2 3

4

2 2

2 2 2

2 8 2

4

4

4

P s h s h s

h s

P s h s

h s

. (59)

The Fredholm equations 57, 58, and 59 can be solved numerically by using the Nystrom routine 

with the N-point Gauss-Legendre rule (Press and others, 1992, section 18.1). 

Finally, solutions of the dimensionless equations 54 give the following expressions for the 

deformation (fig. 22; see tables 13 and 14 for numerical examples):

u a

PU x

r

u a

PU y

r

u a

PU

x r

y r

z z

 

 

 

=

=

= −

∆

∆

∆

μ

μ

μ

( )

( )

( )

; 1

1

1

. (60)
Sill-Like Source 33

Figure 20. Definition of boundary conditions and geometry for a sill-like source (after Fialko 

and others, 2001). A, Dimensional coordinates and boundary-values. x, y, and z, local Cartesian 

coordinates; a, source radius; the normal stress σn and shear stress στ are equal to zero at the 

free surface; the source center is x y z 0 0 0 ( ) , , , and the free surface is at z = 0. B, Dimensionless 

coordinates and boundary-values. r and z, local cylindrical coordinates; the source radius is equal 

to 1; the normal stress σn and shear stress στ are equal to zero at the free surface; the source 

center is at the origin of the coordinate system, and the free surface is at z h = − .

men12-7629_fig20

x (east)

y (north)

z

(x0 y0 0)

σn = σt = 0

z0

A B

(1)

(2)

a

z

(0,0,0)

σn = σt = 0

0 1

0 -h

r

0 = free surface -h = free surface

(0,0,-h)

¤ ¤

(x0 y0 z0)

FIALKO01

3D displacement, tilt, and strain

FIALKODISP

displacement and

volume change

GAULEG

Gauss-Legendre quadrature PSI_PHI

ψ and φ functions

LEGPOL

Legendre’s polynomial

GAULEG

Gauss-Legendre quadrature T function

P function

Figure 21. Matlab function fialko01.m flow chart for calculating and coding the parameters for a 

sill-like source. 
34 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Figure 22. Graphs showing (top) surface deformation due to the inflation of a sill-like (or penny-crack) source; (bottom) internal 

deformation. Comparison between the three-dimensional semianalytical solution (analytical; equation 60), and a numerical finite 

element method model (FEM) of the same source. Source parameters are location x y 0 0 = = 0 , depth z0 = 1 000 , meters, 

radius a = 1 000, meters, ∆P μ = 0 01. , μ = 9 6. Gigapascals, and ν = 0 25. . U: horizontal deformation (x or East component); 

V: horizontal deformation (y or North component); W: vertical deformation (z or Up component). See table 13 for numerical 

examples.

men12-7629_fig22

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.1

0

0.05

0.10

0.15

0.20

0 0

0.05

0.40

0.10

0.15

0.20

0.25

0.30

0.35

0 2 4 6 8 10

U, in meters

0 2 4 6 8 10

z=0, a=1,000m, z0=1,000m 

Radial distance, in kilometers

V, in meters

0 2 4 6 8 10

0

1

2

3

4

5

6

7

W, in meters

0 5 10 15 0 5 10 15

z=500m, a=1,000m, z0=1,000m 

0 5 10 15

0

1

2

3

4

5

6

7

FEM

Analytical

EXPLANATION
Sill-Like Source 35

Volume Change

For a uniformly pressurized crack, the expression for the volume change (fig. 23) is

∆ ∆ V P = − ( ) − a t t dt ∫ 4 1 3

0

1

π ν

μ

φ( ) (61)

with φ( )t given by equation 57.

Ground Tilt

Because of the complexity of Fialko and others’ (2001) expressions for the displacement, 

the ground tilt is computed by using finite-difference derivatives (fig. 24, appendix 1).

Figure 23. Volume change for a sill-like source

(equation 61). ∆ ∆ V a( ) P 3 μ , ∆ is the dimensionless 

volume change; z a 0 the dimensionless source 

depth. The dimensionless pressure change is 

∆P μ = 0 01. .

0 2 4 6 8 10

−2.5

−2.0

−1.5

−1.0

−0.5

0.0

0.5

X, in kilometers

Tilt (north), in radians × 10-3

0 2 4 6 8 10

−5

−4

−3

−2

−1

0

1

Y, in kilometers

Tilt (east), in radians × 10-3

FEM

Analytical

EXPLANATION

Figure 24. Ground tilt caused by the inflation of a sill-like source. Comparison 

between the three-dimensional semianalytical solution (analytical), and a 

numerical finite element method model (FEM) of the same source. Source 

parameters are location x y 0 0 = = 0 , depth z0 = 1 000 , meters, radius 

a = 1 000, meters, ∆P μ = 0 01. , μ = 9 6. Gigapascals, and ν = 0 25. . See 

table 13 for numerical examples.

0

2

4

6

8

10

12

14

0 1 2 3 4

Volume change from Fialko and others (2001)

Volume change from the Finite Element Method 

numerical model

EXPLANATION

∆V/(a3∆P/μ)

z0 /a
36 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Strain

The components of the displacement gradient tensor are computed by finite-difference derivatives (fig. 25, appendix 1).

0 2 4 6 8 10

−5

0

5

10

15

εa, × 10-4

0 2 4 6 8 10

−2

0

2

4

6

8

Radial distance, in kilometers

γ1, × 10-4

0 2 4 6 8 10

−10

−8

−6

−4

−2

0

2

γ2, × 10-4

FEM

Analytical

EXPLANATION

Figure 25. Internal strain caused by the inflation of a sill-like source. Comparison between the three-dimensional 

semianalytical solution described in this publication (analytical), and a numerical finite element method model (FEM) of 

the same source. εa is the areal strain, γ1 and γ2 the shear strains. Source parameters are location x y 0 0 = = 0 , depth 

z0 = 1 000 , meters, radius a = 1 000, meters, ∆P μ = 0 01. , μ = 9 6. Gigapascals. See table 14 for numerical examples.

Examples of Parameters Calculated for a Sill-Like Structure

Table 13. Surface deformation and tilt.

[See equation 22 and appendix 1. Parameters are x0 = 0 meters, y0 = 0 meters, z0 = 1, 000 meters, a = 1 000, meters, ∆P μ = 0 0. 1, μ = 9 6. GPa, ν = 0 25. ; r, 

radial distance from source center; u u x y , , and uz, displacements at the free surface; ∂ ∂ u x z and ∂ ∂ u y z , tilt components] 

r 

(103

 meters)

ux

(meters)

uy

(meters)

uz

(meters)

duz/dx

(X 10-3)

duz

/dy

(X 10-3)

0.0000 0.0000 0.0000 6.4302 0.0000 0.0000

1.0062 0.7493 1.4985 2.6405 -2.1202 -4.2405

2.0125 0.2207 0.4415 0.2980 -0.2867 -0.5734

Table 14. Internal deformation and strain. 

[See equation 22 and appendix 1. Parameters are x0 = 0 meters, y0 = 0 meters, z0 = 1, 000 meters, a = 1 000, meters, ∆P μ = 0 0. 1, μ = 9 6. GPa, ν = 0 25. ; r, 

radial distance from source center; u u x y , , and uz, and internal deformation; ε γ a , 1 and γ 2 , strain] 

r 

(103

 meters)

z

(103 meters)

ux

(meters)

uy

(meters)

uz

(meters)

ea

(X 10-3)

γ1

(X 10-3)

γ2

(X 10-3)

0.0000 0.5000 0.0000 0.0000 7.2591 1.2431 0.0000 0.0000

1.0062 0.5000 0.1759 0.3518 2.6398 -0.2045 0.5917 -0.7890

2.0125 0.5000 0.0415 0.0830 0.2585 0.0229 0.0416 -0.0554
Surface Deformation and Ground Tilt for Rectangular Dikes and Faults 37

Surface Deformation and Ground Tilt for Rectangular Dikes and Faults

The solutions proposed by Okada (1985) present a complete suite of closed analytical expressions (fig. 26) for the surface 

displacements and tilts caused by a strike-slip, dip-slip or tensile rectangular dislocation in a half-space that is homogenous, 

isotropic, flat, and elastic. These expressions are particularly compact and free from mathematical singularities 

The rectangular source geometry is illustrated in the Cartesian coordinate system shown in figure 27A. The points 

( , x y ) i i and ( , x y ) f f represent the initial and end coordinates of the source trace at the surface. The elementary dislocations 

U correspond to the strike-slip (U1), dip-slip (U2), or tensile (U3) component of an arbitrary dislocation. In the convention 

used by Okada (1985), coordinates are assumed to be in kilometers, whereas U and the displacements at the Earth’s surface, 

u = [ , u u , ] u x y z , are in meters. 

Figure 26. Program flow chart for rectangular dislocations (Okada, 1985). 

Figure 27. A, Geometry of the dislocation model in a Cartesian coordinate system; B, Geometry of the source model (after Okada, 

1985). In Okada (1985), input coordinates are assumed to be in kilometers, whereas the dislocations and displacements are in meters. 

Strike-slip dislocation U1

 greater than 0 identifies a right lateral motion, dip-slip dislocation U2 greater than 0 identifies a reverse 

motion, and tensile dislocation: U3 greater than 0 identifies a tensile opening. 

OKADA85

surface displacement and tilt

OK8525

strike-slip dislocation

1. general parameters

2. elementary displacement

3. displacement gradient

OK8526

dip-slip dislocation

1. general parameters

2. elementary displacement

3. displacement gradient

OK8527

tensile dislocation

1. general parameters

2. elementary displacement

3. displacement gradient

z (up)

y (north)

x (east)

(xi,yi)

δ

(xf,yf)

φ

zt

zb L

W

A B

δ

z yp

xp

z=-d

U3

L

W

U1

U2

(x0,y0,z0)
38 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Displacement

The angle φ is the strike of the fault

p y d

q y d

p

p

=

= −

cos sin

cos

δ δ

δ δ

+

 sin , (62)

and δ is the dip angle. The length L and width W of the rectangular dislocation are

L x x y y

W z z

f i f i

b t

= − ( ) + − ( )

= −

2 2

 and

sin . δ

(63)

and the coordinates of the lower corner of the dislocation are (fig. 27A)

x x z z

y y z z

z z

i b t

i b t

b

0

0

0

= + ( ) −

= − ( ) −

=

cot cos

cot sin

δ φ

δ φ . (64)

Elementary expressions for surface deformation and the displacement-gradient matrix can 

be written as functions of nine sets of variables: 

1. the fault type, 

2. the coordinates of each dislocation corner (as given by xp, L, and W), 

3. x x x y y y z d z n n = − 0 0 = − n =− =− 0 , (65)

4. sin cos

 cos sin

pn n

p nn

xxy

y xy

= φ⋅ + φ⋅

= − φ⋅ + φ⋅

, (66)

5. the depth in kilometers of the lower left corner of the rectangular fault (d), 

6. the Lame’s elastic parameters μ and λ,

7. the dislocation U, 

8. the dip angle δ, and 

9. the rotated coordinates p and q

p y d

q y d

p

p

=

= −

cos sin

cos

δ δ

δ δ

+

 sin

. (67)

To compute the displacement, the function okada85.m follows a three-step procedure:
Surface Deformation and Ground Tilt for Rectangular Dikes and Faults 39

1. Compute the three components of the displacement and the six components of the displacement gradient matrix at the four 

corners of the dislocation:

( ) ( )

( ) ( )

1

2

3

 ,,,,, , , ,,,, ,

 ,,,,, , , ,,,, ,

 ,,,

i i ip i

i i ip i

i

uuvvww

uvw f q U f x pq U

xyxyx y

uuvvww

uvw f q U f x p Wq U

xyxyx y

uuvvww

uvw f q

xyxyx y

  ∂∂∂∂∂ ∂ = xη μλ δ = μλ δ     ∂∂∂∂∂ ∂

  ∂∂∂∂∂ ∂ = xη μλ δ = − μλ δ     ∂∂∂∂∂ ∂

  ∂∂∂∂∂ ∂   = xη μ

  ∂∂∂∂∂ ∂ ( ) ( )

( ) ( )

4

,, , , ,,,, ,

 ,,,,, , , ,,,, ,

i ip i

i i ip i

U f x Lpq U

uuvvww

uvw f q U f x Lp Wq U

xyxyx y

λ δ = − μλ δ

  ∂∂∂∂∂ ∂ = xη μλ δ = − − μλ δ     ∂∂∂∂∂ ∂

, (68)

where 

 (u, v, w) is the displacement, [ ]  k identifies the corners of the dislocation ( ) k = 1, , 234, , 

 the subscript i refers to the three different types of dislocation (1=strike, 2=dip, 3=tensile),

 x and η are dummy variables of integration.

2. Combine the displacements and displacement gradients from step 1 to compute the deformation, equation 69, and tilt, 

equation 70, in Okada’s (1985) reference system (fig. 27B).

[ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ]

1234

1234

1234

p

x

p

y

p

z

U uu u u

U vv v v

U ww w w

=−−+

=−−+

=−−+

, (69)

and

[ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ]

1 2 3 4

1 2 3 4

p

z

p

z

U ww w w

xxx x x

U ww w w

yyy y y

∂ ∂∂ ∂ ∂

=−−+

∂∂∂ ∂ ∂

∂ ∂∂ ∂ ∂

=−−+

∂∂∂ ∂ ∂

. (70)

3. Rotate the deformation and tilts back to the original Cartesian reference system (fig. 27A):

East

North

Up

u U U

u U U

u U

x x

p

y

p

y y

p

x

p

z z

p

= ⋅ − ⋅

= ⋅ + ⋅

=

sin cos

sin cos

φ φ

φ φ , (71) 

and 

tilt (East)

tilt (N

∂

∂ = ∂

∂ − ∂

∂

















u

x

U

x

U

y

z z

p

z

p

0 001 . sin φ φ cos

orth) ∂

∂ = ∂

∂

+

∂

∂

















u

y

U

y

U

x

z z

p

z

p

0 001 . sin φ φ cos

. (72)

The scale 0.001 in equation 72 takes into account that, in Okada (1985), coordinates are in kilometers and displacements in 

meters. 
40 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Parameters for Displacements and Displacement Gradients Functions fi

The parameters and variables, which are listed below in the order they are coded, are the 

same for all the functions fi . 

1. Elastic parameter

 12 μ α= = − ν

λ+μ

; (73)

2. Geometric transformations, equation 30 from Okada (1985), where x and η are dummy 

integration variables whose upper and lower limits are defined in equation 68

2222

222

cos sin

sin cos

y q

d q

R yd

X q

=η δ+ δ

=η δ− δ

=x + +

=x +





 

(74)

3. Equation 28 (Okada,1985) 

( ) ( )

( ) ( )

( )

( )

( ) ( )

4

1

5

3 4

1 5

2 3

1 ln sin ln

cos

2 cos sin

tan

cos cos

1 ln tan

cos

1 tan

cos

ln

I Rd R

X q XR X I

R X

y I RI

R d

I I

R d

I RI

−

=α + − δ +η     δ

  η + δ+ + δ

= α   δ x+ δ  

 

= α − +η + δ     δ +

  x =α − − δ     δ +

=α − +η −









; (75)

when ξ = 0 (Okada, 1985, singularity condition ii, p. 1148), I5 = 0;

when R + = η 0 (Okada, 1985, singularity condition iii, p. 1148), replace ln ( ) R + η

with − − ln ( ) R η in equation 74;

if cos =δ 0 , use equation 29 (Okada,1985); 

( )

( ) ( )

5

1

2

4

2

3

0.5

0.

 

 

5 ln R

sin

q I

R d

yq I

R d R d

q

R d

R

I

d

I

x =− α

+

 

  η = α + − +η

+ +  

α

+

x

= −

− δα

+ =





 





(76) 

4. Equation 40 (Okada,1985) 
Surface Deformation and Ground Tilt for Rectangular Dikes and Faults 41

( )

( )

3

1

2 3

1 1

cos

1 sin 1

cos

sin cos 1

q y K

R R RR d

K

RR d R R

q K K

R RR

 

= α −  

δ +η +  

  x δ = α −  

δ + + η  

  δ δ

=α − + −     + η







; (77)

if cos 0 δ = , use equation 41 from Okada (1985): 

( )

( )

1 2

2

3 1

q K

RR d

sin K

R d RR d

αx =

+

  δ x = α −  

+   +  



 

(78)

5. Equation 34 (Okada 1985) 

( )

( )

3

3

2

1 2

2 1

2

4

2

1

1 1

cos

1

tan

tan

1 

cos sin 1 

cos

K

J J

R R

q J J

R R

J

R d RR d

y J K

R

RR d

 

  x =α − −

δ + +  

x = α −

δ +

α

δ

δ

  x =− −     + η

  δ δ =− − −     + η

α

 



 ; (79)

if cos 0 δ = , use equation 35 from Okada,1985) 

( ) ( )

( ) ( )

2

1 2

2

2 2

2 0.5 1

sin 2 0.5 1

q J

R d RR d

q J

R d RR d

  x = α −  

+   +  

  x δ = α −  

+   +  

 

 

(80)

6. Equation 36 (Okada,1985) 

( ) ( ) 3 3 2 2

2 2 c nR R A A

RR RR

+ x + η = = + x + η

(81)

7. If R + = η 0 (Okada, 1985, singularity condition iii, p. 1148), replace ln ( ) R + η with

− − ln ( ) R η and set 1 0

R + = η .
42 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Elementary Function f

1 for a Strike-Slip Fault

Elementary Displacement (Okada, 1985, equation 25)

u

U q

R R qR

I

v

U yq

R R

= − ++





 



 +





 





= − +

1 −1

1

1

2

1

2

1

π

ξ

η

ξη δ

π η

tan sin

 +

+

+





 





= − ++

+

+





q RI

w

U dq

R R q R

I

cos sin

sin sin

δ

η

δ

π η

δ

η

δ

1

2

1 1

2

1

4



 





; (82)

if q = 0 (Okada, 1985, singularity condition i, p. 1148), then

1

1

1 sin

2

U q u I

R R

  x =− + δ   π +η  

. (83)

Displacement Gradient (Okada, 1985, equation 37)

1 2

2

1

3 1

2

3 3

cos sin

sin

cos cos si

2

n

2

n

n

w q q A K

x R

w dq yq qA K

y R R

U

R

U  

  π  

π

∂ x   = −x δ + − δ   ∂  

∂ δ     = δ+ x δ− + − δ     ∂    

 

(84)

Elementary Function f2 for a Dip-Slip Fault

Elementary Displacement (Okada, 1985, equation 26)

2

3

2 1

1

2 1

5

sin cos

2

1 cos tan sin cos

2

1 sin tan sin cos

2

U q u I

R

U yq v I

R R qR

U dq w I

R R qR

−

−

  =− − δ δ   π  

  xη =− + δ − δ δ   π +x  

  xη =− + δ − δ δ   π +x  





; (85)

if q = 0 (Okada, 1985, singularity condition i, p. 1148), then

2

1

2

5

1 sin cos

2

1 sin cos

2

U yq v I

R R

U dq w I

R R

  =− − δ δ   π +x  

  =− − δ δ   π +x  





(86)
Surface Deformation and Ground Tilt for Rectangular Dikes and Faults 43

Displacement Gradient (Okada, 1985, equation 38)

3

2

3

2

1

2

2 1 sin sin sin co

sin 1 sin o

s

s

1

2

c

c

U

U d ydqA

w dq q

K

R R

K

x R R R

w

y R R

∂ δ   = + + δδ   ∂   + η

∂ δ

π

    x = − +     δ +

π

δ δ

∂   + x   +   η

  



(87)

Elementary Function f3 for a Tensile Crack

Elementary Displacement (Okada, 1985, equation 27)

2

3 2

3

3 1 2

1

3 1 2

5

sin

2

1 sin tan sin

2

1

cos tan sin

1

2

1

1

U q u I

U dq q v I

R R

R

R

R R r qr

U yq q w I

R R r qr

−

−

  = −δ  

π  

    x   xη = − −δ − − δ       π +x      

    x   xη = +δ − − δ    

+ η

+

  π +x       

η

+ η 





; (88)

if q = 0 (Okada, 1985, singularity condition i, p. 1148), then

v

U dq

R R

q

r I

w

U yq

R

R = − + − 



 



 − 



 





=

+

3

1

2

3

2

1

2

1

1

π ξ

δ ξ δ

π

η





sin sin

R

q

r I+ R

+





 



 − 



 



 ξ +

δ ξ δ

η

cos sin 1

5

2

(89)

Displacement Gradient (Okada, 1985, equation 38)

∂

∂ = − − +

















∂

∂ = −

w

x

q

R

q A K

w

y

y

U

U

n 3

3

2

3

3

3

2

2

2

π

δ

π

δ δ sin

cos sin

sin δ δ ξ δ δ

ξ

+ + − δ

+ − 



 



 



 



 

d q A q A q

R R

cos s c n in cos s K in 2 2

1

2 1 2

. (90)
44 Modeling Crustal Deformation near Active Faults and Volcanic Center—A Catalog of Deformation Models

Verification

Matlab solutions (function okada85.m) based on the analytical expressions by Okada (1985) for strike-slip (fig. 28), dip-slip 

(fig. 29) and tensile (figs. 30 and 31) dislocations are compared to solutions from COULOMB 3.33 and numerical FEM models. 

Strike-Slip Fault

Figure 28. Vector (top) and profile plots (bottom) of the displacement field of a strike-slip fault (thick red line). The profile 

was taken along the green line. Okada85, Matlab script of Okada’s (1985) expressions; Coulomb, solution from COULOMB 3.3. 

Fault parameters are xi

 = 4,000 meters, yi = 4,000 meters, xf = 25,000 meters, yf = 20,000 meters, δ =60°, μ = 1 Gigapascal, 

v = 0.25, U = 1 meter, zt

 =1000 meters, and zb = 10,000 meters. See table 15 for numerical examples.

−40 −20 0 20 40

Okada85

X, in kilometers

Y, in kilometers

−40 −20 0 20 40

−50

−40

−30

−20

−10

0

10

20

30

40

−50

−40

−30

−20

−10

0

10

20

30

40

50 50

Coulomb

−50 0 50

−0.10

−0.08

−0.06

−0.04

−0.02

0.00

0.02

0.04

U, in meters − east

−50 0 50

−0.25

−0.20

−0.15

−0.00

−0.05

0.00

0.05

0.10

V, in meters − north

Coulomb

Okada85

−50 0 50

−0.02

0.00

0.02

0.04

0.06

0.08

0.10

0.12

W, in meters − up

X, in kilometers

EXPLANATION

 3COULOMB, a free software package developed by the USGS Earthquake Hazards Program, covers many of the same deformation sources as our Matlab 

codes, is based on Okada’s (1992) equations (earthquake.usgs.gov/research/modeling/coulomb/overview.php), and is widely used.
Surface Deformation and Ground Tilt for Rectangular Dikes and Faults 45

Dip-Slip Fault

men12-7629fig29

EXPLANATION

Okada85

X, in kilometers

Y, in kilometer

Coulomb

−50 0 50

−0.08

−0.07

−0.06

−0.05

−0.04

−0.03

−0.02

−0.01

0.00

0.01

U (east), in meters

−50 0 50

−0.10

−0.08

−0.06

−0.04

−0.02

0.00

0.02

0.04

V (north), in meters

X, in kilometers

−50 0 50

−0.02

0

0.02

0.04

0.06

0.08

0.

0.12

0.14

0.16

W (up), in meters

Coulomb

Okada85

−40 −20 0 20 40 −40 −20 0 20 40

−50

−40

−30

−20

−10

0

10

20

30

40

50

−50

−40

−30

−20

−10

0

10

20

30

40

50

Figure 29. Vector (top) and profile plots (bottom) of the displacement field of a dip-slip fault (thick red line). The profile was 

taken along the green line. Okada85, Matlab script of Okada’s (1985) expressions; Coulomb, solution from COULOMB 3.3. 

Fault parameters are: xi

 = 4,000 meters, yi = 4,000 meters, xf = 25,000 meters; yf = 20,000 meters; δ =60°; μ = 1 GPa; v = 0.25; U = 

1 meters; zt =1,000 meters; zb = 10,000 meters. See table 16 for numerical examples.
46 Modeling Crustal Deformation near Active Faults and Volcanic Center—A Catalog of Deformation Models

men127629fig30

EXPLANATION

−10 −5 0 5 10

−10

−5

0

5

10

Okada85

X, in kilometers

Y, in kilometers

−10 −5 0 5 10

−10

−5

0

5

10

Finte Element Method numerical solution

−40 −20 0 20 40

−0.05

0.00

0.05

U (east), in meters

−40 −20 0 20 40

−0.4

−0.3

−0.2

−0.1

0.0

0.1

0.2

0.3

V (north), in meters

X, in kilometers

−40 −20 0 20 40

−0.10

−0.05

0.00

0.05

0.10

0.15

w (up), in meters

FEM

Okada85

Figure 30. Vector (top) and profile plots (bottom) of the displacement field of a tensile crack (thick red line). The profile was 

taken along the green line. Okada85, Matlab script of Okada’s (1985) expressions; FEM, numerical solution from FEM model. 

Fault parameters are: xi

 = -5,000 meters; yi = 0 meters, xf = 5,000 meters; yf = 0 meters; δ =90°; μ = 1 GPa; v = 0.25; U = 1 meters; 

zt

 =1,000 meters; zb = 6,000 meters. See table 17 for numerical examples.

Tensile Crack
Surface Deformation and Ground Tilt for Rectangular Dikes and Faults 47

Figure 31. Profile plots showing tilt caused by a tensile crack. The profile was 

taken along the green line shown in figure 30. Okada85, Matlab script of Okada’s 

(1985) expressions; FEM, numerical solution from FEM model. Fault parameters: 

xi

 = -5,000 meters, yi = 0 meters, xf = 5,000 meters; yf = 0 meters; δ =90°; μ = 1 

GPa; v = 0.25; U = 1 meters; zt

 =1,000 meters; zb = 6,000 meters. See table 17 for 

numerical examples.

men12-7629_fig31

−40 −20 0 20 40

−2.0

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

2.0

Tilt (east), in radians × 10-5

−40 −20 0 20 40

−2.0

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

2.0

2.5

X, in kilometers

Okada85

FEM

Tilt (north), in radians × 10-4

EXPLANATION
48 Modeling Crustal Deformation near Active Faults and Volcanic Center—A Catalog of Deformation Models

Table 17. Tensile crack. 

[Vector and profile graphs shown in figure 30; graphs of tilt profiles shown in figure 31. Parameters are xi = -5,000 meters, yi = 0 meters, xf = 5,000 meters, yf = 0 

meters, δ = 90 degrees, μ= 1 Gigapascal,ν = 0.25, U = 1 meter, zt = 1,000 meters, zb = 6,000 meters; x and y, local coordinates; u u x y , and uz, displacements at 

the free surface; ∂ ∂ u x z and ∂ ∂ u y z , tilt components]

x

(103

 meters)

y

 (103 meters)

ux

(meters)

uy

(meters)

uz

(meters)

duz

/dx

(X 10-5)

duz

/dy

(X 10-5)

 0.2010 0.2010 -0.0019 -0.0001 -0.0609 0.0128 11.4199

 2.6131 2.6131 0.0159 0.2129 0.1378 -1.0750 -1.0556

 5.0251 5.0251 0.0467 0.1428 0.0618 -1.2701 -1.2241

Table 15. Strike-slip fault. 

[Vector and profile graphs shown in figure 28. Parameters are xi = 4,000 meters, yi = 4,000 meters, xf = 25,000 meters, yf = 20,000 meters, δ = 60 degrees, μ= 1 

Gigapascal, ν = 0.25, U = 1 meter, zt = 1,000 meters, zb = 10,000 meters; x and y, local coordinates; u u x y , and uz , displacements at the free surface; ∂ ∂ u x z

and ∂ ∂ u y z , tilt components] 

x

(103

 meters)

y

 (103 meters)

ux

(meters)

uy

(meters)

uz

(meters)

duz

/dx

(X 10-5)

duz

/dy

(X 10-5)

-5.0000 0.0000 0.0274 -0.0141 0.0073 0.0837 -0.0250

0.0000 0.0000 -0.0132 -0.0529 0.0241 0.8793 -0.1957

5.0000 0.0000 -0.0837 -0.2261 0.1190 1.3779 2.0268

Table 16. Dip-slip fault.

[Vector and profile graphs shown in figure 29. Parameters are xi = 4,000 meters, yi = 4,000 meters, xf = 25,000 meters, yf = 20,000 meters, δ = 60 degrees, μ= 

1 Gigapascal, ν = 0.25, U = 1 meter, zt = 1,000 meters, zb = 10,000 meters; km, kilometers; x and y, local coordinates; u u x y , and uz , displacements at the free 

surface; ∂ ∂ u x z and ∂ ∂ u y z , tilt components]

x 

(km)

y

(10 km) 

ux

(meters)

uy

(meters)

uz

(meters)

duz

/dx

(X 10-5)

duz

/dy

(X 10-5)

-5.0000 0.0000 0.0035 -0.0057 -0.0126 -0.0321 -0.1479

0.0000 0.0000 -0.0269 -0.0175 -0.0025 0.7715 -0.4180

5.0000 0.0000 -0.0791 -0.0824 0.1071 3.1341 3.0545

Examples of Surface Deformation and Ground Tilt for Three Types of Faults
Internal Deformation and Strain for Rectangular Dikes and Faults 49

Internal Deformation and Strain for Rectangular Dikes and Faults 

The solutions proposed in Okada (1992) present a complete suite of closed analytical expressions for the internal 

deformation and strain caused by a strike-slip, dip-slip, or tensile rectangular dislocation (fig. 32) in a half-space that is 

homogeneous, isotropic, flat, and elastic. These expressions are particularly compact and free from mathematical singularities. 

Source Geometry

The rectangular source geometry and the Cartesian coordinate system for Okada (1992) are identical to those in the coding 

of Okada (1985, see section, “Surface Displacement and Tilt”; in this report, fig. 27, and equations 62−66).

Displacement

Elementary expressions for internal deformation and the displacement-gradient matrix can be written as functions of the 

fault type (fault), the coordinates of each fault corner (as given by xp, L, and W), the depth in kilometers of the lower left corner 

of the rectangular fault (c), the Lame’s elastic parameters μ and λ, the dislocation U, the dip angle δ, the rotated coordinate y p, 

and the depth of the point of internal deformation z. 

The displacement and displacement gradient matrix are given by the superposition of four different contributions (Okada, 

1992, table 6):

1. infinite medium at depth z (Okada, 1992, p. 1031) designated by A in equation 91

2. the mirror image of uA at depth –z (Okada, 1992, p. 1031) designated by A and ~ in equation 92

3. the free-surface deformation (Okada, 1992, equation 18) designated by B in equation 93

4. a depth-multiplied term (Okada, 1992, p. 1031) designated by C in equation 94

To compute the displacement, follow a four-step procedure:

1. Compute the value of the elementary displacement for each contribution at each fault corner a – d.

2. Combine the values computed in step 1 for each contribution.

3. Compute the deformation in Okada’s (1992) coordinate system (as shown in fig. 27B).

4. Rotate the deformation back to the Cartesian coordinate system (fig. 27A).

Figure 32. Program flow chart for rectangular dislocations (Okada, 1992). 

OKADA92

internal deformation and strain

OK92FA

infinite medium

1. general parameters

2. elementary displacement

3. displacement derivatives

OK92FB

free surface

1. general parameters

2. elementary displacement

3. displacement derivatives

OK92FC

depth multiplied term

1. general parameters

2. elementary displacement

3. displacement derivatives
50 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Step 1—Elementary Displacement and Gradient

In the equations below, u is the displacement; j the x-derivative and k the y-derivative; and 

1, 2, and 3 represent the strike-slip, dip-slip, and tensile components, respectively.

1. For an infinite medium at depth z (equations for f

A are in tables 18–21, p. 58–59), 

( )

( )

1 2 31 23 1 2 3

1 2 31 23 1 2 3

1 2 31 23 1 2 3

 , ,0, , , , , , ,

 , , ,,,, ,, ,

 , ,0, ,

A A AAAA A A A

Ap p a

A A AAAA A A A

Ap p b

A A AAAA A A A

A p c

u u u j j j k k k f fault x c U y z

u u u j j j k k k f fault x W c U y z

u u u j j j k k k f fault x L c

  = μλ δ  

  = μλ δ  

  = −   ( )

1 2 31 23 1 2 3 ( )

,, ,, ,

 , , ,,,, ,, ,

p

A A AAAA A A A

A p p d

U yz

u u u j j j k k k f fault x L W c U y z

μλ δ

  = − μλ δ  

. (91)

2. For an infinite medium at depth –z (equations for f

A are in tables 18–21, p. 58–59),

( )

( )

1 2 31 23 1 2 3

1 2 31 23 1 2 3

1 2 31 23 1

 , ,0, , , , , , ,

 , , ,,,, ,, ,

 

A A AAAA A A A

Ap p a

A A AAAA A A A

Ap p b

A A AAAA A

u u u j j j k k k f fault x c U y z

u u u j j j k k k f fault x W c U y z

uuujjjkk

  = μλ δ −  

  = μλ δ −  

  

  

    ( )

( )

2 3

1 2 31 23 1 2 3

 , ,0, , , , , , ,

 , , ,,,, ,, ,

A A

A p p c

A A AAAA A A A

A p p d

k f fault x L c U y z

u u u j j j k k k f fault x L W c U y z

  = − μλ δ −  

  = − μλ δ −  



  

. (92)

3. The free-surface deformation (equations for f

B are in tables 22–25, p. 60–61)

( )

( )

1 2 31 23 1 2 3

1 2 31 23 1 2 3

1 2 31 23 1 2 3

 , ,0, , , , , , ,

 , , ,,,, ,, ,

 , ,0

B B BBBB B B B

Bp p a

B B BBBB B B B

Bp p b

B B BBBB B B B

B p c

u u u j j j k k k f fault x c U y z

u u u j j j k k k f fault x W c U y z

u u u j j j k k k f fault x L

  = μλ δ  

  = μλ δ  

  = −   ( )

1 2 31 23 1 2 3 ( )

,,,, ,, ,

 , , ,,,, ,, ,

p

B B BBBB B B B

B p p d

c U yz

u u u j j j k k k f fault x L W c U y z

μλ δ

  = − μλ δ  

. (93)

4. The depth-multiplied term (equations for f

C are in tables 26–29, p. 62–64)

( )

( )

1 2 31 23 1 2 3

1 2 31 23 1 2 3

1 2 31 23 1 2 3

 , ,0, , , , , , ,

 , , ,,,, ,, ,

 , ,0

C C CCCC C C C

Cp p a

C C CCCC C C C

Cp p b

C C CCCC C C C

C p c

u u u j j j k k k f fault x c U y z

u u u j j j k k k f fault x W c U y z

u u u j j j k k k f fault x L

  = μλ δ  

  = μλ δ  

  = −   ( )

1 2 31 23 1 2 3 ( )

,,,, ,, ,

 , , ,,,, ,, ,

p

C C CCCC C C C

C p p d

c U yz

u u u j j j k k k f fault x L W c U y z

μλ δ

  = − μλ δ  

. (94)
Internal Deformation and Strain for Rectangular Dikes and Faults 51

Step 2—Combination of Values

1. For, A u , an infinite medium at depth z,

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

AA A A A

abcd

AA A A A

abcd

AA A A A

abcd

uu u u u

uu u u u

uu u u u

= −−+        

= −−+        

= −−+        

. (95)

2. For, A u , an infinite medium at depth –z, 

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

AA A A A

abcd

AA A A A

abcd

AA A A A

abcd

uu u u u

uu u u u

uu u u u

= −−+        

= −−+        

= −−+        

   

   

   

. (96)

3. For, , the free-surface deformation, 

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

BB B B B

abcd

BB B B B

abcd

BB B B B

abcd

uu u u u

uu u u u

uu u u u

= −−+        

= −−+        

= −−+        

. (97)

4. For, C u , the depth-multiplied term,

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

CC C C C

abcd

CC C C C

abcd

CC C C C

abcd

uu u u u

uu u u u

uu u u u

= −−+        

= −−+        

= −−+        

. (98)

Step 3—Computation of Deformation in Okada’s Coordinate System (1992)

( ) ( )

( ) ( )

111 1

222 2 333 3

222 2 333 3

cos sin

sin cos

pAAB C

x

p AAB C AAB C

y

p AAB C AAB C

z

U u u u zu

U u u u zu u u u zu

U u u u zu u u u zu

= − + +⋅

= − + + ⋅ δ− − + + ⋅ δ

= − + − ⋅ δ+ − + − ⋅ δ



 

 

. (99)

Step 4—Rotate Deformation Back to Cartesian Coordinate System

sin cos

sin cos

p p

xx y

p p

yy x

p

z z

uU U

uU U

u U

= φ⋅ − φ⋅

= φ⋅ + φ⋅

=

. (100)
52 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Displacement-Gradient Matrix

To compute the displacement-gradient matrix, follow a procedure similar to the one 

described in the previous section (equations 95–98).

Combination of Values for X-Derivatives 

These equations are based Okada’s equations (1992, table 7).

For, jA, an infinite medium at depth z

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

AA A A A

abcd

AA A A A

abcd

AA A A A

abcd

jj j j j

jj j j j

jj j j j

= −−+        

= −−+        

= −−+        

. (101)

For, jA, an infinite medium at depth –z (mirror image),

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

AA A A A

abcd

AA A A A

abcd

AA A A A

abcd

jj j j j

jj j j j

jj j j j

= −−+        

= −−+        

= −−+        

   

   

   

. (102)

For, jB, the free-surface deformation,

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

BB B B B

abcd

BB B B B

abcd

BB B B B

abcd

jj j j j

jj j j j

jj j j j

= −−+        

= −−+        

= −−+        

. (103)

For, j

C , the depth-multiplied term,

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

CC C C C

abcd

CC C C C

abcd

CC C C C

abcd

jj j j j

jj j j j

jj j j j

= −−+        

= −−+        

= −−+        

. (104)

Combination of Values for Y-Derivatives (Okada 1992, table 8)

For, kA, an infinite medium at depth z,

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

AA A A A

abcd

AA A A A

abcd

AA A A A

abcd

kk k k k

kk k k k

kk k k k

= −−+        

= −−+        

= −−+        

. (105)
Internal Deformation and Strain for Rectangular Dikes and Faults 53

 For, 

kA , an infinite medium at depth –z (mirror image),

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

AA A A A

abcd

AA A A A

abcd

AA A A A

abcd

kk k k k

kk k k k

kk k k k

= −−+        

= −−+        

= −−+        

   

   

   

. (106)

For, kB , the free-surface deformation,

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

BB B B B

abcd

BB B B B

abcd

BB B B B

abcd

kk k k k

kk k k k

kk k k k

= −−+        

= −−+        

= −−+        

. (107)

For, kC, the depth-multiplied term,

11 1 1 1

22 2 2 2

33 3 3 3

 

 

 

CC C C C

abcd

CC C C C

abcd

CC C C C

abcd

kk k k k

kk k k k

kk k k k

= −−+        

= −−+        

= −−+        

. (108)

Computation of X and Y Derivatives

( ) ( )

( ) ( )

111 1

222 2 333 3

222 2 333 3

cos sin

sin cos

x AAB C

y AAB C AAB C

z AAB C AAB C

Up j j j zj

x

Up j j j zj j j j zj

x

Up j j j zj j j j zj

x

∂

= − + +⋅

∂

∂

= − + + ⋅ δ− − + + ⋅ δ

∂

∂ = − + − ⋅ δ+ − + − ⋅ δ

∂



 

 

. (109)

( ) ( )

( ) ( )

111 1

222 2 333 3

222 2 333 3

cos sin

sin cos

x AAB C

y AAB C AAB C

z AAB C AAB C

Up k k k zk

y

Up k k k zk k k k zk

y

Up k k k zk k k k zk

y

∂

= − + +⋅

∂

∂

= − + + ⋅ δ− − + + ⋅ δ

∂

∂ = − + − ⋅ δ+ − + − ⋅ δ

∂



 

 

. (110)

Rotation of the Displacement-Gradient Matrix Back to Cartesian Coordinate System

2 2

2 2

2 2

2

sin ( )sin cos cos

sin ( )sin cos cos

sin ( )sin cos cos

sin ( )sin cos

x x y y

x x y y

y y x x

y y x

u Up Up Up Up

xx y x y

u Up Up Up Up

yy x y x

v Up Up Up Up

xx x y y

v Up Up Up Up

yy y x

φ φφ φ

φ φφ φ

φ

∂ ∂ ∂ ∂ ∂

= −+ +

∂∂ ∂ ∂ ∂

∂ ∂ ∂ ∂ ∂

= +− −

∂∂ ∂ ∂ ∂

∂ ∂ ∂ ∂ ∂

= +− −

∂∂ ∂ ∂ ∂

∂ ∂ ∂

φ φ

∂ ∂

= ++ +

∂∂ ∂ ∂

φ

φ φ φ 2 cos x

∂x

φ

. (111)
54 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Strain

The three strain components that describe the horizontal strain tensor are the areal strain 

εa, and γ1 and γ2 shear strains: 

ε ε ε

γ ε ε

γ ε

a xx yy

xx yy

xy

= + ( ) = − ( ) = ⋅

0 001

0 001

0 001 2

1

2

.

.

.

, (112)

where 

ε ε ε ε xx yy xy yx

u

x

v

y

u

y

v

x = ∂

∂ = ∂

∂ = = ∂

∂

+

∂

∂





 



 , and 1

2 . (113)

The components of the displacement-gradient tensor are given in equation 111. The scale 

0.001 in equation 112 takes into account that, in Okada (1992), coordinates are in kilometers 

and displacements in meters. 

General Parameters

The parameters and variables listed below are the same for fA, fB and fC. They are listed 

in the order in which they would be coded: 

Okada (1992, table 6)

α λ μ

λ μ ν

α

α

α

δ δ δ δ

η

= +

+ = − = −

= + = −

= −

2

1

2

1

1

1

1 ,

p y cos s d q in y d sin cos

p W

R

p p

2 2 2 2 = + +

= + = −

= − = +

ξ η

η δ δ η δ δ

q

y q d q

d c z c d z

 

 

cos sin , sin cos

, 

(114)

where x and η are dummy integration variables used the first time in equations 68; their lower 

and limits are defined in equations 91−94. For example, f f fault c U y z A A = ( ) p , , ξ η, ,,, μ λ , , δ , .

Okada (1992, equation 14; subscripts from tables 6–8) 

X

R R

Y

R

X R

R R

X R R

R R 11

11

32 3 2 53

2 2

5 3

1

1

2 8 9 3 = ( ) +

=

= +

+ = + +

ξ +

ξ

ξ

ξ ξ

ξ , ( ) , ( )

R

h q z

Y R

R R

Y R R

R R

Z

( ) +

= −

= +

+ = + +

η +

δ

η

η

η η

η ,

cos

( ) , ( ) 32 3 2 53

2 2

5 3

3

2 8 9 3

2 3 32 53 5 53

0 11

2

32 0 32

2

53

3 = − = −

= − = −

sin , sin

,

δ δ

ξ ξ

R

hY Z

R

hY

Y Y Y Z Z Z

(115)
Internal Deformation and Strain for Rectangular Dikes and Faults 55

Okada (1992, footnote, table 6)

F

c q

c q

c F c

qR

0

1

2

1 0 2

1

0

0

0

=

= ≤

= >

= +





 



 −

,

,

tan ,

 and

Θ ξη

(116)

where 

 q is a vector; 

 F0 is a vector of zeros the same size of q ; 

 c1 and c2 are logical vectors (a logical vector is equal to 1 if a logical 

condition is satisfied and is equal to 0 otherwise).

Singularity Condition iii (Okada,1992, p. 1034)

F

c R

c R

R c F c

R

R c R c R

0

1

2

1 0 2

1 2

0

0

0

1

=

= + ≤

= + >

= +

+

= − ( ) − + + (

ξ

ξ

ξ

ξ ξ

ξ

ξ ln ln ln );

, (117)

where F0 is a dummy vector of zeros the same size as R + ξ, and c1 and c2 are dummy logical 

vectors.

Singularity Condition iv (Okada,1992, p. 1034)

F

c R

c R

c F c

R

R c R c

0

1

2

1 0 2

1 2

0

0

0

1

=

= + ≤

= + >

= +

+

= − ( ) − +

η

η

η

η

η

η

R

ln ln ln ; ( ) R + η

(118)

where F0 is a dummy vector of zeros the same size as R + η , and c1 and c2 are dummy logical 

vectors.
56 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

I Parameters (Okada, 1992, footnote, table 6)

Singularity Condition ii (Okada,1992, p. 1034)

F

c

c

X q

0

1

2

2 2

0

0

0

=

= ≤

= >

= +

ξ

ξ

ξ

, (119)

where F0 is a dummy vector of zeros the same size of R + η , and c1 and c2 are dummy logical 

vectors. 

I y

R d

R R d

I c F c

3 2

4 1 0 2

1 1 = +

− +  ( ) − +  



= +

cos cos

ln sin ln( )

sin

δ δ η δ 

 

δ

δ

ξ

δ

η δ δ

cos cos ξ δ

tan ( cos ) ( )sin

R d ( ) cos

X q X R X

+ R X

+

+ + +

+





  −



2

2

1



 



 





= − + −

= + ( ) +

I

R

I

I R I

d

d

1 4

2 3

ξ δ δ

δ





cos sin

log sin

 

 

; (120)

if cos δ = 0 , then

I

R d

yq

R d

R

I c F c y

R

3 2

4 1 0 2

0 5

0 5

= ++

( ) +

− + ( )

= + ⋅





















. ln

.

η η

ξ









( ) + 

d

2

. (121)

J and K Parameters (Okada, 1992, footnote, table 7)

D

R R d

J y

R d

D

J d y

R d

D

11

2 11

5

2

11

1 =

( ) +

= +

= − +

+























 



ξ , (122)
Internal Deformation and Strain for Rectangular Dikes and Faults 57

11 11

1

11 11

3

1 2

3

3 5

6

sin cos

cos

sin cos

sin

cos

D Y

K

qY yD K

K J J

K J J

− δ = x δ

− = δ

− δ = δ

− δ = δ



, (123)

15 6

4 11 2 3

2 3

4 11 1

cos sin

cos sin

1 sin

cos sin

JJ J

J YJ J

K K

R

KY K

= δ− δ

= −x − δ + δ

=+ δ

=x δ− δ

; (124)

if cos δ = 0 ,

( )

( )

11

1

2

11

3

2

11

3 2

2

11

6 2

.* 1 sin

0.5

0.5

qD K

R d

D

K

R d

q D J

R d

D

J y

R d

x = +

x − = δ

+

− = −x

+

x − = −

+







 

. (125)

Parameters E–Q (Okada, 1992, footnote, table 8)

E

R

yq

R

F d

R

Y

y

d

G X qX

H qX

= −

= +

= −

= +

sin

sin

sin

δ

ξ δ

δ

ξ









3

3

2

32

11 32

32

2 

qY

P

R

qY

Q cd

R

zY Z Z

32

3 32

5 32 32 0 3

sin

cos sin

sin

δ

δ δ

δ

= +

= − ( ) + +  

(126)
58 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Table 18. Strike-slip fault (internal deformation and strain; function fA )

[Okada (1992), type strike—for displacement, see table 6; for X-derivative, see table 7, and 

for Y-derivative, see table 8]

Displacement 

1

A

u ( ) 11 0.5 Θ + αxqY

u

A

2 0.5 q α R

u

A

3 ( ) ( ) 2

11 0.5 1 ln 0.5 −α +η − α R qY

X-derivative 

j

A

1 − − 0 5( ) 1 0 − 5 11

2

32 . . α α qY qY ξ

2

A

j 3 − αx 0.5 q R

3

A

j ( ) 2

11 32 0.5 1 0.5 − α x + αx Y qY

Y-derivative 

1

A k ( ) 11 11 0.5 1 sin 0.5 0.5 − α x δ + + αx Y dX F 

2

A k 0.5 E α

3

A k 0.5 1 cos sin 0.5 ( ) −α δ + δ − α ( R qY11 ) qF

Table 19. Dip-slip fault (internal deformation and strain; function fA )

[Okada (1992), type dip—for displacement, see table 6; for X-derivative, see table 7, and 

for Y-derivative, see table 8] 

Displacement

1

A

u 0.5αq R

2

A

u 11 0.5 0.5 qX Θ + αη

3

A

u 0 5 1 0 5 2

11 . l ( ) − α ξ n . ( ) R q + − α X

X-derivative

j

A

1 −0 5 3 . αξ q R

j

A

2 − − 05 05 11

3 . . qY αηq R

j

A

3 0 5 1 0 5 2 3 . . ( ) − α α R q + R

Y-derivative

k A

1

0 5. αE

k A

2 0 5 1 0 5 0 5 11 11 . . ( ) − α ξ dX + + Y G sin . δ αη

k A

3 0 5 1 0 5 11 . . ( ) − α α yX − qG

Function fA (Infinite Medium)
Internal Deformation and Strain for Rectangular Dikes and Faults 59

Table 21. Scaling used to compute the 

dimensional displacement (function fA ).

[Okada (1992), i=1: strike-slip; i=2: dip-slip; i=3: tensile 

e—for displacement, see table 6; for X-derivative, see 

table 7, and for Y-derivative, see table 8]

Displacement (Okada, 1992; Table 6)

u

A

1 0 5 1 . U u i

A ( ) π

u

A

2 0 5 2 . U u i

A ( ) π

u

A

3 0 5 3 . U u i

A ( ) π

X-Derivative (Okada, 1992; Table 7)

j

A

1 0 5 1 . U j i

A ( ) π

j

A

2 0 5 2 . U j i

A ( ) π

j

A

3 0 5 3 . U j i

A ( ) π

Y-Derivative (Okada, 1992; Table 8)

k A

1 0 5 1 . U k i

A ( ) π

k A

2 0 5 2 . U k i

A ( ) π

k A

3 0 5 3 . U k i

A ( ) π

Table 20. Tensile crack (internal deformation and strain; function fA )

[Okada (1992), type tensile—for displacement, see table 6; for X-derivative, see table 7, and 

for Y-derivative, see table 8]

Displacement

u

A

1 − − 0 5( ) 1 0 ( ) + − 5 211 . l α η n . R qα Y

u

A

2 − − 0 5( ) 1 0 ( ) + − 5 211 . l α ξ n . R qα X

u

A

3 ( ) 11 11 0.5 0.5 Θ− α η +x qX Y

X-derivative

j

A

1 − − 0 5( ) 1 0 + 5 11

2

32 . . α ξY q αξ Y

j

A

2 − − 0 5( ) 1 0 + 5 2 3 . . α α R q R

j

A

3 − − 0 5( ) 1 0 − 5 11

3

32 . . α α qY q Y

Y-derivative

k A

1 − − 0 5( ) 1 0 ( ) + − 5 11 . c α δ os R qY q sin . δ α F

k A

2 ( ) 11 0.5 1 0.5 − −α − α yX qG 

k A

3 0 5 1 0 5 11 11 . s ( ) − α ξ ( ) + in δ α + . 

dX Y qH
60 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Table 22. Strike-slip fault (internal deformation and strain; function fB ).

[Okada (1992), type strike—for displacement, see table 6; for X-derivative, see table 7, 

and for Y-derivative, see table 8] 

Displacement

u

B

1 − − ξ α qY I − δ 11 Θ 1 1 sin

u

B

2 − + q R α δ y R( ) + d 1   sin

u

B

3 q Y I 2

11 − α1 2 sin δ

X-derivative

j

B

1 ξ α δ 232 1 1 qY J − sin

j

B

2 ξ α q R J δ 3 − 1 2 sin

j

B

3 − − ξ α q Y J δ 232 1 3 sin

Y-derivative 

k B

1 − − ξ α F dX Y + + ( ) ξ δ J  11 1 11 4 sin

k B

2 − +α + δ E RJ 1 5 ( ) 1 sin

k B

3 qF qY J −α − δ 1 11 6 ( )sin

Table 23. Dip-slip fault (internal deformation and strain; function fB ). 

[Okada (1992), type dip—for displacement, see table 6; for X-derivative, see table 7, and for 

Y-derivative, see table 8] 

Displacement

1

B

u ( ) ( ) 3 − + −α α δ δ q R 1 / sin cos I

2

B

u −η − Θ − − α α x + δ δ qX11 ( ) ( ) 1 ( ) ( ) R d sin cos 

3

B

u ( ) ( ) 2

11 4 q X + −α α δ δ 1 sin cos I

X-derivative

1

B

j 3

1 4 x +α δ δ qR J sin cos

2

B

j 3

11 1 5 η + +α δ δ q R qY J sin cos

3

B

j 2 3

1 6 − +α δ δ qR J sin cos

Y-derivative

1

B k 1 1 − +α δ δ E J sin cos 

2

B k 11 1 2 sin sin cos −η − x α + α δ δ GY J

3

B k 1 3 qG J sin cos +α δ δ

Function fB (Free Surface Deformation)
Internal Deformation and Strain for Rectangular Dikes and Faults 61

Table 25. Scaling (internal deformation and 

strain; function fB )

[Okada (1992), i=1: strike-slip; i=2: dip-slip; i=3: tensile 

e—for displacement, see table 6; for X-derivative, see 

table 7, and for Y-derivative, see table 8]

Displacement

u

B

1 0 5 1 . U u i

B ( ) π

u

B

2 0 5 2 . U u i

B ( ) π

u

B

3 0 5 3 . U u i

B ( ) π

X-derivative

j

B

1 0 5 1 . U j i

B ( ) π

j

B

2 0 5 2 . U j i

B ( ) π

j

B

3 0 5 3 . U j i

B ( ) π

Y-derivative

k B

1 0 5 1 . U k i

B ( ) π

k B

2 0 5 2 . U k i

B ( ) π

k B

3 0 5 3 . U k i

B ( ) π

Table 24. Tensile crack (internal deformation and strain; function fB ). 

[Okada (1992), type tensile—for displacement, see table 6; for X-derivative, see table 7, 

and for Y-derivative, see table 8] 

Displacement

u

B

1 q Y I 2

11 1 3

2 − α sin δ

u

B

2 q X R dtilde

2

11 1

2 + + α ξ( ) ( ) sin δ

u

B

3 q Xη ξY I α δ 11 11 1 4

2 ( ) + − − Θ sin

X-derivative

j

B

1 − − ξ α q Y J δ 232 1 4

2 sin

j

B

2 − − q R J 2 31 5

2 α δ sin

j

B

3 q Y J 332 1 6

2 − α sin δ

Y-derivative

k B

1 qF J − α δ 1 1

2 sin

k B

2 qG J − α δ 1 2

2 sin

k B

3 − − qH J α δ 1 3

2 sin
62 Modeling Crustal Deformation near Active Faults and Volcanic Center—A Catalog of Deformation Models

Table 26. Strike-slip fault (internal deformation and strain; function fC ). 

[Okada (1992), type strike—for displacement, see table 6; for X-derivative, see table 7, and for Y-derivative, see table 8]

Displacement

uC

1 1( ) − α ξ 11 − 32 Y q cos δ αξ Z

uC

2 1 2 11

3 ( ) − α δ ( ) cos s R q + Y c in δ α−  q R

uC

3 1 11

3

11

2 ( ) − α δ − − α η ( ) + 32 qY cos c R  zY ξ Z

X-derivative

j

C

1 1( ) − α δ Y q 0 0 cos − α Z

j

C

2 − − ( ) 1 2 ( ) + + 3 3

32

5 α ξ cos s δ δ R qY c in α ξ q R

j

C

3 − − ( ) 1 3 32 + − ( ) − − 5 α ξ 32 32 0 qY cos δ αξ ηc R  zY Z Z

Y-derivative

kC

1 − − ( ) 1 α ξP Q cos δ α− ξ

kC

2 2 1 3 3

0

3 3 3 ( ) − α δ ( ) dR −Y y sin sin δ δ − −   cos s R c α δ ( ) + d R  in − − η R cy q R5 

 

kC

3 − − ( ) 1 3 + − ( ) + + ( ) + − 3 3

0

3 5 α δ 0 q R y R   Y c δ α d R δ cd q R Y    cos sin  cos / ( ) cos δ δ +  

 qZ0 sin

Function fC (Depth Multiplied Term)
Internal Deformation and Strain for Rectangular Dikes and Faults 63

Table 27. Dip-slip fault (internal deformation and strain; function fC ). 

[Okada (1992), type dip—for displacement, see table 6; for X-derivative, see table 7, and for Y-derivative, see 

table 8]

Displacement

uC

1 1 11

3 ( ) − α δ cos s R q − − Y c in δ α q R

uC

2 1( ) − α α 11 − η 32 yX c qX

uC

3 − − − − ( ) 

dX11 11 Y c X q 11 X2

32 ξ δ sin α

X-derivative

j

C

1 − − ( ) 1 3 + + 3

32

5 α ξcos s δ ξ R qY c in δ α ξ q R

j

C

2 − − ( ) 1 3 +3 5 α α y R  c q η R

j

C

3 

d R Y c q R R 3

0

2 2 3 − + sin δ α ( ) 1 3 −

Y-derivative

kC

1 − − ( ) + − ( ) + −  

 1 3 3

0

2 3 5 α η R Y sin s δ α δ c d   in R cy q R

kC

2 1 2 11

2 ( ) − ( ) − 32 − + ( ) 32 53 −  

 α α X y  X c d q  cos δ η X y qX

kC

3 ξ δ P y sin s + + dX α δ c y ( ) + q X in − yq X  

      32 32

2

53 2
64 Modeling Crustal Deformation near Active Faults and Volcanic Center—A Catalog of Deformation Models

Table 29. Scaling (internal deformation and strain; function fC ).

[Okada (1992), i=1: strike-slip; i=2: dip-slip; i=3 tensile e—for displacement, 

see table 6; for X-derivative, see table 7, and for Y-derivative, see table 8]

Displacement

uC

1 0 5 1 . U u i

C ( ) π

uC

2 0 5 2 . U u i

C ( ) π

uC

3 0 5 3 . U u i

C ( ) π

X-derivative

j

C

1 0 5 1 . U j i

C ( ) π

j

C

2 0 5 2 . U j i

C ( ) π

j

C

3 0 5 3 . U j i

C ( ) π

Y-derivative

kC

1 0 5 1 . U k i

C ( ) π

kC

2 0 5 2 . U k i

C ( ) π

kC

3 0 5 3 . U k i

C ( ) π

Table 28. Tensile crack (internal deformation and strain; function fC ).

[Okada (1992), type tensile—for displacement, see table 6; for X-derivative, see table 7, and for Y-derivative, see table 8] 

Displacement

uC

1 − − ( ) 1 ( ) + 11 − − ( ) 11

2 α δ 32 sin c R qY z os δ α Y q Z

uC

2 1 2 11 11 11

2 ( ) − α ξ + − ( ) − 32 Y d sin δ α 

X c X q X

uC

3 1( ) − α ξ ( ) 11 11 + + + ( ) 32 32 yX  Y q cos δ α c Xη ξZ

X-derivative

j

C

1 1 3 2 3

32

5 ( ) − α ξ + + ( ) − − 32 0 sin c δ ξ R qY c os δ αη η R Z Z

j

C

2 1 2 1 3 0

3 2 2 3 ( ) − α δ Y d sin − +  R c α( ) − q R R

j

C

3 − − ( ) 1 3 ( ) − − − ( ) 3

0

5 α δ α η 0 y R   Y c cos q R qZ

Y-derivative

kC

1 1 3 3

0

3 5 ( ) − α δ ( ) + + + ( ) − 0 q R Y z sin cos δ α cos s δ δ R cd q R qZ in 

kC

2 − − ( ) − + ( ) + −  

 1 2 2 32 32

2 α ξ 53 P y sin s δ α dX c d   q X in δ yq X

kC

3 − − ( )( ) − + + + ( ) −  

 1 2 11 + 2 α ξ 32 32 53 P X cos c δ α y X  c d q X os δ ηy qX αξ   Q
Internal Deformation and Strain for Rectangular Dikes and Faults 65

Verification and Examples

Strike-Slip Fault (see fig. 33 and table 30)

−10 −5 0 5 10 −10

−8

−6

−4

−2

0

2

4

6

8

10

−10

−8

−6

−4

−2

0

2

4

6

8

10 Okada92

Y, in kilometers

−10 −5 0 5 10

Coulomb

−10 −5 0 5 10 −0.018

−0.016

−0.014

−0.012

−0.010

−0.008

−0.006

−0.004 U (east), in meters

−10 −5 0 5 10

−0.020

−0.015

−0.010

−0.005

0.000

0.005

0.010

0.015

0.020

V (north), in meters

Coulomb

Okada92

−10 −5 0 5 10

−0.010

−0.005

0.000

0.005

0.010

W (up), in meters 

−10 −5 0 5 10

−4

−3

−2

−1

0

1

2

3

4

εa, × 10-6

−10 −5 0 5 10

−8

−6

−4

−2

0

2

4

6

8

X, in kilometers

γ1, × 10-6

Coulomb

FEM

Numerical

Okada92

−10 −5 0 5 10 −1.0

−0.5

0

0.5

1.0

γ2, × 10-5

EXPLANATION

EXPLANATION

Figure 33. Vector (top) and profile plots showing displacement field (center) and strain (bottom) caused by a strike-slip fault 

(thick red line). The profile was taken at a depth of 1 kilometer along the magenta line (top). The difference between the Okada92 

solution for γ2 and that from COULOMB 3.3 is due to an error in COULOMB 3.3 (R. Stein, U.S. Geological Survey, oral commun.). 

Okada92, Matlab script of Okada’s (1992) expressions; Coulomb 3.3, solution from COULOMB 3.3; FEM, numerical solution from a 

finite element method numerical model; numerical, strain from finite-difference derivatives (appendix 1). Fault parameters are xi

 = 0 

meters, yi

 = 0 meters, xf = 3,000 meters, yf = 0 meters, δ =90°, μ = 1 Gigapascals, ν = 0.25, U = 1 meter, zt =2,000 meters, and zb = 

4,000 meters. See table 30 for numerical example.
66 Modeling Crustal Deformation near Active Faults and Volcanic Center—A Catalog of Deformation Models

Table 30. Internal deformation and strain around a strike-slip fault.

[See also figure 33. Parameters are xi= 0 meters, yi = 0 meters, xf = 3,000 meters; yf= 0 meters; zt = 2,000 meters; zb = 4,000 meters; x, y and z, local Cartesian 

coordinates; u u x y , and uz, internal deformation; ε γ a , 1 and γ 2 , strain; μstrain , a strain expressed in terms of parts per million]

x

(kilometers)

y

(kilometers)

z

(kilometers)

ux

(meters)

uy

(meters)

uz

(meters)

εa

(μstrain)

γ1

(μstrain)

γ2

(μstrain)

 -5.0000 -5.0000 -1.0000 0.0103 0.0092 -0.0014 1.6162 1.1973 2.9945

-1.0000 0.0000 -1.0000 -0.0000 0.0186 -0.0000 0.0000 0.0000 -27.5646

3.0000 -1.0000 -1.0000 0.0283 -0.0259 0.0230 2.1363 -5.9239 -14.7870

 7.0000 5.0000 -1.0000 -0.0116 -0.0113 -0.0023 1.9409 0.3598 3.8515
Internal Deformation and Strain for Rectangular Dikes and Faults 67

Dip-Slip Fault (see fig. 34 and table 31)

men12-7629_fig34

−4 −2 0 2 4 6

−5

−4

−3

−2

−1

0

1

2

3

4

5

−5

−4

−3

−2

−1

0

1

2

3

4

5

Okada92

Y, in kilometers

−4 −2 0 2 4 6

Coulomb

−5 0 5 10

−0.010

−0.005

0.000

0.005

0.010

U (east), in meters

−5 0 5 10

−0.030

−0.025

−0.020

−0.015

−0.010

−0.005

0

V (north), in meters

−5 0 5 10

−0.030

−0.025

−0.020

−0.015

−0.010

−0.005

0.000

W (up), in meters

−10 −5 0 5 10

0

1

2

3

4

5

εa, × 10-6

−10 −5 0 5 10

−20

−15

−10

−5

0

5

X, in kilometers

γ1, × 10-6

−10 −5 0 5 10

−1.0

−0.5

0

0.5

1.0

γ2, × 10-5

Coulomb

FEM

Numerical

Okada92

Coulomb

Okada92

EXPLANATION

EXPLANATION

Figure 34. Vector (top) and profile plots of the displacement field (center) and strain (bottom) caused by a dip-slip fault (thick 

red line). The profile was taken at a depth of 1 kilometer along the magenta line (top). The difference between the Okada92 

solution for γ2 and that from COULOMB 3.3 is due to a error in COULOMB 3.3 (R. Stein, U.S. Geological Survey, oral commun.). 

Okada92, Matlab script of Okada’s (1992) expressions; Coulomb 3.3, solution from COULOMB 3.3; FEM, numerical solution from 

FEM model; numerical, strain from finite-difference derivatives (appendix 1). Fault parameters: xi

 = 0 meters, yi = 0 meters,

xf = 3,000 meters; yf = 0 meters; δ =90°; μ = 1 GPa; ν = 0.25; U = 1 meters; zt =2,000 meters; zb = 4,000 meters.
68 Modeling Crustal Deformation near Active Faults and Volcanic Center—A Catalog of Deformation Models

Table 31. Internal deformation and strain around a dip-slip fault. 

[See figure 34. Parameters are xi = 0 meters, yi = 0 meters, xf = 3,000 meters; yf = 0 meters; δ =90°; μ= 1 GPa; ν = 0.25; U = 1 meters; zt =2,000 meters; zb = 

4,000 meters; x, y and z, local Cartesian coordinates; u u x y , and uz, internal deformation; ε γ a , 1 and γ 2 , strain; μstrain , a strain expressed in terms of parts per 

million]

x

(kilometers)

y

 (kilometers)

z

(kilometers )

ux

(meters)

uy

(meters)

uz

(meters)

εa

(μstrain)

γ1

(μstrain)

γ2

(μstrain)

 -5.0000 -5.0000 -1.0000 0.0041 -0.0038 0.0024 -0.7843 -1.0461 -1.9073

-1.0000 0.0000 -1.0000 0.0000 -0.0195 0.0000 -0.0000 0.0000 5.9463

3.0000 -1.0000 -1.0000 0.0177 -0.0486 0.0653 5.7981 0.4789 26.5485

 7.0000 5.0000 -1.0000 0.0051 -0.0056 -0.0036 1.0556 0.8267 2.9135
Internal Deformation and Strain for Rectangular Dikes and Faults 69

Tensile Crack (see fig. 35 and table 32)

men12-7629_fig35

−4 −2 0 2 4 6

−5

−4

−3

−2

−1

0

1

2

3

4

5

Okada92

Y, in kilometers

−4 −2 0 2 4 6

−5

−4

−3

−2

−1

0

1

2

3

4

5

Coulomb

−5 0 5 10

−0.010

−0.005

0.000

0.005

0.010

U (east), in meters

−5 0 5 10

0

0.01

0.02

0.03

0.04

0.05

0.06

V (north), in meters

−5 0 5 10

0

0.005

0.010

0.015

0.020

0.025

0.030

W (up), in meters

−5 0 5 10

−4

−3

−2

−1

0

1

× 10−6

εa

−5 0 5 10

−5

0

5

10

15

20

× 10−6

X, in kilometers

γ1

Coulomb

FEM

Numerical

Okada92

Coulomb

Okada92

−5 0 5 10

−1.5

−1.0

−0.5

0.0

0.5

1.0

1.5

× 10−5

γ

EXPLANATION 2

EXPLANATION

Figure 35. Vector (top) and profile plots of the displacement field (center) and strain (bottom) caused by a tensile crack

(thick red line). The profile was taken at a depth of 1 kilometer along the magenta line (top). The difference between the 

Okada92 solution for γ2 and that from COULOMB 3.3 is due to a error in COULOMB 3.3 (R. Stein, U.S. Geological Survey, oral 

commun.). Okada92, Matlab script of Okada’s (1992) expressions; Coulomb 3.3, solution from COULOMB 3.3; FEM, numerical 

solution from FEM model; numerical, strain from finite-difference derivatives (appendix 1). Fault parameters: xi

 = 0 meters, 

yi

 = 0 meters, xf = 3,000 meters; yf = 0 meters; δ =90°; μ = 1 GPa; ν = 0.25; U = 1 m; zt =2,000 meters; zb = 4,000 meters.
70 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Table 32. Internal deformation and strain around a tensile crack. 

[See figure 35. Parameters are xi = 0 meters, yi = 0 meters, xf = 3,000 meters; yf = 0 meters; δ=90°; μ=1 Gigapascal; 

ν= 0.25, U=1 meter, zt=2,000 meters; zb = 4,000 meters. x, y and z, local Cartesian coordinates; u u x y , and uz , internal 

deformation; ε γ a , 1 and γ 2 , strain; μstrain , a strain expressed in terms of parts per million]

x 

(kilometers)

y

 (kilometers)

z

(kilometers)

ux

(meters)

uy

(meters)

uz

(meters)

εa

(μstrain)

γ1

(μstrain)

γ2

(μstrain)

 -5.0000 -5.0000 -1.0000 -0.0041 -0.0079 0.0041 -0.0138 -2.6669 -1.9129

-1.0000 0.0000 -1.0000 0.0152 -0.0000 -0.0058 14.8157 -7.2473 -0.0000

3.0000 -1.0000 -1.0000 -0.0013 -0.0296 0.0141 27.4109 -32.0105 2.6374

 7.0000 5.0000 -1.0000 0.0056 0.0114 0.0055 -0.2245 -2.9155 -3.4452
Discretized Faults 71

Discretized Faults

Inverse modeling of surface-displacement measurements is required to estimate spatial slip variations along the faults 

that ruptured in an earthquake. Appropriate analytical solutions can relate surface displacements to offsets across many small 

rectangular dislocations (referred to as segments, subsegments, and tiles4 in the model; see fig. 36) that may or may not slip in a 

halfspace that is homogeneous, isotropic, flat, and elastic.

4

Tiles may also be called patches in the literature.

Figure 36. Representation of the fault geometry by segments, subsegments, and tiles. The fault above has 5 segments, 13 

subsegments and 78 tiles. The fault is identified by its initial x y z i i t ( ) , , and final x y z f f t ( ) , , coordinates; zt, depth of the 

top edge of the fault.

5

4 6 8 10 12 14 16 18 20 22 24

10

15

20

0

1

2

3

4

5

6

7

8

9

10

X, in kilometers

Depth, in kilometers

segment

tile

(xi, yi, zt)

(xf

, yf, zt)

Y, in kilometers

subsegment

subsegment

subsegment

This segment is divided 

in 3 subsegments and 

3 X 6 = 18 tiles 

(patches)
72 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

The chart below (fig. 37) describes the process for coding a discretized fault.

men12-7629_fig37

SEGMENTS

divide a fault into segments, 

and each segment into several subsegment

FAULT

define fault geometry; use matrix of tiles from PATCHES and OKADA Green’s function to 

compute displacement, tilt, and strain

PATCHES

approximate a fault by a matrix of rectangular tiles

OKADA

compute the Green's function

Figure 37. The Matlab script for coding a discretized fault includes four major steps: (1) the 

function SEGMENTS divides a fault into several rectangular dislocations (segments and 

subsegments); (2) each segment/subsegment is divided into several tiles by the function 

PATCHES, approximating the fault with a matrix of tiles; (3) the function FAULT uses the function 

OKADA to compute the displacement of each tile, and then (4) add the displacements of all of 

the tiles to compute the crustal deformation caused by the slip along the fault. 
Discretized Faults 73

Fault Segments and Subsegments

Estimating the location, size, dip, orientation, and slip across multiple rectangular 

dislocations that may or may not slip is a numerically challenging task given the number of 

data points and the number of parameters required to estimate these dislocations. Information 

from the mapped fault trace can be used to constrain several parameters of the fault geometry, 

including the number of segments that could reasonably represent the fault trace at the surface: 

L x x y y n

i

n

f

n

i

n

f

n

N

= − ( ) + − ( ) =

∑ )

2 2

1

(127)

and

N = 



 



 length of fault trace

width of segment , (128)

where 

L is the length of the fault trace, 

( , x y ) n

i

n

i and 

( , x y ) n

f

n

f are the initial and final points of the segment,

n is index of the nth segment, 

N is the number of segments, and 

  is the ceiling function (fig. 36). The (x, y) coordinates of points of the grid 

along the upper edge of each subsegment (figs. 36 and 38) are defined by 

( )

( )

0:1:

f i

i n n n

f

f fi n n nn f i

n n

x x xx x S

S

xx x

yy y y y

x x

− = +

− =+ − −

, (129)

where S is the total number of subsegments, and each subsegment’s initial and final point

( ) ( )

( ) ( )

 1: , 1: 

2: 1 , 2: 1

i i

s s

f f

s s

x xx S y yy S

x xx S y yy S

= =

=+=+

. (130)
74 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Tiles

Each subsegment (figs. 36 and 38) can be divided into several tiles following the steps described below. First, define three 

parameters:

1. the subsegment at the initial (lower left) corner

 x x z z

y y z z

i

s

i

s

b

s

t

i

s

i

s

b

s

t

0

0

= + ( ) −

= − ( ) −

cot cos

cot sin

δ φ

δ φ

, (131)

men12-7629_fig38

zt

zb

(xi

n, yin)

(xf

n, yfn)

X 

Y 

Z (depth )

tile p

(xi

p, yip, zbp) (xf

p, yfp, zbp)

(xi

p, yip, zt

p) (xf

p, yfp, zt

p)

(xi

s, yis, zt

s) (xf

s, yfs, zt

s)

(xi

s, yis, zbs) (xf

s, yfs, zbs)

subsegment s

segment n

Figure 38. Segment, subsegment, and tile geometry. x y

i i ( ) , , horizontal coordinates for the left corner of segments, 

subsegments, or tiles; x y f f ( ) , , horizontal coordinates for the right corner segments, subsegments, or tiles; zt , depth of 

the top edge of segments, subsegments or tiles; zb , depth of the bottom edge of segments, subsegments, or tiles; n, s or p, 

subindex for segments, subsegments, or tiles.
Discretized Faults 75

2. the subsegment at the final (lower right) corner

 x x z z

y y z z

f

s

f

s

b

s

t

f

s

f

s

b

s

t

0

0

= + ( ) −

= − ( ) −

cot cos

cot sin

δ φ

δ φ

, and (132)

3. the number of tiles P for each subsegment (the symbol   is the ceiling function)

P = 



 



 height of subsegment

height of tile . (133)

Then, compute the coordinates of the grid of tiles at the initial (left) and final (right) side of 

each segment by using the formulas below:

1. grid along the depth

zz z P s z z P t

s

b

s

t = + ( ) 0 1: : ( ) − , (134)

2. grid along the left side of the segment

xx x x x zz z z z

yy y y y zz z

i

i

s

i i

s

b

s

t

s

b

i

i

s

i i

s

b

= + ( ) − ( ) − ( ) −

= + ( ) − ( −

0 0

0 0 ) ( ) z z s −t

s

b

 and (135)

3. grid along the right side of the segment

 ( )( ) ( )

( )( ) ( )

0 0

0 0

f f f b tb

f s s ss

f f f b tb

f s s ss

xx x x x zz z z z

yy y y y zz z z z

=+ − − −

=+ − − −

. (136)

Finally, we define the tile at the top left corner (initial point), the tile at the top right corner 

(final point), and the tiles at the top and bottom as

( ) ( )

( ) ( )

( ) ( )

1: 1:

1: 1:

1: 2: 1

i i

pi pi

f f

pf pf

t b

p p

x xx P y yy P

x xx P y yy P

z zz P z zz P

= =

= =

= = +

. (137)

Verification of the Discretized Fault Model

We approximate a fault by a matrix of rectangular tiles. The model computes the 

deformation of a single tile by using Green’s function for a rectangular dislocation (Okada, 

1985, 1992) and then adds the contributions to calculate the total deformation (see tables 33, 34

and 35 for numerical examples). We compared our Matlab coding of the analytical expressions 

for discretized strike-slip faults (fig. 39), dip-slip faults (fig. 40), and tensile cracks (fig. 41) 

faults against solutions from COULOMB 3.3.
76 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Figure 39. Fault geometry and vector plot (top row) and profile plots (bottom row) of the surface displacement caused by a rightlateral strike-slip fault (thick red line). The profile is along the magenta line (top row, right). Top row—Okada85, Matlab script of Okada’s 

(1985) expressions; COULOMB, solution from COULOMB 3.3. Bottom row—“Discretized fault” labels the solution obtained by applying 

the Okada (1985) equations to a fault composed of 144 tiles (top row, left). “Single fault” labels the solution from COULOMB 3.3 for a 

fault divided into five segments, each of which exhibited uniform slip. See table 33 for numerical example.

Strike-Slip Fault

Table 33. Parameters for strike-slip fault segments. 

[Parameters given here are for right-lateral strike-slip fault segments. See figure 39. xi, yi, xf, and yf, coordinates of the initial and final positions along the top 

edges of the fault segments; U1, final right-lateral dislocation of segment; zt and zb, coordinates of the top and bottom of the segment]

Segment xi

 (kilometers)

yi

(kilometers)

xf

(kilometers)

yf

(kilometers)

U1

(meters)

Dip angle

(degrees)

zt

(kilometers)

zb

(kilometers)

1 5 6 9 8 1 85 2 10

2 9 8 11 13 0.5 90 2 10

3 11 13 15 16 1.5 75 2 10

4 15 16 18 18 0.5 80 2 10

5 18 18 21 20 0.75 70 2 10

men12-7629fig39

5 10 15 20 25 5

10

15

20

0

2

4

6

8

10

Strike-slip fault discretized into 144 tiles

X, in kilometers

Y, in kilometers Depth, in kilometers

0 5 10 15 20 25

0

5

10

15

20

25

Discretized fault

(Okada85)

Y, in kilometers

0 5 10 15 20 25

0

5

10

15

20

25

Single fault

(Coulomb)

X, in kilometers

0 5 10 15 20 25

−0.08

−0.06

−0.04

−0.02

0.00

0.02

0.04

U (east) , in meters

0 5 10 15 20 25

−0.12

−0.00

−0.08

−0.06

−0.04

−0.02

0.00

X, in kilometers

V (north) , in meters

Single fault 

(Coulomb)

Discretized fault 

(Okada85)

0 5 10 15 20 25

−0.01

0.00

0.01

0.02

0.03

0.04

0.05

0.06

W (up) , in meters

EXPLANATION
Discretized Faults 77

Figure 40. Fault geometry and vector plot (top row) and profile plots (bottom row) of the surface displacement caused by a reverse 

dip-slip fault. The profile is along the magenta line (top row, right). Top row—Okada85, Matlab script of Okada’s (1985) expressions; 

COULOMB, solution from COULOMB 3.3. Bottom row—“Discretized fault” labels the solution obtained by applying the Okada (1985) 

equations to a fault composed of 78 tiles (top row, left). “Single fault” labels the solution from COULOMB 3.3 for a fault divided into five 

segments, each of which exhibited uniform slip. See table 34 for numerical example.

men12-7629_fig40

5 10 15 20 25 5

10

15

20

0

2

4

6

8

10

X, in kilometers

Depth, in kilometers

0 5 10 15 20 25

0

5

10

15

20

25

Y, in kilometers

0 5 10 15 20 25

0

5

10

15

20

25

X, in kilometers

0 5 10 15 20 25

−0.06

−0.04

−0.02

0

0.02

0.04

0.06

U (east), in meters

0 5 10 15 20 25

−0.12

−0.1

−0.08

−0.06

−0.04

−0.02

0

X, in kilometers

V (north), in meters

0 5 10 15 20 25

−0.05

0

0.05

0.1

0.15

0.2

W (up), in meters

Dip-slip fault discretized into 78 tiles

Discretized fault

(Okada85)

Single fault

(Coulomb)

Y, in kilometers

Single fault

(Coulomb)

Discretized fault

(Okada85)

EXPLANATION

Dip-Slip Fault

Table 34. Parameters for dip-slip fault segments.

[Parameters given are for reverse dip-slip fault segments. See figure 40. xi, yi, xf, and yf, coordinates of the initial and final positions along the top edges of the 

fault segments; U2, final reverse dislocation of segment; zt and zb, coordinates of the top and bottom of the segment.

Segment xi

 (kilometers)

yi

(kilometers)

xf

(kilometers)

yf

(kilometers)

U2

(meters)

Dip angle 

(degrees)

zt

(kilometers)

zb

(kilometers)

1 5 6 9 8 1 85 2 10

2 9 8 11 13 0.1 90 2 10

3 11 13 15 16 0.7 75 2 10

4 15 16 18 18 0.5 80 2 10

5 18 18 21 20 1.5 70 2 10
78 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

men12-7629fig41

5 10 15 20 25 5

10

15

20

0

2

4

6

8

10

0 5 10 15 20 25

0

5

10

15

20

25

0 5 10 15 20 25

0

5

10

15

20

25

Tensile (opening) crack discretized into 78 tiles

X, in kilometers

Y, in kilometers Depth, in kilometers

Discretized fault

(Okada85)

Y, in kilometers

Single fault

(Coulomb)

X, in kilometers

0 5 10 15 20 25

0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0 5 10 15 20 25

−0.14

−0.12

−0.10

−0.08

−0.06

−0.04

−0.02

0.02

0 5 10 15 20 25

0.02

0.04

0.06

0.08

0.10

0.12

0

0

U (east), in meters

X, in kilometers

Single fault

V (north), in meters

W (up), in meters

(Coulomb)

Discretized fault

(Okada85)

EXPLANATION

Figure 41. Fault geometry and vector plot (top row) and profile plots (bottom row) of the surface displacement caused by the opening 

of a tensile crack. The profile is along the magenta line (top row, right). Top row—Okada85, Matlab script of Okada’s (1985) expressions; 

COULOMB, solution from COULOMB 3.3. Bottom row—“Discretized fault” labels the solution obtained by applying the Okada (1985) 

equations to a fault composed of 78 tiles (top row, left). “Single fault” labels the solution from COULOMB 3.3 for a fault divided into five 

segments, each of which exhibited uniform slip. See table 35 for numerical example.

Table 35. Parameters for tensile-crack fault segments.

[Parameters given are for the tensile-crack shown in figure 41. xi, yi, xf, and yf, coordinates of the initial and final positions along the top edges of the fault 

segments; U3, final opening dislocation of segment; zt and zb, coordinates of the top and bottom of the segment]

Segment xi

 (kilometers)

yi

(kilometers)

xf

(kilometers)

yf

(kilometers)

U3

(meters) §

Dip angle

(degrees)

zt

(kilometers)

zb

(kilometers)

1 5 6 9 8 1 85 2 10

2 9 8 11 13 0.1 90 2 10

3 11 13 15 16 0.7 75 2 10

4 15 16 18 18 0.5 80 2 10

5 18 18 21 20 1.5 70 2 10

Tensile Crack
Smoothing Operator D 79

Smoothing Operator D

The surface deformation can be inverted onto fault planes that are discretized with evenly 

sized dislocations; fault tiles in the model may or may not slip (see Jonsson and others, 2002; 

Barnhart and Lohman, 2010). The observed surface displacements d can be described as a 

function g(m) of the fault model parameters m:

d m = g( ) , (138)

where the function g describes how slip across the rectangular dislocations causes displacement 

at the ground surface (Okada, 1985). The surface displacements are nonlinear functions of the 

fault- geometry parameters (for example, length, width, location, dip, and strike) and linear 

functions of the slip along the fault 

g( ) m G= s. (139)

where G is a matrix of data kernels (for example, Green’s functions from Okada, 1985), which 

depend on the assumed fault-geometry parameters and relate the surface displacement to 

the model; and sis the vector of slip values for every tile. By substituting equation 139 into 

equation 138, the surface displacement d as a linear system in matrix form

d G= s. (140)

Finally, the displacement d is weighted by using the data covariance matrix Σ

Wd WGs W W = T -1 = and Σ , (141)

where W is the weight matrix from the Cholesky decomposition of the data covariance matrix.

The system of equations relating slip on each tile to the surface displacement d consists 

of both a matrix of Green’s functions G and a smoothing matrix D

Wd

0

WG

D

s 



 



 = 



 



 κ2 , (142)

The smoothness constraint on the slip s minimizes the two-dimensional second derivative 

(Laplacian) of the fault slip. The Lagrange multiplier κ2 determines the weight of smoothing, 

and D is the second-order finite difference operator

∇ = 2s Ds . (143)

Building the Smoothing Operator

We can build the smoothing operator D for a nonuniform grid in three steps:

1. Compute the location of the center of each tile.

2. Compute the vertical and horizontal distances between neighboring tile centers assuming 

that the vertical distances ∆ ∆ h h 2 4 ( ) and are constant (all the rows of the discretized fault 

are of equal width), whereas the horizontal distances ∆ ∆ h h 1 3 ( ) and may vary (the width of 

the fault columns can change), figure 42.

3. Compute the D matrix.
80 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Computing the Location of the Center Point of Each Tile

y y z z

x x z z

z z z

f

m

m b

m

t

f

m

m b

m

t

b

m

t b

= − ( ) −

= + ( ) −

= +

cot sin

cot cos

.

δ φ

δ φ

0 5( ) − zt

, and (144)

where

x x x x

y y y y

x x

x x

m i f i

m f f i

m f

f i

= + ( ) −

= + ( ) − −

−

0 5.

, (145)

( , x y ) i i and ( , x y ) f f are the (top) starting and ending coordinates of the tile (see also fig. 43), 

δ is the dip angle, and φ is the strike

φ = −

−

















− tan 1 x x

y y

f i

f i

. (146)

men12-7629_fig42

(xm, ym, zm ∆h4 ) ∆h2

∆h1

(xi, yi) (xf, yf)

f f f

∆h3

∆h1 and ∆h3 −̶ column width

∆h2 and ∆h4 − row height

(xm, ym, zm) − center point of tile f f f

(xi, yi

) − starting coordinates 

(xf, yf

) − ending coordinates

Figure 42. Grid and tile geometry for smoothing matrix D . 
Smoothing Operator D 81

Computing Distances Between the Center Points of Tiles

1. The vertical distance is

∆h x x y y z z f

m

f

m

f

m

f

m

b

m

b

m

2

2 2 2

= − ( ) ( ) 2 1( ) + − ( ) ( ) 2 1( ) + − ( ) ( ) 2 1( ) , (147)

where (1) and (2) indicate the component of the vectors; 

2. The horizontal distance is

a. at the left edge (first column, where n = 1),

∆h P x P x y P y z P z f

m

f

m

f

m

f

m

b

m

1 b

2 2

( : 1 1 ) ( = + ( ) ) ( − 1 1 ) ( + + ( ) ) ( − 1 1 ) ( + + ) − m ( )1

2

( ) ; (148)

b. at the right edge (last column, where n N = ),

∆h NP P NP

x N f P P x NP y NP P y NP m

f

m

f

m

f

m

3

2 2

( : 1 )

( ) ( ) ( ) ( )

+ − =

= − ( ) − + − ( ) − + z N b P P z NP m

b

m ( ) ( ) − − ( )

2

, (149)

where P is the number of rows of tiles, N is the number of columns of tiles, and D is a square matrix with the number of rows 

number of columns PN; and

c. from the central columns (where 2 1 ≤ ≤ n N − ),

∆h n P nP

x n f P P x nP y nP P y nP m

f

m

f

m

f

m

1

2

( ) 1 1 + − =

= − ( ) − + − ( ) −

( ) :

( ) ( ) ( ) ( )

2 2

3 1 1

+ − ( ) −

( ) + − =

= + −

z nP P z nP

h n P nP

x nP P x

b

m

b

m

f

m

f

m

( ) ( )

( ) :

( )

 and

∆

( ) nP y n( ) P P y n( ) P z ( ) nP P z ( ) nP . f

m

f

m

b

m

b

m ( ) + + ( ) − + + ( ) − 2 2 2

(150)

men12-7629_fig43

C. Matrix representation of 

 smoothing operator D

5 6 7 8 9 10 11 12 13 14 −2

0

2

4

0

2

4

6

8

10

X, in kilometers

Depth, in kilometers

A B. Fault discretized into 18 tiles . Grid

2 4 6 8 10 12 14 16 18

2

4

6

8

10

12

14

16

18

(xm, ym, zm ∆h4 ) ∆h2

∆h1

(xi, yi) (xf, yf)

f f f

∆h3

∆h1 and ∆h3 ̶− column width

∆h2 and ∆h4 − row height

(xm, ym, zm) − center point of tile f f f

(xi, yi) − starting coordinates 

(xf, yf) − ending coordinates Y, in kilometers

Columns

Rows

Figure 43. A, Grid used to compute the second-order finite-difference operator D; B, fault divided into P=2 rows, N=9 columns, and 

18 tiles of different width; C, smoothing operator D for the fault in B. D has ( ) NP = 2 324 elements, but the plot shows only the nonzero 

elements of the matrix. Blue, central diagonal; red, first diagonal above and below the central diagonal; green, Pth diagonal above and 

below the central diagonal.
82 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Computing the Smoothing Operator D

The smoothing operator D is a second-order finite difference operator (equation 143). 

As the slip is determined on discrete tiles (figs. 42 and 43), ∇2

s = Ds is a two-dimensional, 

second-order finite-difference sum

s s s

h h

s s s

h

i j − + − + i j i j ij ij ij − +

( ) +

+

1 1 − +

1 3

2

1 1

2

2

2 2

4

, , , , , ,

∆ ∆ ∆

. (151)

The matrix implementing D is always zero, except along a limited number of diagonals 

(fig. 43C): the central diagonal and the first and Pth diagonals above and below the central 

diagonal, where P is the number of rows of tiles on the fault (fig. 43B). 

Case 1—Fault does not Break the Surface

1. First, compute the central diagonal (0)

d

h h h = − +















 2 1 1

1 3 2

2 ∆ ∆ ∆

, (152)

2. diagonal (+1)

dp h 1 dp N P dp dp NP

2

2 1 1 1

1 = ( ) ( ) 1 0 = = 1 1− ∆

, : , ( : ) , (153)

3. diagonal (-1)

dm

h 1 dm N P dm dm NP

2

2 1 1 1

1

= + ( ) 1 1 = = 0 2

∆ , ( : ) ( : ) , (154)

4. diagonal (+P) 

dp 2

h 4 4 = dp dp4 1 NP P ( ) + = − ( ) ∆ ∆ h h 3 1 ∆ 3

, : , and finally (155)

5. diagonal (-P) 

( ) 4 4 4 ( )

11 3

2 dm , 1 : dm dm P NP

hh h = = +∆ ∆ +∆ . (156)
Smoothing Operator D 83

Case 2—Fault Breaks the Surface

1. First, compute the central diagonal (0)

( ) (( ) )

2 2 (( ) )

1 1 1 , 1: 1 1 (1) 1: 1 1

d dNP

h h NP = − +=

∆ ∆ −+ , (157)

2. diagonal (+1)

( ) (( ) ) (( ) ) 1 1

2 2

1 1 1 , 1: 1 1 (1) 1: 1 1 dp dp N P h h NP = − − + =−

∆ ∆ −+ , (158)

3. diagonal (-1)

dm N P 1 ( ) ( ) 1: 1 0 − = , (159)

4. diagonal (+P) — this condition allows slip in the upper row of tiles to be nonzero when the 

fault breaks the surface (that is, produces measurable surface offsets)

dp dp N P 4 4 ( ) ( ) 1 0, 1: 2 1 0 = ( − +=) , and finally (160)

5. diagonal (-P) — see note for step 4.

dm4 4 ( ) 1 0, 1: 2 1 0 = dm N P (( ) − +=) . (161)

Case 3—Dike Opening at the Base of the Dislocation

1. First, compute the central diagonal (0)

( ) ( )

2 ( ) ( )

1 1: 

1:

d NP

h NP = ∆ , (162)

2. diagonal (+1)

 1: 1 0 dp N P 1 ( ) ( ) − = , (163)

3. diagonal (-1)

( ) ( )

(( ) ) (( ) )

1

2

1

2

1 1 , 1

1 1: 1 1

1: 1 1

dm P

h P

dm N P P

h N PP

− =− ∆ −

− + − =−

∆ − +−

(164)
84 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

4. diagonal (+P) — this condition allows slip in the lower row of tiles to be nonzero when the dike opens at the bottom 

dp N P 4 ( ) ( ) 1: 1 0 − = , and finally (165)

5. compute diagonal (-P) — see note for step 4.

dm N P 4 ( ) ( ) 1: 1 0 − = . (166)

D Matrix

( ) ( ) ( ) ( ) ( ) 11 4 4 D = + + −+ + − diag d diag dp diag dm diag dp P diag dm P ,1 , 1 , , , (167)

where diag v k ( ) , is a function that returns a square matrix of order n k + with the elements of v on the kth diagonal when v is 

a vector of n components; k = 0 represents the main diagonal, k > 0 a diagonal above the main diagonal, and k < 0 a diagonal 

below the main diagonal. 

Verification

The algorithm for the smoothing operator D was tested by using the fault geometry shown in figure 44 and table 36. 

Following equations 143 and 151, we applied D to the vector of fault slips and compared the results to 2 ∇ s (fig. 45), which is 

computed numerically by using finite differences (appendix 1).

4 6 8 10 12 14 16 18 20 22

5

10

15

20

0

2

4

6

8

10

Strike-slip fault discretized into 368 tiles

X, in kilometers

Depth, in kilometers

Y, in kilometers

Figure 44. Fault geometry used to validate the algorithm for the 

smoothing operator D. See table 36 for parameters.
Smoothing Operator D 85

Table 36. Values of fault parameters. 

[Fault geometry is for the right-lateral strike-slip fault shown in figure 44. Segments were divided into tiles (or patches) by using dL = [2.0 0.5 3.0 1.0 1.5]

kilometers and dW = 0.5 kilometers. xi, yi, xf, and yf, coordinates of the initial and final positions along the top edges of the fault segments; U1, final slip of 

segment; zt and zb, coordinates of the top and bottom of the segment] 

Segment xi

(kilometers) 

yi

(kilometers)

xf

(kilometers)

yf

(kilometers)

U1

(meters)b

Dip angle

(degrees)

zt

(kilometers)

zb

(kilometers)

1 5 6 9 8 1 90 2 10

2 9 8 11 13 0.5 90 2 10

3 11 13 15 16 1.5 90 2 10

4 15 16 18 20 0.5 90 2 10

5 18 20 21 20 0.75 90 2 10

−5

5 10 15

6

8

10

12

14

16

18

20

Slip distribution

0.5

1

1.5

Depth, in kilometers zt

zb

xi xf

X, in kilometers

xi xf xi xf

5 10 15

6

8

10

12

14

16

18

20 −6

−5

−4

−3

−2

−1

0

Second derivative x 10

5 10 15

−8

−6

−4

−2

0

2

4

6

8

x 10 −22

6

8

10

12

14

16

18

20

∆2

s - DU ∆2s

Figure 45. Color maps for the smoothing operator, slip distribution, and finite difference second derivative 2 ∇ s (appendix 1) for the 

fault geometry shown in figure 44. The smoothing operator D is a diagonal matrix with a number of elements equal to the square of the 

number of tiles (fig. 44). The difference between the finite difference second derivative 2 ∇ s and DU is approximately 10-22. xi and xf , 

horizontal east coordinate for the left and right corner of the fault; zt , depth of the top edge of the fault; zb , depth of the bottom edge of 

the fault.
86 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Example for a Right-Lateral Strike-Slip Fault

The example is based on the fault geometry shown below (fig. 46 and table 37). Numerical values of the smoothing 

operator D, slip vector U, product DU, and finite-difference second derivative (grad)2s are in tables 38–41.

Figure 46. Fault geometry employed in the example. See table 37

for fault parameters.

4 6 8 10 12 14 −2

0

2

4

0

2

4

6

8

10

Fault discretized into 10 tiles

X, in kilometers

Depth, in kilometers

Y, in kilometers

Table 37. Values of fault parameters.

[Fault geometry is for the right-lateral strike-slip fault shown in figure 46. Segments were divided into tiles by using dL = [ ] 3.0 1.0 3.0 kilometers and dW = 5.0

kilometers. xi, yi, xf, and yf, coordinates of the initial and final positions along the top edges of the fault segments; U3, final opening dislocation of segment; zt and 

zb, coordinates of the top and bottom of the segment]

Segmenta xi

(kilometers) 

yi

(kilometers)

xf

 (kilometers)

yf

(kilometers)

U1

(meters) b

Dip angle

(degrees)

zt

(kilometers)

zb

(kilometers)

1 5 4 9 4 1 60 2 10

2 9 4 11 4 0.5 60 2 10

3 11 4 14 4 1.5 60 2 10
Smoothing Operator D 87

Table 38. Values of the smoothing operator D. 

[The reader should multiply the elements below by 10-5. Value of D is in (kilometers)-2]

1 2 3 4 5 6 7 8 9 10

1 -0.0594 0.0047 0.0250 0 0 0 0 0 0 0

2 0.0047 -0.0594 0 0.0250 0 0 0 0 0 0

3 0.0286 0 -0.0760 0.0047 0.0381 0 0 0 0 0

4 0 0.0286 0.0047 -0.0760 0 0.0381 0 0 0 0

5 0 0 0.0533 0 -0.1427 0.0047 0.0800 0 0 0

6 0 0 0 0.0533 0.0047 -0.1427 0 0.0800 0 0

7 0 0 0 0 0.0667 0 -0.1094 0.0047 0.0333 0

8 0 0 0 0 0 0.0667 0.0047 -0.1094 0 0.0333

9 0 0 0 0 0 0 0.0250 0 -0.0594 0.0047

10 0 0 0 0 0 0 0 0.0250 0.0047 -0.0594

Table 39. Values of the slip vector U.

[Values of U are in meters] 

1 2 3 4 5 6 7 8 9 10

1 1.0 1.0 1.0 1.0 0.5 0.5 0.5 0.5 1.5 1.5

Table 40. Values of the product DU. 

[Values of DU are in meters per square kilometer]. 

1 2 3 4 5

1 -0.2969 -0.2374 0.2432 0.3099 -0.6953

2 -0.2969 -0.2374 0.2432 0.3099 -0.6953

Table 41. Values of the finite-difference second derivative 2 ∇ s .

[Values 2 ∇ s of are in meters per square kilometer. See Appendix 1.] 

1 2 3 4 5

1 -0.2969 -0.2374 0.2432 0.3099 -0.6953

2 -0.2969 -0.2374 0.2432 0.3099 -0.6953
88 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Summary

This report is a practical reference for readers interested in applying mathematical models 

to investigate volcano and earthquake physics. The report emphasizes the application of 

analytical models to fit and interpret ground deformation. Although analytical models are based 

on numerous simplifications (for example, the assumption that the crust is a homogenous, 

isotropic, elastic, flat half-space) that make the set of differential equations describing the 

problem tractable, they can take into account a vast array of source geometries. The careful use 

of analytical models, together with high-quality data sets, can yield valuable insights into the 

nature of the deformation source. 

The paper introduces analytical formulation that describe the deformation caused by 

magma chambers, magma conduit, sill, dikes and faults, and provides tables with numerical 

examples and Matlab scripts to aid in the use and coding of the formulas. All formulas have 

been checked for errors and verified against numerical models (COULOMB 3.3; http://

earthquake.usgs.gov/research/modeling/coulomb/; FEM models). A zipped file with MATLAB 

scripts can be downloaded from http://pubs.usgs.gov/tm/13/b1. 

Suggestions about new models that could be included in future revisions, reports of 

existing bugs in the Matlab scripts (or errors in the formulas) and comments on how to improve 

this publication will be greatly appreciated. 

References Cited

Amoruso, Antonella, and Crescentini, Luca, 2009, Shape and volume change of pressurized 

ellipsoidal cavities from deformation and seismic data: Journal of Geophysical Research, 

v. 114, B02210 (also available at http://dx.doi.org/ 10.1029/2008JB005946).

Barnhart, W.D., and Lohman, R.B., 2010, Automated fault model discretization for inversions 

for coseismic slip distributions: Journal of Geophysical Research, v. 115, B10419, 17 p. (also 

available at http://dx.doi.org/ 10.1029/2010JB007545).

Cervelli, P.F., Fournier, T.J., Freymueller, J.T., and Power, J.A., 2006, Ground deformation 

associated with the precursory unrest and early phases of the January 2006 eruption of 

Augustine Volcano, Alaska: Geophysical Research Letters, v. 33, L18304 (also available at 

http://dx.doi.org/ 10.1029/2006GL027219).

Dawson, John, and Woods, Alex, 2010, ITRF to GDA94 Coordinate Transformations: 

Journal of Applied Geodesy, v. 4, 189–199 (also available at http://dx.doi.org/ 10.1515/

JAG.2010.019)

Defense Mapping Agency (DMA), 1989, The Universal Grids—Universal Transverse Mercator 

(UTM) and Universal Polar Stereographic (UPS): Fairfax, Va., DMA Technical Manual 

DMATM 8358.2.

Dzurisin, Daniel, 2003, A comprehensive approach to monitoring volcano deformation as a 

window on the eruption cycle: Review of Geophysics, v. 41, no. 1, 1001 p. (also available at 

http://dx.doi.org/10.1029/2001RG000107).

Fialko, Y., Khazan, Y., and Simons, M., 2001, Deformation due to a pressurized horizontal 

circular crack in an elastic half-space, with applications to volcano geodesy: Geophysical 

Journal International, v. 146, no.1, p. 181–190.

Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., 1997, Global positioning system—

Theory and practice (4th revised ed.): New York and Vienna, Springer, 389 p.
References Cited 89

Jónsson, S., Zebker, H., Segall, P., and Amelung, F., 2002, Fault slip distribution of the 

1999 M7.2 Hector Mine earthquake, California, estimated from satellite radar and GPS 

measurements: Bulletin of the Seismological Society of America, v. 92, no. 4, p. 1377–1389.

Masterlark, Timothy, 2007, Magma instrusion and deformation predictions—Sensitivities to 

the Mogi assumptions: Journal of Geophysical Research, v. 112, no. B06419, 17 p. (also 

available at http://dx.doi.org/10.1029/2006JB004860).

McTigue, D.F., 1987, Elastic stress and deformation near a finite spherical magma body—

Resolution of the point source paradox: Journal of Geophysical Research, v. 92, no. B12, p. 

12931–12940 (also available at http://dx.doi.org/10.1029/JB092iB12p12931).

National Imagery and Mapping Agency (NIMA), 2000, Department of Defense World Geodetic 

System 1984—Its definition and relationship with local geodetic systems: Technical Report 

TR8350.2.

National Imagery and Mapping Agency (NIMA), 2004, Implementation of the World Geodetic 

System 1984 (WGS 84) Reference Frame G1150. Addendum to NIMA TR 8350.2.

Newman, A.V., Dixon, T.H., and Gourmelen, N., 2006, A four-dimensional viscoelastic 

deformation model for Long Valley Caldera, California, between 1995 and 2000: Journal of 

Volcanology and Geothermal Research, v. 150, p. 244–269.

Okada, Yoshimitsu, 1985, Surface deformation due to shear and tensile faults in a half-space: 

Bulletin of the Seismological Society of America, v. 75, no. 4, p. 1135–1154.

Okada, Yoshimitsu, 1992, Internal deformation due to shear and tensile faults in a half-space: 

Bulletin of the Seismological Society of America, v. 82, no. 2, p. 1018–1040.

Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 1992, Numerical recipes 

in Fortran 77—The Art of Scientific Computing (2d ed.), v. 1: Cambridge, Cambridge 

University Press, 921 p.

Segall, Paul, 2010, Earthquake and volcano deformation: Princeton, N.J., Princeton University 

Press, 458 p.

Tiampo, K.F., Rundle, J.B., Fernandez, J., and Langbein, J.O., 2000, Spherical and ellipsoidal 

volcanic sources at Long Valley caldera, California, using a genetic algorithm inversion 

technique: Journal of Volcanology and Geothermal Research, v. 102, no. 3, p. 189–206.

U.S. Army, 1987, Standard UTM-to-GEO/GEO-to-UTM Coordinate Conversion Algorithm: 

Technical Manual TM 5–241–8 (accessed on August 2, 2012 at http://ntrs.nasa.gov/archive/

nasa/casi.ntrs.nasa.gov/19890002875_1989002875.pdf).

Yang, X.-M., Davis, P.M., and Dieterich, J.H., 1988, Deformation from inflation of a dipping 

finite prolate spheroid in an elastic half-space as a model for volcanic stressing: Journal of 

Geophysical Research, v. 93, no. B5, p. 4249–4257.
90 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Appendix 1. Mathematical Methods for 

Computing Displacement

Displacement Gradient Tensor

The partial differentiation of the displacement vector u with respect to the spatial 

coordinates yields the displacement gradient tensor ∇u. The finite strains measured by a 

Gladwin Tensor Strainmeter (GTSM) can be expressed as a function of the displacementgradient tensor. We can compute the components of the displacement gradient tensor using the 

finite-difference derivatives (Press and others, 1992, p. 180):

∂

∂ = ( ) + − − ( ) +

∂

∂ = ( ) + −

u

x

u x h y z u x h y z

h

o h

u

y

u x y h z u x

x x x

x x x

,, ,, ( )

, ,

2

2

, , ( )

,, ,, ( )

y h z

h

o h

u

x

u x h y z u x h y z

h

o h

u

y y y

y

( ) − +

∂

∂ = ( ) + − − ( ) +

∂

∂

2

2

2

2

y

u x y h z u x y h z

h

o h y y = ( ) + − − ( ) + , , , , ( ) 2

2

(1-1)

The horizontal derivatives of the vertical deformation correspond to the east and north 

components of tilt:

τ

τ

E z z z

N z z

u

x

u x h y z u x h y z

h

o h

u

y

u x y h z

= ∂

∂ = ( ) + − − ( ) +

= ∂

∂ = +

,, ,, ( )

, ,

2

2

( ) − − ( ) +

u x y h z

h

o h z , , ( ) 2

2

. (1-2)

Length Change to an Electronic-Distance-Meter Base Line 

The change in length of an electronic-distance-meter (EDM) baseline (δrab) is given 

by the sum of the projections of the radial deformations at the two benchmarks along the 

baseline (fig. 1-1)

δ θ r u ab = + ra cos c a r u b b osθ , (1-3)

where

cos c θ θ os a a ab b

a ab

b

b ab a

b ab

r r r

r r

r r r

r r = + − = + − 2 2 2 2 2 2

2 2

and . (1-4)

It is also possible to write equation 1-3 in the vector notation

δr u r

r u r

r ab ra ab

ab

rb

ab

ab

= ⋅ + ⋅        

. (1-5)
Appendix 1 91

Second Derivative of the Slip Vector

We use this formulation for second derivative of the slip vector s (after Zwillinger, 1997, 

p. 661) to test the smoothing operator D. The second derivative of s is defined as

∇ = ∂

∂

+

∂

∂

+

∂

∂

















2

2

2

2

2

2

2 s s

x y z , (1-6)

where 

 s is a matrix describing the slip distribution along the fault. 

The irregular grid and index convention used are shown in figure 1-2 below.

Figure 1-1. Geometry of a length change to an electronicdistance-meter (EDM) baseline. x y a a ( ) , , Cartesian coordinates

of the EDM base station; x y b b ( ) , , Cartesian coordinates of the

EDM target station; x y s s ( ) , , Cartesian coordinates of source

location; rab, EDM baseline length; ra, Cartesian distance between 

EDM base and source; rb , Cartesian distance between EDM target and 

source; θa and θb, angles; ura and urb , radial deformation.

(xa,ya) uracosθa

(xs,ys)

(xb,yb)

ura urb

ra

θa θb



urbcosθb

rb

EDM baseline = rab

∆h4 (i, j) ∆h2

∆h1

∆h3

∆h1 and ∆h3 − column width

∆h2 and ∆h4 − row height

(i+1, j)

(i-1, j)

(i, j-1) (i, j+1)

i is the row index

j is the column index

Figure 1-2. Irregular grid geometry, after Zwillinger (1997, p. 661).
92 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

The Laplacian is computed in three steps:

Step 1. Compute the Location of the Center Point of Each Tile

( )

( )

( )

cot cos , East

cot sin , North

0.5 , Up

m m

f m bt

m m

f m bt

m

b t bt

x x zz

y y zz

z z zz

=+ − δφ

= − − δφ

=+ −

, (1-7)

where 

x x x x y y y y

x x

x x

m i f i m f f i

m f

f i

= + ( ) − = + ( ) − −

− 0 5. , . (1-8)

Step 2. Build Matrices with East, North, and Up Components of the Center 

and the Slip of the Tile

( ) ( )

( ) ( )

( ) ( )

( ) ( )

index

– index

1 : , 1 : East

1 : , 1 : North

1 : , 1 : Up

1 : , 1 : Slip

m

f

m

f

m

b

J JP

K KP

MX J K n x J K

MY J K n y J K

MZ J K n z J K

MU J K n U J K

= +

=

+ =+

+ =+

+ =+

+ =+

, (1-9)

where n N = 1,, .

Step 3. Compute the Laplacian L

Internal Tiles (i P = − 2 1 , ,  and j N = − 2 1 , ,  ),

∆

∆

h MX MX MY MY MZ MZ

h

4 1 i j i j i j i j i j i j

2

1

2

1

2

3

= − ( ) + − ( ) + − ( )

=

, , − − , , , , −

MX MX MY MY MZ MZ

h MX

ij ij ij ij ij ij

i

( ) , , + + − + − ( ) , , + − ( ) , , +

=

1

2

1

2

1

2

∆ 2 ( ) + + − + − ( ) + − ( ) +

= −

1

2

1

2

1

2

, , , , , ,

,

j i j i j i j i j i j

i j

MX MY MY MZ MZ

∆h M( ) X MX M i j , , − − 1 + − ( ) Y M ij ij Y M , , + − ( ) Z M ij ij Z , −

2

1

2

1

2

, (1-10)

and

( ) ( ) ( )

( )

, 1 1, , 1

,

11 3 2 2 4 31 3

1,

,

42 4 13 24

2 

1 1 

i j i j i j

i j

i j

i j

MU MU MU

L

hh h h h h hh h

MU

MU

h h h hh h h

−+ +

−



= + ++ 

∆ ∆ +∆ ∆ ∆ +∆ ∆ ∆ +∆ 

 

+− +   ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆ 





. (1-11)
Appendix 1 93

Top edge (i = 1) where the Fault does not Break the Earth’s Surface

Left tile ( j =1)

∆

∆

h MX MX MY MY MZ MZ

h

2 1 i j i j i j i j i j i j

2

1

2

1

2

4

= − ( ) + − ( ) + − ( ) =

+ + , , , , + , ,

∆

∆

h

h MX M ij ij ij ij ij ij X MY MY MZ MZ

2

3 1

2

1

2

= − ( ) , , + + + − ( ) , , + − ( ) , , +1

2

∆ ∆ 1 3 h h =

, (1-12)

and

 ( ) ( )

1, , 1

, ,

2 2 4 31 3 13 2 2

1 1 2 i j i j

i j i j

MU MU

L MU

h h h hh h hh h h

+ +    

= +− +     ∆ ∆ +∆ ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆  

. (1-13)

Internal tiles ( j N = − 2, , 1  )

∆

∆

h MX MX MY MY MZ MZ

h

2 1 i j i j i j i j i j i j

2

1

2

1

2

4

= − ( ) + − ( ) + − ( ) =

+ + , , , , + , ,

∆

∆

∆

h

h MX MX MY MY MZ MZ

h

ij ij ij ij ij ij

2

3 1

2

1

2 2

= − ( ) , , + + + − ( ) , , + − ( ) , ,

1 1

2

1

2

1

2

= − ( ) MX MX MY MY MZ MZ ij ij , , − − + − ( ) ij ij , , + − ( ) ij ij , , −

’ (1-14)

and

( ) ( ) ( )

, 1 1, , 1

, ,

11 3 2 2 4 31 3 13 2 2

1 1 2 i j i j i j

i j i j

MU MU MU

L MU

hh h h h h hh h hh h h

−+ +    

= + +− +     ∆ ∆ +∆ ∆ ∆ +∆ ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆  

. (1-15)

 Right tile ( j N= )

∆

∆

h MX MX MY MY MZ MZ

h

2 1 i j i j i j i j i j i j

2

1

2

1

2

4

= − ( ) + − ( ) + − ( ) =

+ + , , , , + , ,

∆

∆

h

h MX M ij ij X MY M ij ij Y MZ M ij ij Z

2

1 1

2

1

2

= − ( ) , , − − + − ( ) , , + − ( ) , , −1

2

∆ ∆ 3 1 h h =

, (1-16)

and

( ) ( )

, 1 1,

, ,

11 3 2 2 4 13 2 4

1 1 2 i j i j

i j i j

MU MU

L MU

hh h h h h hh h h

− +    

= +− +     ∆ ∆ +∆ ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆  

. (1-17)
94 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Top Edge (i = 1) where the Fault Breaks the Earth’s Surface

If the fault breaks the surface, we must replace equations 1-12 to 1-17 with

∆h MX MX MY MY MZ MZ

L

i j i j i j i j i j i j

i j

2 1

2

1

2

1

2

= − ( ) + + , , + − ( ) , , + − ( ) + , ,

, = − = MU MU + j N h

i j 1 i j

2

1 , , , , , ∆ 

. (1-18)

Bottom Edge (i P = ),Tensile Crack (Dike) is not Opening

Left tile (j = 1)

∆

∆

h MX MX MY MY MZ MZ

h

4 1 i j i j i j i j i j i j

2

1

2

1

2

3

= − ( ) + − ( ) + − ( )

=

, , − − , , , , −

MX MX MY MY MZ MZ

h h

( ) ij ij ij ij ij ij , , + + − + − ( ) , , + − ( ) , , +

=

1

2

1

2

1

2

∆ ∆ 2 4

∆ ∆ h h 1 3 =

, (1-19)

( ) ( )

, 1 1,

, ,

31 3 4 2 4 13 24

1 1 2 i j i j

i j i j

MU MU

L MU

hh h h h h hh hh

+ −    

= + −+     ∆ ∆ +∆ ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆  

. (1-20)

Internal tiles ( j N = − 2 1 , ,  )

∆

∆

h MX MX MY MY MZ MZ

h

4 1 i j i j i j i j i j i j

2

1

2

1

2

3

= − ( ) + − ( ) + − ( )

=

, , − − , , , , −

MX MX MY MY MZ MZ

h h

( ) ij ij ij ij ij ij , , + + − + − ( ) , , + − ( ) , , +

=

1

2

1

2

1

2

∆ ∆ 2 4

1 1

2

1

2

1

2

∆h M = − ( ) X M ij ij , , X M − − + − ( ) Y M ij ij , , Y M + − ( ) Z M ij ij , , Z −

, (1-21)

( ) ( ) ( )

, 1 , 1 1,

, ,

11 3 31 3 4 2 4 13 24

1 1 2 i j i j i j

i j i j

MU MU MU

L MU

hh h hh h h h h hh h h

− +−    

= ++ − +     ∆ ∆ +∆ ∆ ∆ +∆ ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆  

. (1-22)
Appendix 1 95

Right tile ( j N = )

∆

∆

h MX MX MY MY MZ MZ

h

4 1 i j i j i j i j i j i j

2

1

2

1

2

1

= − ( ) + − ( ) + − ( )

=

, , − − , , , , −

MX MX MY MY MZ MZ

h h

( ) ij ij , , − + − ( ) ij ij , , + − ( ) ij ij , ,

=

− − 1 −

2

1

2

1

2

∆ ∆ 3 1

∆ ∆ h h 2 4 =

, (1–23)

( ) ( )

, 1 1,

, ,

11 3 4 2 4 13 2 4

1 1 2 i j i j

i j i j

MU MU

L MU

hh h h h h hh h h

− −    

= +− +     ∆ ∆ +∆ ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆  

. (1-24)

Bottom Edge (i P = ), the Tensile Crack (Dike) is Opening

If the crack is opening, we must replace all the equations above with

∆h MX MX MY MY MZ MZ

L

i j i j i j i j i j i j

i j

4 1

2

1

2

1

2

= − ( ) , , − − + − ( ) , , + − ( ) , , −

, = − = ... MU MU − j N h

i j , , i j , , , 1

4

1

∆

. (1-25)

Left Edge (i P = − 2 1 , ,  and j = 1)

∆

∆

h MX MX MY MY MZ MZ

h

4 1 i j i j i j i j i j i j

2

1

2

1

2

3

= − ( ) + − ( ) + − ( )

=

, , − − , , , , −

MX MX MY MY MZ MZ

h MX

ij ij ij ij ij ij

i

( ) , , + + − + − ( ) , , + − ( ) , , +

=

1

2

1

2

1

2

∆ 2 ( ) + + − + − ( ) + − ( ) +

=

1

2

1

2

1

2

1 3

, , j i MX MY MY MZ MZ j i , , j i j i , , j i j

∆ ∆ h h 

, (1-26)

( ) ( ) ( )

1, , 1 1,

, ,

2 2 4 31 3 4 2 4 13 24

1 1 2 i j i j i j

i j i j

MU MU MU

L MU

h h h hh h h h h hh h h

+ +−    

= ++ − +     ∆ ∆ +∆ ∆ ∆ +∆ ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆  

(1-27)
96 Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models

Right Edge (i P = − 2 1 , ,  and j N = )

∆

∆

h MX MX MY MY MZ MZ

h

1 1 ij ij ij ij ij ij

2

1

2

1

2

2

= − ( ) + − ( ) + − ( )

=

, , − − , , , , −

MX MX MY MY MZ MZ

h h

( ) i j + + − i j + − ( ) i j i j + − ( ) i j + i j

=

1

2

1

2

1

2

3 1

, , , , , ,

∆ ∆

∆h M 4 1 X M i j X M i j Y M i j Y M i j Z M i j Zi j

2

1

2

1

2

= − ( ) , , − − + − ( ) , , + − ( ) , , −

, (1-28)

( ) ( )

( )

, 1 1,

,

11 3 2 2 4

1,

,

42 4 13 24

2

1 1

i j i j

i j

i j

i j

MU MU

L

hh h h h h

MU

MU

h h h hh h h

− +

−



= ++ 

∆ ∆ +∆ ∆ ∆ +∆ 

 

+− +   ∆ ∆ +∆   ∆ ⋅∆ ∆ ⋅∆ 

. (1-29)

References Cited

Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 1992, Numerical recipes in 

Fortran 77—The Art of Scientific Computing (2nd ed.), v. 1: Cambridge University Press, 

New York, N.Y., p. 921.

Zwillinger, Daniel, 1997, Handbook of differential equations (3d ed.): Academic Press, p. 801.
Publishing support provided by the U.S. Geological Survey

Publishing Network, Tacoma Publishing Service Center

For additional information contact:

Director, Menlo Park Science Center

U.S. Geological Survey

345 Middlefield Road MS 977

Menlo Park, California 94025

http://earthquake.usgs.gov/regional/nca/
Battaglia and others— Modeling Crustal Deformation near Active Faults and Volcanic Centers—A Catalog of Deformation Models—TM 13-B1