Köpping et al., 2023 - Intrusion tip velocities control emplacement mechanisms of sheet intrusions
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
## change font for plots
import matplotlib
from matplotlib import rc
matplotlib.rcParams['pdf.fonttype'] = 42 # to import text in illustrator
matplotlib.rcParams['ps.fonttype'] = 42
rc('font',**{'family':'serif','serif':['Arial']})
SMALL_SIZE = 8
MEDIUM_SIZE = 10
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
## define cm2inch for plot size
def cm2inch(*tupl): # cm to inch for figure size
inch = 2.54
if isinstance(tupl[0], tuple):
return tuple(i/inch for i in tupl[0])
else:
return tuple(i/inch for i in tupl)
Below, we approximate the temperature change ahead of a propagating intrusion tip, due to heat transfer by thermal diffusivity.
Unknown parameter
$T$ = Temperature (ºC) ahead of a propagating intrusion tip at location $L$
Known parameters
$T_{m}$ = Magma temperature (ºC); 1000–1200 ºC (Cas and Wright, 1987)
$T_{\infty}$ = Background temperature (ºC); 52.5 ºC at 1500 m depth, assuming a thermal gradient of 2.5 ºC per 100 m and a surface temperature of 15 ºC
$U$ = Intrusion tip velocity $(m s^{-1})$
$x$ = Length (m) ahead of a propagating intrusion tip
$L_{d}$ = Characteristic length (m) ahead of a propagating intrusion tip that is affected by heat diffusion
$L_{adv}$ = Distance (m) travelled by the intrusion tip moving at velocity $U$ after time $t$
$\kappa$ = Thermal diffusivity $(m^{2} s^{-1})$
An average $\kappa$ for a porous host rock including pore fluids can be estimated by:
$\kappa = \phi * \kappa_{water} + (1-\phi) * \kappa_{rock}$,
with $\phi$ = porosity, ranging from 0 to 1.
Heat from the intruding magma diffuses into host rocks at rates higher than tip propagation when $L_{d}>L{adv}$. In this scenario, pore fluids ahead of the propagating tip may be heated to high enough temperatures for flash boiling to occur.
A temperature change over time in one dimension can be described by:
$\frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2}$, with $\frac{\partial T}{\partial t} = -U \frac{\partial T}{\partial x}$ (Turcotte and Schubert, 2002; Eq. 4.68)
Substitution leads to a heat equation in the frame moving with the propagating tip velocity ($U$), where $x$ equals the distance ahead of the intrusion tip:
$-U \frac{\partial T}{\partial x} = \kappa \frac{\partial^2 T}{\partial x^2}$ (Eq. 1)
Eq. 1 can be solved for $T$ using integration, which results in:
$T = \hat{\hat{A}} + Be^{-\frac{U}{\kappa}x}$, with $\hat{\hat{A}} = T_{\infty}$, and $B=(T_m-T_{\infty})$.
The temperature ahead of an igneous intrusion propagating at velocity ($U$) is therefore described by:
$T = T_{\infty} + (T_{m}-T_{\infty})e^{-\frac{U}{\kappa}x}$ (Eq. 2)
# thermal diffusivity
K_w = 0.143E-6 # (m^2/s) for water at 25 deg C
K_r = 1.3E-6 # (m^2/s); 1.3e-6 sandstone, 0.1e-6 coal, 0.8e-6 shale
phi = 0.15 # porosity
K = phi * K_w + (1-phi) * K_r # (m^2/s)
# background temperature; surface temp 15 deg C + 2.5 deg C per 100m depth
T_inf = 52.5 # 52.5 deg C for 1500 m depth
# intrusion temperature
T_int = [1000, 1200] # (deg C), mafic intrusion 1000–1200; (Cas & Wright, 1988)
# intrusion velocity
U = [1e-6, 1e-5, 1e-4] # (m/s)
# distance ahead of intrusion
step = 1e-3
L = np.arange(0, 2 + step, step)
c = ['gainsboro', 'grey', 'dimgrey']
ls = ['dashed', 'dotted', 'solid']
label = [r'$10^{-6}$ $m$ $s^{-1}$',
r'$10^{-5}$ $m$ $s^{-1}$',
r'$10^{-4}$ $m$ $s^{-1}$']
fig, ax = plt.subplots(figsize=cm2inch(16,12))
for i in range(len(U)): # loop over tip velocity range
# calculate temperature at position X
# for upper and lower temperature limit (1000–1200 deg C) using Eq. 2
Tx_lower = T_inf + (T_int[0]-T_inf)*np.exp(-(U[i]/K)*L)
Tx_upper = T_inf + (T_int[1]-T_inf)*np.exp(-(U[i]/K)*L)
ax.fill_between(L,Tx_lower, Tx_upper, fc=c[i],
ec='k', ls=ls[i], linewidth = 1, label = label[i])
ax.set(xlabel='x (m)')
ax.set(ylabel='T (ºC)')
ax.hlines(350, 0, 2, linestyles='dashed', linewidth=1, color='r')
ax.hlines(100, 0, 2, linestyles='dashed', linewidth=1, color='r')
ax.set_xlim([0, 2])
ax.set_ylim([0,1220])
ax.legend(loc='upper right')
fig.tight_layout()
plt.show()
Assuming sheet intrusions propagate via a linear elastic hydraulic fracture, their tip velocities (V) may be approximated by aperture ($w$) measurements in the near tip region.
The aperture ($w$) of a radial fracture is described by:
$w \sim 2 \cdot 3^{\frac{7}{6}} \cdot (\frac{\mu V}{E'})^{\frac{1}{3}} \cdot s^{\frac{2}{3}}$ (Eq. 3) Desroches et al. (1994), Xing et al. (2017)
with $w$ = fracture opening, $s$ = distance from the fracture tip, $\mu$ = fluid viscosity, $E'$= plain strain modulus, $V$ = fracture propagation velocity.
By inverting Eq. 3 for unknown $V$, the tip velocity of a radial fracture is described by:
$V \sim \frac{E'w^3}{216 \cdot 3^\frac{1}{2} \mu s^2}$ (Eq. 4)
The aperture ($w$) at a given distance from the crack tip ($s$) is constrained based on field examples of dykes and sills as presented in Galland et al. (2018) (Fig. 5.6), Poppe et al. (2020) (Figs. 3E, 5A, 5C), Schmiedel et al. (2021) (Fig. 4), and Walker et al. (2021) (Figs. 1C, 1D, 1F, 1G).
We note that fracture tip velocities may be influenced by local changes in host rock properties (e.g., fracture toughness and/or elastic moduli), which have not been considered in our model.
Galland et al., (2018), Storage and Transport of Magma in the Layered Crust—Formation of Sills and Related Flat-Lying Intrusions. Volcanic and Igneous Plumbing Systems, 113-138.
Poppe et al., (2020), Structural and Geochemical Interactions Between Magma and Sedimentary Host Rock: The Hovedøya Case, Oslo Rift, Norway. Geochemistry, Geophysics, Geosystems, 21, 1-22.
Schmiedel et al., (2021), Emplacement and Segment Geometry of Large, High-Viscosity Magmatic Sheets. Minerals, 11, 1113.
Walker et al., (2021), Segment tip geometry of sheet intrusions, I: Theory and numerical models for the role of tip shape in controlling propagation pathways. Volcanica, 4, 189-201.
Equations 3 and 4 are established in:
Desroches et al., (1994), The crack tip region in hydraulic fracturing. Proceedings of the Royal Society of London, 447, 39-48.
Xing et al., (2017), Laboratory measurement of tip and global behaviour for zero-toughness hydraulic fractures with circular and blade-shaped (PKN) geometry. Journal of the Mechanics and Physics of Solids, 104, 172-186.
## Import t and s measurements as collected from field photographs
# import data (t and s measurements from literature)
df_input = pd.read_excel('/PATH/2_Supplemental_material_s-t-measurements_literature.xlsx')
# create empty dataframe
df_V = pd.DataFrame(columns=('ID', 'Host_rock', 'mu', 'E', 'v', 'E_p', 's', 'w',
'V_rad', 'KIC', 'ell_rad'))
# create array with range of viscosities (10^1 – 10^15 Pa s)
mu1 = np.array(())
for i in range (1,10):
temp = i * np.logspace(1, 14, 14)
mu1 = np.append(mu1, temp)
mu1 = np.append(mu1, 1e15)
mu1 = np.sort(mu1)
# copy entries of viscosity array multiple times for all t values
mu = np.array(())
for i in range (len(df_input)):
mu = np.append(mu, mu1)
## write ID, Host rock, mu, s, t to dataframe
for i in range (len(df_input)):
for j in range (len(mu1)):
df_temp = pd.DataFrame(data={'ID':[df_input['ID'][i]],
'Host_rock':[df_input['Host_rock'][i]],
'mu':[mu1[j]],
'E':[np.nan],
'v':[np.nan],
'E_p':[np.nan],
's':[df_input['s'][i]],
'w':[df_input['t'][i]],
'V_rad':[np.nan],
'KIC':[np.nan],
'ell_rad':[np.nan]})
df_V = pd.concat([df_V, df_temp], ignore_index=True)
## Define Young's modulus (E) and Poisson's ration (v) for the different host rocks
# Data from: Gudmundsson (2011), Rock Fractures in Geological Processes
df_V.loc[df_V['Host_rock']=='Mudstone', ('E')] = 15e9 # Małkowski et al., 2018
df_V.loc[df_V['Host_rock']=='Gabbro', ('E')] = 70e9
df_V.loc[df_V['Host_rock']=='Sandstone', ('E')] = 35e9
df_V.loc[df_V['Host_rock']=='Granite', ('E')] = 35e9
df_V.loc[df_V['Host_rock']=='Lava_flows', ('E')] = 5.4e9 # Heap et al., 2020
df_V.loc[df_V['Host_rock']=='Mudstone', ('v')] = 0.2
df_V.loc[df_V['Host_rock']=='Gabbro', ('v')] = 0.2
df_V.loc[df_V['Host_rock']=='Sandstone', ('v')] = 0.2
df_V.loc[df_V['Host_rock']=='Granite', ('v')] = 0.2
df_V.loc[df_V['Host_rock']=='Lava_flows', ('v')] = 0.2
# Determine plane strain modulus (E')
df_V['E_p'] = df_V['E']/(1 - df_V['v']**2)
# Fracture propagation velocity (Eq. 4)
df_V['V_rad'] = ((df_V['E_p'] * df_V['w']**3) /
(216 * np.sqrt(3) * df_V['mu'] * df_V['s']**2))
## Figure
alpha = 0.4 # transparency
lw = 1 # line width
# axis limits
x_min = 1e4
x_max = 1e14
y_max = 1e0
y_min = 1e-8
unique_ID = df_V['ID'].unique()
fig, ax = plt.subplots(figsize=cm2inch(13.5,10.5))
## Scatter plot
# ax.scatter(df_V['mu'], df_V['V_rad'], marker='.', s=20, c='k', edgecolor='none', alpha=alpha)
## Line plot
for i in range(len(unique_ID)):
df_ = df_V[df_V['ID']==unique_ID[i]]
df_.reset_index(inplace=True)
ax.plot(df_['mu'], df_['V_rad'], lw=lw, c='k', alpha=alpha, zorder=1)
ax.set(xlabel='$\mu$ $(Pa$ $s)$')
ax.set(ylabel='$V$ $(m$ $s^{-1})$')
ax.set_xscale('log')
ax.set_yscale('log')
# axis grid
ax.grid(visible=True, which='both', axis='both',
color='grey', linestyle='--', linewidth=1, alpha=0.3)
# horizontal line at V=10^-5 m s^-1
ax.hlines(1e-5, x_min, x_max, linestyles='dashed', linewidth=0.5, color='r')
ax.set_xlim([x_min, x_max])
ax.set_ylim([y_min,y_max])
fig.tight_layout()
plt.show()
Here we check if the model of viscous limited propagation holds at the velocity and viscosity combinations used above.
For this, we assume that there is a region nearest to the fracture tip that is governed by LEFM, and then farther from the tip it could be governed by viscous flow. If the size of the region governed by LEFM (referred to as tip characteristic length "$\ell$") is smaller than the distance "$s$" used in the models above, fracture growth was governed by viscous flow such that the modelled results may be interpreted as meaningful.
To check this, we calculate the tip characteristic length ($\ell$):
$\ell = \frac{K_{IC}^6}{E'^4 \mu^2 V^2}$ Eq. 5 Bunger & Detournay (2008)
Bunger & Detournay (2008), Experimental validation of the tip asymptotics for a fluid-driven crack, Journal of the Mechanics and Physics of Solids, 56, 3101–3115.
## Determine ell
df_V['KIC'] = 1e7 # in the range of 1e5 – 1e7
df_V['ell_rad'] = (( df_V['KIC']**6) /
((df_V['E_p']**4) * df_V['mu']**2 * df_V['V_rad']**2))
len(df_V[df_V['s']<df_V['ell_rad']])
0
For $K_{IC}$ values $\leq 10^7 Pa$ $m^{-2}$, $\ell < s$, suggesting that fracture growth was governed by viscous flow. Our modelled data presented above may thus be considered as meaningful.