Density

Methods

The density equation of state can be specified in a number of ways

GeoParams.MaterialParameters.Density.ConstantDensity Type

ConstantDensity(ρ=2900kg/m^3)

Set a constant density:

$$\rho =cst$$

where $\rho$ is the density [$kg/m^3$].

source

GeoParams.MaterialParameters.Density.PT_Density Type

PT_Density(ρ0=2900kg/m^3, α=3e-5/K, β=1e-9/Pa, T0=0C, P=0MPa)

Set a pressure and temperature-dependent density:

$$\rho ={\rho }_0(1.0-\alpha (T-T_0)+\beta (P-P_0))$$

where ${\rho }_0$ is the density [$kg/m^3$] at reference temperature $T_0$ and pressure $P_0$, $\alpha$ is the temperature dependence of density and $\beta$ the pressure dependence.

source

GeoParams.MaterialParameters.Density.T_Density Type

T_Density(ρ0=2900kg/m^3, α=3e-5/K, T₀=273.15K)

Set a temperature-dependent density:

$$\rho ={\rho }_0(1-\alpha \ast (T-T\mathrm{\_}0))$$

where ${\rho }_0$ is the density [$kg/m^3$] at reference temperature $T_0$ and $\alpha$ the temperature dependence.

source

GeoParams.MaterialParameters.Density.Compressible_Density Type

Compressible_Density(ρ0=2900kg/m^3, β=1e-9/Pa, P₀=0MPa)

Set a pressure-dependent density:

$$\rho ={\rho }_0\exp (\beta \ast (P-P\mathrm{\_}0))$$

where ${\rho }_0$ is the density [$kg/m^3$] at reference pressure $P_0$ and $\beta$ the pressure dependence.

source

GeoParams.MaterialParameters.Density.MeltDependent_Density Type

MeltDependent_Density(ρsolid=ConstantDensity(), ρmelt=ConstantDensity())

If we use a single phase code the average density of a partially molten rock is

$$\rho =\varphi {\rho }_{\text{melt}}+(1-\varphi ){\rho }_{\text{solid}}$$

where $\rho$ is the average density [$kg/m^3$], ${\rho }_{\text{melt}}$ the melt density, ${\rho }_{\text{solid}}$ the solid density and $\varphi$ the melt fraction.

Note that any density formulation can be used for melt and solid.

source

GeoParams.MaterialParameters.Density.BubbleFlow_Density Type

