Melting Parameterizations

Methods

A number of melting parameterisations are implemented, which can be set with:

GeoParams.MeltingParam.MeltingParam_Caricchi Type

MeltingParam_Caricchi()

Implements the T-dependent melting parameterisation used by Caricchi, Simpson et al. (as for example described in Simpson)

$$\theta =\frac{a-(T+c)}{b}$$$${\varphi }_{melt}=\frac{1.0}{1.0+e^{\theta }}$$

Note that T is in Kelvin. As default parameters we employ:

$$b=23K,a=800K,c=273.15K$$

Which gives a reasonable fit to experimental data of granodioritic composition (Piwinskii and Wyllie, 1968):

References

• Simpson G. (2017) Practical finite element modelling in Earth Sciences Using MATLAB.

source

GeoParams.MeltingParam.MeltingParam_Smooth3rdOrder Type

MeltingParam_Smooth3rdOrder()

Implements the a smooth 3rd order T-dependent melting parameterisation (as used by Melnik and coworkers)

$$x=\frac{T-273.15}{1000.0}$$$$\theta =a+b\ast x+c\ast x^2+d\ast x^3$$$${\varphi }_{melt}=\frac{1.0}{1.0+e^{\theta }}$$

Note that T is in Kelvin.

As default parameters we employ:

$$a=517.9,b=-1619.0,c=1699.0,d=-597.4$$

which gives a reasonable fit to experimental data for basalt.

Data for rhyolite are:

$$a=3043.0,b=-10552.0,c=12204.9,d=-4709.0$$

Red: Rhyolite, Blue: Basalt

References

source

GeoParams.MeltingParam.MeltingParam_5thOrder Type

MeltingParam_5thOrder(a,b,c,d,e,f,T_s,T_l)

Uses a 5th order polynomial to describe the melt fraction phi between solidus temperature T_s and liquidus temperature T_l

$$\varphi =a{T}^5+b{T}^4+c{T}^3+d{T}^2+eT+f\text{ for }{T}_s\le T\le T_l$$$$\varphi =1\text{ if }T>T_l$$$$\varphi =0\text{ if }T<T_s$$

Temperature T is in Kelvin.

The default values are for a composite liquid-line-of-descent:

