Equations governing movement of partial molten rock#
Following McKenzie (1984), the two-phase flow of a partially molten rock can be described by the following equations:
Mass conservation
Solid phase:
Melt phase:
where, \(\phi\) is the porosity (or melt fraction), \(\rho_s\) and \(\rho_f\) are the densities of the solid and melt phases, \(\mathbf{u}_s\) and \(\mathbf{u}_f\) are the velocities of the matrix and melt phases, and \(\Gamma\) is the rate at which melt is transferred from matrix to melt.
Momentum conservation
The terms on the left-hand side represent the deviatoric stresses (\(\tau_s\)) in the solid matrix, the compaction pressure due to compaction viscosity (\(\xi\)) of the matrix as the melt is removed, and the darcy drag force experienced by the melt of viscosity \(\eta_f\) as it flows through a porous media with permeability \(k\), respectively. The right-hand side represents the buoyancy force due to the density difference between the solid and melt phases.
Parameterization of melt production#
The melt fraction, \(\phi\), can be parameterized as a function of the temperature, \(T\), pressure, \(P\), residual modal cpx, \(M_{cpx}\), and weight fraction of water dissolved in melt, \(X_{H2O}\) for a peridotite mantle composition following Katz et al. (2003).
Based on the Bowen’s series, cpx melts first followed by opx. Degree of melting, \(F_{cpx}\), is given by:
where, liquidus of lherzolite and solidus temperatures include pressure dependence through empirical relationship, \(T_{liquidus}^{lzerz}(P) = B_1 + B_2P + B_3P^2\) and \(T_{solidus}(P) = A_1 + A_2P + A_3P^2\), respectively.
For \( F > F_{cpx}\), melting reaction changes such that mostly opx is consumed and the degree of melting is given by:
where, \(T_{cpx-out}\) is the temperature at which cpx from peridotite is completely consumed into the melt, and can be calculated from equation 1 assuming a modal cpx content in the peridotite composition.
Plotting the melt fraction#
import numpy as np
import matplotlib.pyplot as plt
#-----parameterized constants
# Constants for expressing liquidus temperature of the model peridotite system
C1 = 1780 # degree C
C2 = 45 # degree C GPa^-1
C3 = -2.0 # degree C GPa^-2
# Constants for expressing solidus temperature of the model peridotite system
A1 = 1085.7 # degree C
A2 = 132.9 # degree C GPa^-1
A3 = -5.1 # degree C GPa^-2
# Constants for expressing liquidus temperature of lherzolite
B1 = 1475 # degree C
B2 = 80 # degree C GPa^-1
B3 = -3.2 # degree C GPa^-2
# Constants related to the melt fraction
beta1 = 1.50
beta2 = 1.50
# Constants related to reaction coefficient of cpx
r0 = 0.5 # cpx/melt
r1 = 0.08 # cpx/melt/GPa
def melting_opx(wt_fraction_cpx, temperatures, pressure):
'''
Calculate the amount of melt fraction in a peridotite system following Katz (2002)
parameterization of a modal peridotite. The equations are described in detail in
section 2.1 of Katz et al., (2002)
Parameters:
- modal_cpx: float
Modal abundance of clinopyroxene in the residue (weight fraction of melt).
- temperature: numpy 1-D array
Temperature in degree K.
- pressure: numpy 1-D array
Pressure in GPa.
Returns:
- melt_fraction: numpy 1-D array
Melt fraction in weight percent.
'''
temperature_lherz_liq = B1 + B2*pressure + B3*pressure**2
temperature_solidus = A1 + A2*pressure + A3*pressure**2
temperature_liquidus = C1 + C2*pressure + C3*pressure**2
modal_cpx = wt_fraction_cpx/(r0 + r1 * pressure)
temperature_cpx_out = modal_cpx**(1/beta1)*(temperature_lherz_liq - temperature_solidus) + temperature_solidus
melt_fraction = np.zeros_like(temperatures)
for i in range(len(temperatures)):
if (temperatures[i] <= temperature_cpx_out) and (temperatures[i] >= temperature_solidus):
melt_fraction[i] = ((temperatures[i] - temperature_solidus) / (temperature_lherz_liq - temperature_solidus))**beta1
elif (temperatures[i] > temperature_cpx_out):
melt_fraction[i] = modal_cpx + (1 - modal_cpx) * ((temperatures[i] - temperature_cpx_out) / (temperature_liquidus - temperature_cpx_out))**beta2
return melt_fraction
# plot the melting curve as a function of temperature and pressure
pressures = np.array([0, 1, 2, 3]) # GPa
temperatures = np.linspace(1000, 2000, 50) # degree C
cpx_fraction = 0.15 # weight fraction of cpx in rock
fig, ax = plt.subplots(figsize=(5, 5))
for pressure in pressures:
melt_fraction = melting_opx(cpx_fraction, temperatures, pressure)
ax.plot(temperatures, melt_fraction, label=f'{pressure} GPa')
ax.set_xlabel('Temperature (°C)')
ax.set_ylabel('Melt Fraction')
ax.set_ylim([0, 1])
ax.legend()
The melting curves above are at a constant pressure specified in the legend. The kink in the plot occurs at the temperature where cpx is completely consumed. The melt fraction increases nonlinearly with increasing temperature and pressure, as expected.
References#
McKenzie, D. A. N. “The generation and compaction of partially molten rock.” Journal of petrology 25.3 (1984): 713-765.
Katz, R. F., Spiegelman, M., & Langmuir, C. H. (2003). A new parameterization of hydrous mantle melting. Geochemistry, Geophysics, Geosystems, 4(9)