BubbleFlow_Density(ρmelt=ConstantDensity(), ρgas=ConstantDensity(), c0=0e0, a=0.0041MPa^(-1//2))

Defines the BubbleFlow_Density as described in Slezin (2003) with a default gas solubility constant of 0.0041MPa${}^{-1/2}$ used in e.g. Sparks et al. (1978)

$$\rho =\frac{1}{\frac{c_0-c}{{\rho }_g}+\frac{1-(c_0-c)}{{\rho }_m}}$$

with

$$c=\{\begin{array}{ll}a{P}^{1/2}& \text{for }P<\frac{c_0^2}{a^2}\\ c_0& \text{for }P\ge \frac{c_0^2}{a^2}\end{array}$$

Arguments

• ρmelt: Density of the melt

• ρgas: Density of the gas

• c0: Total volatile content

• a: Gas solubility constant (default: 4.1e-6Pa${}^{-1/2}$) (after Sparks et al., 1978)

Possible values for a are 3.2e-6-6.4e-6Pa${}^{-1/2}$ where the lower value corresponds to mafic magmas at rather large pressures (400-600MPa) and the higher value to felsic magmas at low pressures (0 to 100-200MPa) (after Slezin (2003))

Example

julia> rheology = SetMaterialParams(;
           Phase = 1,
           CreepLaws = (PowerlawViscous(), LinearViscous(; η = 1e21Pa * s)),
           Gravity = ConstantGravity(; g = 9.81m / s^2),
           Density = BubbleFlow_Density(;
               ρmelt = ConstantDensity(ρ = 2900kg / m^3),
               ρgas = ConstantDensity(ρ = 1kg / m^3),
               c0 = 0.0,
               a = 0.0041MPa^(-1 // 2),
           ),
       );

julia> rheology.Density[1]
Bubble flow density: ρ = 1/((c0-c)/ρgas + (1-(c0-c))/ρmelt); ρmelt=Constant density: ρ=2900.0 kg m⁻³·⁰; ρgas=Constant density: ρ=1.0 kg m⁻³·⁰; c0=0.0; a=0.0041 MPa⁻⁰·⁵

References

• Slezin, Yu. B. (2003), The mechanism of volcanic eruptions (a steady state approach), Journal of Volcanology and Geothermal Research, 122, 7-50, https://doi.org/10.1016/S0377-0273(02)00464-X

• Sparks, R. S. J.(1978), The dynamics of bubble formation and growth in magmas: A review and analysis, Journal of Volcanology and Geothermal Research, 3, 1-37, https://doi.org/10.1016/0377-0273(78)90002-1

source

GeoParams.MaterialParameters.Density.GasPyroclast_Density Type

GasPyroclast_Density(ρmelt=ConstantDensity(), ρgas=ConstantDensity(), δ=0e0)

Defines the GasPyroclast_Density as described in Slezin (2003) with a default volume fraction of free gas in the flow of 0.0 This is also used to model partly destroyed foam in the conduit.

$$\rho ={\rho }_g\delta +{\rho }_p(1-\delta )$$

with

$${\rho }_p={\rho }_m(1-\beta )+{\rho }_g\beta \approx {\rho }_l(1-\beta )$$

Arguments

• ρmelt: Density of the melt

• ρgas: Density of the gas

• δ: Volume fraction of free gas in the flow

• β: Gas volume fraction enclosed within the particles

Example

julia> rheology = SetMaterialParams(;
           Phase = 1,
           CreepLaws = (PowerlawViscous(), LinearViscous(; η = 1e21Pa * s)),
           Gravity = ConstantGravity(; g = 9.81m / s^2),
           Density = GasPyroclast_Density(;
               ρmelt = ConstantDensity(ρ = 2900kg / m^3),
               ρgas = ConstantDensity(ρ = 1kg / m^3),
               δ = 0.0,
               β = 0.0,
           ),
       );

julia> rheology.Density[1]
Gas-Pyroclast mixture density: ρ = ρgas*δ + ρmelt*(1-β); ρmelt=Constant density: ρ=2900.0 kg m⁻³·⁰; ρgas=Constant density: ρ=1.0 kg m⁻³·⁰; δ=0.0; β=0.0

References

• Slezin, Yu. B. (2003), The mechanism of volcanic eruptions (a steady state approach), Journal of Volcanology and Geothermal Research, 122, 7-50, https://doi.org/10.1016/S0377-0273(02)00464-X

source

GeoParams.MaterialParameters.Density.Melt_DensityX Type

Melt_DensityX(; oxd_wt = oxd_wt)

Set a density depending on the oxide composition after the python script by Iacovino K & Till C (2019)

Arguments

• oxd_wt::NTuple{9, T}: Melt composition as 9-element Tuple containing concentrations in [wt%] of the following oxides ordered in the exact sequence

         (SiO2 TiO2 Al2O3 FeO MgO CaO Na2O K2O H2O) 
  
         Default values are for a hydrous N-MORB composition

The callable also accepts an mH2O keyword (melt water content, mass fraction) that overrides the struct's own oxd_wt[9] for that call, recomputing only the water-dependent aggregates. mH2O is accepted as given: GeoParams does not check it for consistency with any Solubility or other closure's dissolved-water output — that is the caller's/solver's responsibility.

References

• Iacovino K & Till C (2019). DensityX: A program for calculating the densities of magmatic liquids up to 1,627 °C and 30 kbar. Volcanica 2(1), p 1-10. doi:10.30909/vol.02.01.0110

source

GeoParams.MaterialParameters.Density.RedlichKwong_Density Type

RedlichKwong_Density(; coeffs=(-112.528, 127.811, 112.04), T0=273.15K, Tref=1K, Pref=1e5Pa, ρref=1e3kg/m^3)

Modified Redlich–Kwong gas density of Huber et al. (2010), as used by Degruyter & Huber (2014). Fitted for H2O gas over roughly $873<T<1173$ K and $30<P<400$ MPa:

$${\rho }_g={\rho }_{ref}(a_1{\tau }^{-0.381}+a_2{\varpi }^{-1.135}+a_3{\tau }^{-0.411}{\varpi }^{0.033}),\tau =\frac{T-T_0}{T_{ref}},\varpi =\frac{P}{P_{ref}}$$

The fit's specific units (°C, bar) enter through the reference quantities: with T0=273.15K, Tref=1K the group $\tau$ equals $T$ in °C, and with Pref=1e5Pa (1 bar) $\varpi$ equals $P$ in bar. Because $\tau$ and $\varpi$ are ratios of like-dimensioned quantities, the expression is dimensionally homogeneous and nondimensionalizes cleanly: only the reference GeoUnits scale, the fitted dimensionless coeffs do not. The parameterisation ignores gas composition (pseudo-pure H2O), matching the reference.

Returns NaN for T ≤ 273.15K or P ≤ 0 (τ ≤ 0 or ω ≤ 0), since τ^(-0.381)/τ^(-0.411) and ω^(-1.135) require a positive base: raising to those exponents would otherwise throw an opaque low-level DomainError about complex exponentiation.

References

• Huber C., Bachmann O. , Manga M., Two Competing Effects of Volatiles on Heat Transfer in Crystal-rich Magmas: Thermal Insulation vs Defrosting, Journal of Petrology, 847–867, https://doi.org/10.1093/petrology/egq003

• 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.MaterialParameters.Density.IdealGas_Density Type

IdealGas_Density(; Rs=461.5J/kg/K)

Ideal-gas density

$${\rho }_g=\frac{P}{R_{s}T}$$

with specific gas constant Rs (default: water vapor, $R/M_w=8.314/0.01802$). Unlike RedlichKwong_Density this is dimensionally homogeneous, so it nondimensionalizes cleanly and serves as the analytic derivative check for the Redlich–Kwong fit.

source

GeoParams.MaterialParameters.Density.ThreePhase_Density Type

ThreePhase_Density(; ρmelt=ConstantDensity(ρ=2300kg/m^3), ρsolid=ConstantDensity(ρ=2600kg/m^3), ρgas=RedlichKwong_Density(), ρ=2400kg/m^3)

Volume-fraction mixture density of a melt + crystal + gas magma (Degruyter & Huber 2014):

$$\rho ={\epsilon }_m{\rho }_m+{\epsilon }_g{\rho }_g+{\epsilon }_x{\rho }_x,{\epsilon }_m=1-{\epsilon }_x-{\epsilon }_g$$

The gas and crystal volume fractions ϕ_gas, ϕ_x are inputs (from the solver's ODE state / a melting law), passed through args; the sub-densities receive the same args (e.g. P, T) so an equation-of-state gas density is evaluated in place. ρ is a dimensional-tracking sentinel.

ρgas is only evaluated when ϕ_gas != 0: a cell with no exsolved gas never depends on the gas closure at all, so a strict EOS like RedlichKwong_Density (which returns NaN outside its calibration window) cannot poison a ϕ_gas=0 cell regardless of P, T.

Arguments

• ρmelt: melt-phase density

• ρsolid: crystal-phase density

• ρgas: gas-phase density (e.g. RedlichKwong_Density)

References

• 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

Computational routines

To evaluate density within a user routine, use this:

GeoParams.MaterialParameters.Density.compute_density Function

compute_density(s::Vector_Density; index::Int64, kwargs...)

Pointwise calculation of density from a vector where index is the index of the point

source

compute_density(P,T, s::AbstractPhaseDiagramsStruct)

Interpolates density as a function of T,P from a lookup table

source

GeoParams.MaterialParameters.Density.compute_density! Function

compute_density!(rho::AbstractArray{_T, ndim}, MatParam::NTuple{N,AbstractMaterialParamsStruct}, Phases::AbstractArray{_I, ndim}; P=nothing, T=nothing) where {ndim,N,_T,_I<:Integer}

In-place computation of density rho for the whole domain and all phases, in case a vector with phase properties MatParam is provided, along with P and T arrays. This assumes that the Phase of every point is specified as an Integer in the Phases array.

Example

julia> MatParam = (
           SetMaterialParams(Name = "Mantle", Phase = 1,
               CreepLaws = (PowerlawViscous(), LinearViscous(η = 1e23Pa * s)),
               Density = PT_Density()),
           SetMaterialParams(Name = "Crust", Phase = 2,
               CreepLaws = (PowerlawViscous(), LinearViscous(η = 1e23Pa * s)),
               Density = ConstantDensity(ρ = 2900kg / m^3)),
       );

julia> Phases = ones(Int64, 2, 4); Phases[:, 3:end] .= 2;

julia> rho = zeros(size(Phases));

julia> args = (P = ones(size(Phases)) * 10, T = ones(size(Phases)));

julia> compute_density!(rho, MatParam, Phases, args)

julia> rho
2×4 Matrix{Float64}:
 2899.91  2899.91  2900.0  2900.0
 2899.91  2899.91  2900.0  2900.0

The routine is made to minimize allocations:

using BenchmarkTools
@btime compute_density!($rho, $MatParam, $Phases, P = $P, T = $T)

---

compute_density!(rho::AbstractArray{_T, N}, MatParam::NTuple{K,AbstractMaterialParamsStruct}, PhaseRatios::AbstractArray{_T, M}, P=nothing, T=nothing)

In-place computation of density rho for the whole domain and all phases, in case a vector with phase properties MatParam is provided, along with P and T arrays. This assumes that the PhaseRatio of every point is specified as an Integer in the PhaseRatios array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)

source

Note that density values are usually not used in itself in the governing PDE's, but usually in combination with other parameters, such as $\rho g$ or $\rho c_p$. The non-dimensional value of $\rho$ may thus have very large or small values, but multiplied with the other values one often obtains numbers that are closer to one.