• the upper part is for Andesite from: (Blatter, D. L. & Carmichael, I. S. (2001) Hydrous phase equilibria of a Mexican highsilica andesite: a candidate for a mantle origin? Geochim. Cosmochim. Acta 65, 4043–4065

• the lower part is extrapolated to the granitic minimum using the Marxer & Ulmer LLD for Andesite (Marxer, F. & Ulmer, P. (2019) Crystallisation and zircon saturation of calc-alkaline tonalite from the Adamello Batholith at upper crustal conditions: an experimental study. Contributions Mineral. Petrol. 174, 84)

source

GeoParams.MeltingParam.MeltingParam_4thOrder Type

MeltingParam_4thOrder(b,c,d,e,f,T_s,T_l)

Uses a 4th order polynomial to describe the melt fraction phi between solidus temperature T_s and liquidus temperature T_l

$$\varphi =b{T}^4+c{T}^3+d{T}^2+eT+f\text{ for }{T}_s\le T\le T_l$$$$\varphi =1\text{ if }T>T_l$$$$\varphi =0\text{ if }T<T_s$$

Temperature T is in Kelvin.

The default values are for Tonalite experiments from Marxer and Ulmer (2019):

• Marxer, F. & Ulmer, P. (2019) Crystallisation and zircon saturation of calc-alkaline tonalite from the Adamello Batholith at upper crustal conditions: an experimental study. Contributions Mineral. Petrol. 174, 84

source

GeoParams.MeltingParam.MeltingParam_Quadratic Type

MeltingParam_Quadratic(T_s,T_l)

Quadratic melt fraction parameterisation where melt fraction $\varphi$ depends only on solidus ($T_s$) and liquidus ($T_l$) temperature:

$$\varphi =1.0-{\left(\frac{T_l-T}{T_l-T_s}\right)}^2$$$$\varphi =1.0\text{ if }T>T_l$$$$\varphi =0.0\text{ if }T<T_s$$

Temperature T is in Kelvin.

This was used, among others, in Tierney et al. (2016) Geology

source

GeoParams.MeltingParam.MeltingParam_Assimilation Type

MeltingParam_Assimilation(T_s,T_l,a)

Melt fraction parameterisation that takes the assimilation of crustal host rocks into account, as used by Tierney et al. (2016) based upon a parameterisation of Spera and Bohrson (2001)

Here, the fraction of molten and assimilated host rocks $\varphi$ depends on the solidus ($T_s$) and liquidus ($T_l$) temperatures of the rocks, as well as on a parameter $a=0.005$

$$X=\frac{T-T_s}{T_l-T_s}$$$$\varphi =a\cdot ({\exp }^{2ln(100)X}-1.0)\text{ if }X\le 0.5$$$$\varphi =1-a\cdot {\exp }^{2ln(100)(1-X)}\text{ if }X>0.5$$$$\varphi =1.0\text{ if }T>T_l$$$$\varphi =0.0\text{ if }T<T_s$$

Temperature T is in Kelvin.

This was used, among others, in Tierney et al. (2016), who employed as default parameters:

$$T_s=973.15,T_l=1173.15,a=0.005$$

References

• Spera, F.J., and Bohrson, W.A., 2001, Energy-Constrained Open-System Magmatic Processes I: General Model and Energy-Constrained Assimilation and Fractional Crystallization (EC- AFC) Formulation: Journal of Petrology, v. 42, p. 999–1018.

• Tierney, C.R., Schmitt, A.K., Lovera, O.M., de Silva, S.L., 2016. Voluminous plutonism during volcanic quiescence revealed by thermochemical modeling of zircon. Geology 44, 683–686. https://doi.org/10.1130/G37968.1

source

GeoParams.MeltingParam.MeltingParam_Volatile Type

MeltingParam_Volatile(; a_coeffs, b_coeffs, c_coeffs, T0=273K, Tref=1K, Pref=1e6Pa)

Volatile-dependent silicic crystallinity-temperature curve fitted by Scholz et al. (2023) to rhyolite-MELTS runs, for use in Degruyter & Huber (2014)-type magma-chamber box models. The crystal fraction follows a complementary error function whose amplitude a, width b, and centre c depend on the dissolved water and CO2 contents and pressure:

$${\epsilon }_x=a\mathrm{erfc}(b(T_C-c)),\varphi =1-{\epsilon }_x$$

with $T_C=(T-T_0)/T_{ref}$ (numerically °C) and each of $a,b,c$ a quadratic polynomial in $x=100m_{H_{2}O}$, $y=100m_{C{O}_2}$ (dissolved mass fractions, in wt%) and $z=P/P_{ref}$ (numerically MPa). Dissolved water depresses the liquidus, so more water shifts the curve to lower temperature. The returned quantity is the melt fraction $\varphi$; the exsolved-gas fraction ${\epsilon }_g$ is a solver state and is not subtracted here.

The fit is calibrated only over crystal fractions from 0 up to ~0.5-0.6 (the "eruptible" range these box models track); 2a is the crystal fraction at the cold end of the fit and is bounded to [0.5, 1] by construction, so ϕ has a positive floor 1 - 2a rather than reaching 0 — this is intentional, not a missing clamp.

The dimensionless polynomial coefficients are stored as plain NTuples; only the reference scales T0, Tref, Pref are GeoUnits, so the parameterisation nondimensionalizes. Dissolved contents mH2O, mCO2 are dimensionless and pass through args.

References

• Scholz, K., Townsend, M., Huber, C., Troch, J., Bachmann, O., Coonin, A.N. (2023), Investigating the impact of an exsolved H2O-CO2 phase on magma chamber growth and longevity: A thermomechanical model, G-cubed 24, e2023GC011151, https://doi.org/10.1029/2023GC011151 (Eq. 9-12, silicic crystallinity-temperature curve)

• Degruyter, W., Huber, C. (2014), A model for the eruption frequency of upper crustal silicic magma chambers, EPSL 403, 117-130, https://doi.org/10.1016/j.epsl.2014.06.047

source

GeoParams.MeltingParam.MeltingParam_MaficVolatile Type

MeltingParam_MaficVolatile(; a_coeffs, b_coeffs, T0=273K, Tref=1K, Pref=1e6Pa)

Volatile-dependent mafic crystallinity-temperature curve fitted by Scholz et al. (2023) to rhyolite-MELTS runs, for use in Degruyter & Huber (2014)-type magma-chamber box models: linear in temperature, ε_x = a*T_C + b, with slope a and intercept b each a degree-2 polynomial in (100*mH2O, 100*mCO2, P/Pref) — the same polynomial form as MeltingParam_Volatile's a/b/c, just combined linearly instead of through erfc. ϕ = 1 - ε_x (ε_g is added by the solver, not here).

A linear model has no saturation, so ε_x leaves [0,1] outside a narrow temperature window (~100-200 K wide; at 200 MPa: 1337-1443 K dry, 1221-1436 K at 1 wt% H2O, 1010-1417 K at 2 wt%). compute_meltfraction clamps ε_x to [0,1] there instead of returning it as computed:

• colder than the window, the fit predicts ${\epsilon }_x>1$, clamped to $1$ ($\varphi =0$, fully solid);

• hotter than the window, it predicts ${\epsilon }_x<0$, clamped to $0$ ($\varphi =1$, fully liquid).

Both clamps are the physically sensible reading, but the value is still an extrapolation past the fit's calibration, not a measurement — no warning is raised when this triggers.

References

• Scholz, K., Townsend, M., Huber, C., Troch, J., Bachmann, O., Coonin, A.N. (2023), Investigating the impact of an exsolved H2O-CO2 phase on magma chamber growth and longevity: A thermomechanical model, G-cubed 24, e2023GC011151, https://doi.org/10.1029/2023GC011151 (Eq. 14-16, mafic crystallinity-temperature curve)

• Degruyter, W., Huber, C. (2014), A model for the eruption frequency of upper crustal silicic magma chambers, EPSL 403, 117-130, https://doi.org/10.1016/j.epsl.2014.06.047

source

GeoParams.MeltingParam.SmoothMelting Type

SmoothMelting(; p=MeltingParam_4thOrder(), k_sol=0.2/K,  k_liq=0.2/K)

This smoothens the melting parameterisation $p$ around the solidus $T_{sol}$ and liquidus $T_{liq}$ using a smoothened Heaviside step functions for the solidus:

$$H_{sol}=\frac{1.0}{1+\exp (-2k_{sol}(T-T_{sol}-\frac{2}{k_{sol}}))}$$

and liquidus:

$$H_{liq}=1.0-\frac{1.0}{1+\exp (-2k_{liq}(T-T_{liq}+\frac{2}{k_{liq}}))}$$

The resulting melt fraction $\varphi$ is computed from the original melt fraction ${\varphi }_0$ (computed using one of the methods above) as:

$$\varphi ={\varphi }_0{H}_{sol}{H}_{liq}+1.0-H_{liq}$$

The width of the smoothening zones is controlled by $k_{sol},k_{liq}$ (larger values = sharper boundary).

This is important, as jumps in the derivative $d\varphi /dT$ can cause numerical instabilities in latent heat computations, which is prevented with this smoothening.

Example

Smoothening a 4th order parameterisation, with a sharper liquidus than solidus:

julia> p_s = SmoothMelting(p = MeltingParam_4thOrder(), k_liq = 0.21/K)
4th order polynomial melting curve: phi = -7.594512597174117e-10T^4 + 3.469192091489447e-6T^3 + -0.00592352980926T^2 + 4.482855645604745T + -1268.730161921053  963.15 K ≤ T ≤ 1270.15 K with smooth Heaviside function smoothening using k_sol=0.2 K⁻¹·⁰, k_liq=0.21 K⁻¹·⁰

With a Makie backend loaded, PlotMeltFraction returns the curve and its derivative for both the original and the smoothened parameterisation:

using GLMakie, GeoParams
p = MeltingParam_4thOrder()
T = collect(650.0:1:1050.0) .+ 273.15
T, phi, dϕdT = PlotMeltFraction(p, T = T)
T_s, phi_s, dϕdT_s = PlotMeltFraction(p_s, T = T)

The derivative no longer has a jump now:

source

Computational routines

To compute the melt fraction at given T and P, use:

GeoParams.MeltingParam.compute_meltfraction! Function

compute_meltfraction!(ϕ::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T:AbstractArray{<:AbstractFloat}, p::AbstractPhaseDiagramsStruct)

In-place computation of melt fraction in case we use a phase diagram lookup table. The table should have the column :meltFrac specified.

source

compute_meltfraction(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties MatParam is provided, along with P and T arrays.

source

GeoParams.MeltingParam.compute_meltfraction Function

compute_meltfraction(P,T, p::AbstractPhaseDiagramsStruct)

Computes melt fraction in case we use a phase diagram lookup table. The table should have the column :meltFrac specified.

source

ϕ = compute_meltfraction(Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

Computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties MatParam is provided, along with P and T arrays.

source

You can also obtain the derivative of melt fraction versus temperature with (useful to compute latent heat effects):

GeoParams.MeltingParam.compute_dϕdT! Function

compute_dϕdT!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

Computes the derivative of melt fraction ϕ versus temperature T, \frac{\partial \phi}{\partial T} for the whole domain and all phases, in case an array with phase properties MatParam is provided, along with P and T arrays. This is employed, for example, in computing latent heat terms in an implicit manner.

source

GeoParams.MeltingParam.compute_dϕdT Function

compute_dϕdT(P,T, p::AbstractPhaseDiagramsStruct)

Computes derivative of melt fraction vs T in case we use a phase diagram lookup table. The table should have the column :meltFrac specified. The derivative is computed by finite differencing.

source

ϕ = compute_dϕdT(Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

Computates the derivative of melt fraction ϕ versus temperature T for the whole domain and all phases, in case an array with phase properties MatParam is provided, along with P and T arrays. This is employed in computing latent heat terms in an implicit manner, for example

source

Also note that phase diagrams can be imported using PerpleX_LaMEM_Diagram, which may also have melt content information. The computational routines work with that as well.

Plotting routines

You can use the routine PlotMeltFraction to create a plot, provided that a Makie backend has been loaded.