InjectSills.jl

This collects analytical solutions that give the displacement around sills, dikes, and magmatic intrusions embedded in an elastic half-space. It can be combined with GeophysicalModelGenerator to plot surface displacement caused by intrusion events, and with JustPIC to advect particles during sill injection.

Available source types

All source types are subtypes of AbstractSill{N,T} and share a common interface: hostrock_displacement, inside, and isdimensional. Parameters are stored as GeoParams units and can be non-dimensionalized.

PennyShapedSill — penny-shaped tensile crack (2D / 3D)

Displacement field around a pressurized penny-shaped crack (Sun, 1969). A 2D sill is the axisymmetric section of the penny: the 3D solution in the vertical plane through the sill axis.

sill = PennyShapedSill(
    Center = Point3(0.0, 0.0, -5000.0)*m,
    R      = 2000.0m,     # radius
    H      = 100.0m,      # maximum opening thickness
    E      = 1.5e10Pa,
    ν      = 0.3*NoUnits,
    Angle  = Vec2(10.0, 0.0),   # dip, strike [degrees]
)
d = hostrock_displacement(sill, Point3(x, y, z))

Reference: Sun, R.J. (1969): Theoretical size of hydraulically induced horizontal fractures and corresponding surface uplift in an idealized medium. J. Geophys. Res. 74, 5995–6011. https://doi.org/10.1029/JB074i025p05995

---

PlaneStrainSill — pressurized crack in plane strain (2D)

Sill in a Cartesian 2D model: a crack of half-length R with the elliptical opening H√(1 - x²/R²), where H = 4(1-ν²) ΔP R/E. Its cross-sectional area is Q = π R H/2; half of it crosses every line parallel to the sill on either side. Takes the same keywords as the 2D PennyShapedSill, with Q an area.

sill = PlaneStrainSill(
    Center = Point2(0.0, -5000.0)*m,
    R      = 2000.0m,     # half-length
    H      = 10.0m,       # maximum opening
    E      = 1.5e10Pa,
    ν      = 0.3*NoUnits,
    Angle  = Vec1(10.0),  # dip [degrees]
)
d = hostrock_displacement(sill, Point2(x, z))

References: Westergaard, H.M. (1939): Bearing pressures and cracks. J. Appl. Mech. 6, A49–A53. Pollard, D.D. & Segall, P. (1987): Theoretical displacements and stresses near fractures in rock. In Atkinson, B.K. (ed.), Fracture Mechanics of Rock, Academic Press, 277–349.

---

MogiSphere — Mogi point-source sphere (2D / 3D)

First-order solution for a pressurized spherical cavity; displacement is proportional to (p − Center) / |p − Center|³.

src = MogiSphere(
    Center = Point3(0.0, 0.0, -5000.0)*m,
    r      = 1500.0m,
    ΔP     = 10e6Pa,
    G      = 10e9Pa,
    ν      = 0.25*NoUnits,
)
d = hostrock_displacement(src, Point3(x, y, z))

Reference: Mogi, K. (1958): Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bull. Earthq. Res. Inst. Univ. Tokyo 36, 99–134.

---

McTigueSphere — McTigue finite-radius sphere (2D / 3D)

Same geometry as MogiSphere with a higher-order (r/d)³ correction for finite sphere size.

src = McTigueSphere(
    Center = Point3(0.0, 0.0, -5000.0)*m,
    r      = 1500.0m,
    ΔP     = 10e6Pa,
    G      = 10e9Pa,
    ν      = 0.25*NoUnits,
)
d = hostrock_displacement(src, Point3(x, y, z))

Reference: McTigue, D.F. (1987): Elastic stress and deformation near a finite spherical magma body: resolution of the point source paradox. J. Geophys. Res. 92 (B12), 12931–12940. https://doi.org/10.1029/JB092iB12p12931

---

EllipticalIntrusion — kinematic elliptical intrusion (2D / 3D)

Finite displacement that opens an ellipse (2D) or spheroid (3D) of width W and thickness H by growing similar ellipses/spheroids from the center with the cube rule. It conserves volume in 3D and in an axisymmetric 2D section, but not area in Cartesian 2D. Its far-field footprint scales like depth × W/H; for thin sills under a free surface, use PennyShapedSill (3D) or PlaneStrainSill (Cartesian 2D).

intr = EllipticalIntrusion(Center=Point2(0.0, -5000.0)*m, Angle=Vec1(0.0)*NoUnits, W=2000.0m, H=100.0m)
d = hostrock_displacement(intr, Point2(x, z))

