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

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

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

rheology = SetMaterialParams(;
                      Phase=1,
                      CreepLaws=(PowerlawViscous(), LinearViscous(; η=1e21Pa * s)),
                      Gravity=ConstantGravity(; g=9.81.0m / s^2),
                      Density= GasPyroclast_Density(ρmelt=ConstantDensity(ρ=2900kg/m^3), ρgas=ConstantDensity(ρ=1kg/m^3), δ=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

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

Computational routines

To evaluate density within a user routine, use this:

GeoParams.MaterialParameters.Density.compute_density Function

compute_density(P,T, s::PhaseDiagram_LookupTable)

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(η=1e23Pas)),
                        Density   = ConstantDensity(ρ=2900kg/m^3))
                  );
julia> Phases = ones(Int64,10,10);
julia> Phases[:,5:end] .= 2
julia> rho     = zeros(size(Phases))
julia> T       =  ones(size(Phases))
julia> P       =  ones(size(Phases))*10
julia> args = (P=P, T=T)
julia> compute_density!(rho, MatParam, Phases, args)
julia> rho
10×10 Matrix{Float64}:
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0
 2899.91  2899.91  2899.91  2899.91  2900.0  2900.0  2900.0  2900.0  2900.0  2900.0

The routine is made to minimize allocations:

julia> using BenchmarkTools
julia> @btime compute_density!($rho, $MatParam, $Phases, P=$P, T=$T)
    203.468 μs (0 allocations: 0 bytes)

---

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.