Density

ConstantDensity(ρ=2900kg/m^3)

Set a constant density:

\[ \rho = cst\]

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

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.

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

Set a temperature-dependent density:

\[ \rho = \rho_0 (1 - \alpha * (T - T\_0) )\]

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

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

Set a pressure-dependent density:

\[ \rho = \rho_0 \exp(β*(P - P\_0))\]

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

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

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

\[ \rho = \phi \rho_{\textrm{melt}} + (1-\phi) \rho_{\textrm{solid}}\]

where $\rho$ is the average density [$kg/m^3$], $\rho_{\textrm{melt}}$ the melt density, $\rho_{\textrm{solid}}$ the solid density and $\phi$ the melt fraction.

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

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{cases} aP^{1/2} & \text{for } P < \frac{c_0^2}{a^2} \\ c_0 & \text{for } P \geq \frac{c_0^2}{a^2} \end{cases}\]

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

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

compute_density(P,T, s::PhaseDiagram_LookupTable)

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

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)