---

FiniteEllipsoidalCavity — finite ellipsoidal cavity (3D)

Surface deformation from a uniformly-pressurized triaxial ellipsoidal cavity. Uses an adaptive grid of point Compound Dislocation Models (pCDMs) internally. The array form is efficient (grid built once per struct).

fec = FiniteEllipsoidalCavity(
    Center = Point3(0.0, 0.0, -10000.0)*m,
    ax     = 500.0m,  ay = 500.0m,  az = 2000.0m,
    Angle  = Vec{3}(0.0, 25.0, -10.0)*NoUnits,  # omegaX, omegaY, omegaZ [degrees]
    ΔP     = 10e6Pa,
    mu     = 10e9Pa,
    lambda = 10e9Pa,
)
# Array of surface observation points (efficient):
ue, un, uv, dV, DV, Ns = hostrock_displacement(fec, X, Y)

# Single point:
d = hostrock_displacement(fec, Point3(x, y, 0.0))

References:

• Nikkhoo, M., Rivalta, E. (2023): Surface deformations and gravity changes caused by pressurized finite ellipsoidal cavities. Geophys. J. Int. https://doi.org/10.1093/gji/ggac351
• Nikkhoo, M., Rivalta, E. (2022): Analytical solutions for gravity changes caused by triaxial volumetric sources. Geophys. Res. Lett. https://doi.org/10.1029/2021GL095442
• Nikkhoo, M., Walter, T.R., Lundgren, P.R., Prats-Iraola, P. (2017): Compound dislocation models (CDMs) for volcano deformation analyses. Geophys. J. Int. https://doi.org/10.1093/gji/ggw427

---

Common interface

# Check if a point is inside the source volume
inside(Point3(x, y, z), src)

# Non-dimensionalize all parameters
src_nd = nondimensionalize(src, CharDim)

hostrock_displacement — computing displacement fields

For all source types, hostrock_displacement can be called on a single point or on arrays.

2D — single point

Ux, Uz = hostrock_displacement(src, Point2(x, z))

2D — array

X = range(-10e3, 10e3, length=128)
Z = range(-20e3,    0, length=128)
Ux, Uz = hostrock_displacement(src, X, Z)   # returns two arrays

3D — single point

Ux, Uy, Uz = hostrock_displacement(src, Point3(x, y, z))

3D — array

X = range(-10e3, 10e3, length=64)
Y = range(-10e3, 10e3, length=64)
Z = range(-20e3,    0, length=64)
Ux, Uy, Uz = hostrock_displacement(src, X, Y, Z)   # returns three arrays

For FiniteEllipsoidalCavity the array form additionally returns the volume change, potency, and number of point CDMs:

ue, un, uv, dV, DV, Ns = hostrock_displacement(fec, X, Y)

surface_displacement — displacement on a GeophysicalModelGenerator surface

When GeophysicalModelGenerator is loaded, surface_displacement computes the displacement at every point of a CartData surface. CartData uses km; the conversion to metres is handled internally.

using InjectSills, GeophysicalModelGenerator

# Build a CartData surface (flat, at sea level)
nx, ny, nz = 201, 201, 1
x2D = [xi for xi in range(-100.0, 100.0, length=nx), _ in 1:ny,  _ in 1:nz]  # km
y2D = [yi for _ in 1:nx, yi in range(-100.0, 100.0, length=ny),  _ in 1:nz]
z2D = zeros(nx, ny, nz)   # km; replace with topography if available
surf = CartData(x2D, y2D, z2D, (Elevation = z2D,))

# Define a source
src = MogiSphere(
    Center = Point3(0.0, 0.0, -5000.0)*m,
    r = 1500.0m, ΔP = 10e6Pa, G = 10e9Pa, ν = 0.25*NoUnits,
)

# Option 1 — return displacement arrays (metres), useful for inversion
Ux, Uy, Uz = surface_displacement(src, surf)

# Option 2 — add displacement as a vector field to the CartData
surf_out = surface_displacement(src, surf; add_fields = true)
# surf_out.fields.Displacement_m = (Ux, Uy, Uz) tuple in metres
# components: [1] = Ux, [2] = Uy, [3] = Uz

The add_fields = true form preserves all existing fields in surf and appends :Displacement_m as a (Ux, Uy, Uz) tuple — consistent with the GMG convention for vector fields. The plain array form is preferable for inversion workflows where you want to avoid constructing new CartData objects at each iteration.

A working example for the Galapagos is provided here, which looks like: