List of GeoParams functions

Here an overview of all functions in GeoParams.jl, for a complete list see here:

GeoParams.DislocationCreep_info Constant

SetDislocationCreep["Name of Dislocation Creep"]

Sets predefined dislocation creep data from a dictionary

source

GeoParams.ConstTemp Type

ConstTemp(T=1000)

Sets a constant temperature inside the box

Parameters:

• T : the value

source

GeoParams.HalfspaceCoolTemp Type

HalfspaceCoolTemp(Tsurface=0, Tmantle=1350, Age=60, Adiabat=0)

Sets a halfspace temperature structure in plate

Parameters:

• Tsurface: surface temperature [C]

• Tmantle: mantle temperature [C]

• Age: Thermal Age of plate [Myr]

• Adiabat: Mantle Adiabat [K/km]

source

GeoParams.LinTemp Type

LinTemp(Ttop=0, Tbot=1000)

Sets a linear temperature structure from top to bottom

Parameters:

• Ttop: the value @ the top

• Tbot: the value @ the bottom

source

GeoParams.CompTempStruct Method

CompTempStruct(Z, s::AbstractTempStruct)

Returns a temperature vector that matches the coordinate vector and temperature structure that were passed in

Parameters:

• Z: vector with depth coordinates

• s: Temperature structure (CostTemp, LinTemp, HalfspaceCoolTemp)

source

GeoParams.LithPres Method

LithPres(MatParam, Phases, ρ, T, dz, g)

Iteratively solves a 1D lithostatic pressure profile (compatible with temperature- and pressure-dependent densities)

Parameters:

• MatParam: a tuple of materials (including the following properties: Phase, Density)

• Phases: vector with the distribution of phases

• ρ: density vector for initial guess(can be zeros)

• T: temperature vector

• dz: grid spacing

• g: gravitational acceleration

source

GeoParams.StrengthEnvelopeComp Method

StrengthEnvelopeComp(MatParam::NTuple{N, AbstractMaterialParamsStruct}, Thickness::Vector{U}, TempType::AbstractTempStruct=LinTemp(0C, 800C), ε=1e-15/s, nz::Int64=101) where {N, U}

Calculates a 1D strength envelope. Pressure used for Drucker Prager plasticity is lithostatic. To visualize the results in a GUI, use StrengthEnvelopePlot.

Parameters:

• MatParam: a tuple of materials (including the following properties: Phase, Density, CreepLaws, Plasticity)

• Thickness: a vector listing the thicknesses of the respective layers (should carry units)

• TempType: the type of temperature profile (ConstTemp(), LinTemp(), HalfspaceCoolTemp())

• ε: background strainrate

• nz: optional argument controlling the number of points along the profile (default = 101)

Example:

julia> using GLMakie
julia> MatParam = (SetMaterialParams(Name="UC", Phase=1, Density=ConstantDensity(ρ=2700kg/m^3), CreepLaws = SetDislocationCreep("Wet Quartzite | Ueda et al. (2008)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)),
                   SetMaterialParams(Name="MC", Phase=2, Density=Density=ConstantDensity(ρ=2900kg/m^3), CreepLaws = SetDislocationCreep("Plagioclase An75 | Ji and Zhao (1993)"), Plasticity = DruckerPrager(ϕ=20.0, C=10MPa)),
                   SetMaterialParams(Name="LC", Phase=3, Density=PT_Density(ρ0=2900kg/m^3, α=3e-5/K, β=1e-10/Pa), CreepLaws = SetDislocationCreep("Maryland strong diabase | Mackwell et al. (1998)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)));
julia> Thickness = [15,10,15]*km;

julia> StrengthEnvelopeComp(MatParam, Thickness, LinTemp(), ε=1e-15/s)

source

GeoParams.StrengthEnvelopePlot Method

StrengthEnvelopePlot(MatParam, Thickness; TempType, nz)

Requires GLMakie

Creates a GUI that plots a 1D strength envelope. In the GUI, temperature profile and strain rate can be adjusted. The Drucker-Prager plasticity uses lithostatic pressure.

Parameters:

• MatParam: a tuple of materials (including the following properties: Phase, Density, CreepLaws, Plasticity)

• Thickness: a vector listing the thicknesses of the respective layers (should carry units)

• TempType: the type of temperature profile (LinTemp=default, HalfspaceCoolTemp, ConstTemp)

• nz: optional argument controlling the number of points along the profile (default = 101)

Example:

julia> using GLMakie
julia> MatParam = (SetMaterialParams(Name="UC", Phase=1, Density=ConstantDensity(ρ=2700kg/m^3), CreepLaws = SetDislocationCreep("Wet Quartzite | Ueda et al. (2008)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)),
                   SetMaterialParams(Name="MC", Phase=2, Density=Density=ConstantDensity(ρ=2900kg/m^3), CreepLaws = SetDislocationCreep("Plagioclase An75 | Ji and Zhao (1993)"), Plasticity = DruckerPrager(ϕ=20.0, C=10MPa)),
                   SetMaterialParams(Name="LC", Phase=3, Density=PT_Density(ρ0=2900kg/m^3, α=3e-5/K, β=1e-10/Pa), CreepLaws = SetDislocationCreep("Maryland strong diabase | Mackwell et al. (1998)"), Plasticity = DruckerPrager(ϕ=30.0, C=10MPa)));
julia> Thickness = [15,10,15]*km;

julia> StrengthEnvelopePlot(MatParam, Thickness, LinTemp())

source

GeoParams.compute_p_τij! Method

compute_p_τij!(Txx, Tyy, Txy, Tii, η_vep, P, Exx, Eyy, Exy, P_o, Txx_o, Tyy_o, Txy_o, phase, MatParam, dt)

Computes 2D deviatoric stress components (Txx,Tyy,Txy) and pressure P given deviatoric strainrate components (Exx,Eyy,Exy) and old deviatoric stresses (Txx_o, Tyy_o, Txy_o) and old pressure P_o (only used for viscoelastic cases). Also returned are Tii (second invariant of the deviatoric stress tensor), and η_vep the viscoelastoplastic effective viscosity. Also required as input is MatParam, the material parameters for every phase and phase, an integer array of size(Exx) that indicates the phase of every point.

This function assumes that strainrate points are collocated and that Exx,Eyy,Exy are at the same points.

source

GeoParams.compute_p_τij Method

P,τij, τII = compute_p_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple{N2,T}, P_old::T, args, τij_old::NTuple{3,T}, phase::I)

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a collocated grid case with a single phase phase.

source

GeoParams.compute_p_τij Method

τij, τII, η_eff = compute_p_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple, args, τij_old::NTuple, phases::NTuple)

This computes pressure p and deviatoric stress components τij, their second invariant τII, and effective viscosity η_eff for given deviatoric strainrates εij, old stresses τij_old, phases (integer) for every point and arguments args. It handles staggered grids of various taste

source

GeoParams.compute_p_τij Method

p,τij,τII = compute_p_τij(v, εij::NTuple{n,T}, P_old::T, args,  τij_old::NTuple{N,T})

Computes pressure p and deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a case that all points are collocated and we have a single phase.

source

GeoParams.compute_p_τij Method

p, τij, τII = compute_p_τij(v, εij::NTuple{N,Union{T,NTuple{4,T}}}, P_old::T, args, τij_old::NTuple{3,T})

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a staggered grid case with a single phase.

source

GeoParams.compute_p_τij_stagcenter! Method

compute_τij_stagcenter!(Txx, Tyy, Txy, Tii, η_vep,  P, Exx, Eyy, Exyv, P_o, Txx_o, Tyy_o, Txyv_o, phase_center, phase_vertex, MatParam, dt)

Updates deviatoric stresses on a staggered grid in a centered based manner (averages strainrates/old stresses from vertices -> centers). Take care of the sizes of the input matrixes:

• Txx,Tyy,Txy,Tii,P: 2D matrixes of size (nx,ny) which describe updated deviatoric stress components (x,y,xy, 2nd invariant) at the center point

• η_vep: viscoelastoplastic viscosity @ center (nx,ny)

• Exx,Eyy,P_o: deviatoric strain rate components & old pressure @ center

• Exy: shear strain rate @ vertices (nx+1,ny+1)

• Txx_o, Tyy_o: deviatoric stress normal components of last timestep at center (nx,ny)

• Txy_o: deviatoric stress shear components at vertices (nx+1,ny+1)

• phase_center: integer array with phase at center points (nx,ny)

• phase_vertex: integer array with phase at vertex points (nx+1,ny+1)

• MatParam: material parameter array that includes a CompositeRheology field which specifies the rheology for different phases

• dt: timestep

source

GeoParams.compute_viscosity_εII Method

compute_viscosity_εII(s::AbstractConstitutiveLaw, εII, kwargs...)

Compute effective viscosity given a 2nd invariant of the deviatoric strain rate tensor, extra parameters are passed as a named tuple, e.g., (;T=T)

source

GeoParams.compute_viscosity_τII Method

compute_viscosity_τII(s::AbstractConstitutiveLaw, τII, kwargs...)

Compute effective viscosity given a 2nd invariant of the deviatoric stress tensor and, extra parameters are passed as a named tuple, e.g., (;T=T)

source

GeoParams.compute_τij Method

τij, τII = compute_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple{N2,T}, args, τij_old::NTuple{3,T}, phase::I)

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a collocated grid case with a single phase phase.

source

GeoParams.compute_τij Method

τij, τII, η_eff = compute_τij(v::NTuple{N1,AbstractMaterialParamsStruct}, εij::NTuple, args, τij_old::NTuple, phases::NTuple)

This computes deviatoric stress components τij, their second invariant τII, and effective viscosity η_eff for given deviatoric strainrates εij, old stresses τij_old, phases (integer) for every point and arguments args. This handles various staggered grid arrangements; if staggered components are given as NTuple{4,T}, they will be averaged. Note that the phase of all staggered points should be the same.

source

GeoParams.compute_τij Method

τij,τII = compute_τij(v, εij::NTuple{n,T}, args, τij_old::NTuple{N,T})

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a case that all points are collocated and we have a single phase.

source

GeoParams.compute_τij Method

τij, τII = compute_τij(v, εij::NTuple{N,Union{T,NTuple{4,T}}}, args, τij_old::NTuple{3,T})

Computes deviatoric stress τij for given deviatoric strain rate εij, old stress τij_old, rheology v and arguments args. This is for a staggered grid case with a single phase.

source

GeoParams.compute_τij_stagcenter! Method

compute_τij_stagcenter!(Txx, Tyy, Txy, Tii, η_vep,  Exx, Eyy, Exyv, P, Txx_o, Tyy_o, Txyv_o, phase_center, phase_vertex, MatParam, dt)

Updates deviatoric stresses on a staggered grid in a centered based manner (averages strainrates/old stresses from vertices -> centers). Take care of the sizes of the input matrixes:

• Txx,Tyy,Txy,Tii: 2D matrixes of size (nx,ny) which describe updated deviatoric stress components (x,y,xy, 2nd invariant) at the center point

• η_vep: viscoelastoplastic viscosity @ center (nx,ny)

• Exx,Eyy,P: deviatoric strain rate components & pressure @ center

• Exy: shear strain rate @ vertices (nx+1,ny+1)

• Txx_o, Tyy_o: deviatoric stress normal components of last timestep at center (nx,ny)

• Txy_o: deviatoric stress shear components at vertices (nx+1,ny+1)

• phase_center: integer array with phase at center points (nx,ny)

• phase_vertex: integer array with phase at vertex points (nx+1,ny+1)

• MatParam: material parameter array that includes a CompositeRheology field which specifies the rheology for different phases

• dt: timestep

source

GeoParams.rotate_elastic_stress2D Method

rotate_elastic_stress2D(ω, τ::T, dt) where T

Bi-dimensional rotation of the elastic stress where τ is in the Voig notation and ω = 1/2(dux/dy - duy/dx)

source

GeoParams.second_invariant_staggered Method

second_invariant_staggered(Axx::NTuple{4,T}, Ayy::NTuple{4,T}, Azz::NTuple{4,T}, Aij::NTuple{3,T}) where {T}

Computes the second invariant of the 2D tensor A when its diagonal components need to be mapped from cell center to cell vertex. Axx, Ayy, and Azz are tuples containinig the diagonal terms of A at the cell centers around the i-th vertex., and Aij is a tuple that contains the off-diagonal components at the i-th vertex.

source

GeoParams.second_invariant_staggered Method

second_invariant_staggered(Axx::NTuple{4,T}, Ayy::NTuple{4,T}, Axy::Number) where {T}

Computes the second invariant of the 2D tensor A when its diagonal components need to be mapped from cell center to cell vertex. Axx, and Ayy are tuples containinig the diagonal terms of A at the cell centers around the i-th vertex., and Axy is the xy component at the i-th vertex.

source

GeoParams.second_invariant_staggered Method

second_invariant_staggered(Aii::NTuple{3,T}, Ayz::NTuple{4,T}, Axz::NTuple{4,T}, Axy::NTuple{4,T}) where {T}

Computes the second invariant of the 2D tensor A when its off-diagonal components need to be mapped from cell center to cell vertex. Aii is a tuple containinig the diagonal terms of A at the i-th vertex, and Ayz, Axz, and Axy are tuples that contain the off-diagonal components of the tensor at the cell centers around the i-th vertex.

source

GeoParams.second_invariant_staggered Method

second_invariant_staggered(Aii::NTuple{2,T}, Axy::NTuple{4,T}) where {T}

Computes the second invariant of the 2D tensor A when its off-diagonal components need to be mapped from cell center to cell vertex. Aii is a tuple containinig the diagonal terms of A at the i-th vertex, and Axy is a tuple that contains A_xy at the cell centers around the i-th vertex.

source

GeoParams.time_p_τII_0D! Method

time_p_τII_0D!(P_vec::Vector{T}, τ_vec::Vector{T}, v::CompositeRheology, εII_vec::Vector{T}, εvol_vec::Vector{T}, args, t_vec::AbstractVector{T}) where {T}

Computes stress-time evolution for a 0D (homogeneous) setup with given shear and volumetric strainrate vectors (which can vary with time).

source

GeoParams.time_p_τII_0D Method

t_vec, P_vec, τ_vec = time_p_τII_0D(v::CompositeRheology, εII::Number, εvol::Number, args; t=(0.,100.), τ0=0., nt::Int64=100)

This performs a 0D constant strainrate experiment for a composite rheology structure v, and a given, constant, shear strainrate εII and volumetric strainrate εvol, as well as rheology arguments args. The initial stress τ0, the time range t and the number of timesteps nt can be modified

source

GeoParams.time_τII_0D! Method

time_τII_0D!(τ_vec::Vector{NTuple{N,T}}, τII_vec::Vector{T}, v::Union{CompositeRheology,Tuple, Parallel}, ε_vec::Vector{NTuple{N,T}}, args, t_vec::AbstractVector{T}; verbose=false) where {N,T}

Computes stress-time evolution for a 0D (homogeneous) setup with given strainrate tensor (which can vary with time).

source

GeoParams.time_τII_0D! Method

time_τII_0D!(τ_vec::Vector{T}, v::CompositeRheology, εII_vec::Vector{T}, args, t_vec::AbstractVector{T}) where {T}

Computes stress-time evolution for a 0D (homogeneous) setup with given strainrate vector (which can vary with time).

source

GeoParams.time_τII_0D Method

t_vec, τ_vec = time_τII_0D(v::CompositeRheology, εII::Number, args; t=(0.,100.), τ0=0., nt::Int64=100)

This performs a 0D constant strainrate experiment for a composite rheology structure v, and a given, constant, strainrate εII and rheology arguments args. The initial stress τ0, the time range t and the number of timesteps nt can be modified

source

GeoParams.time_τII_0D Method

t_vec, τ_vec = time_τII_0D(v::CompositeRheology, ε::NTuple{N,T}, args; t=(0.,100.), τ0=NTuple{N,T}, nt::Int64=100)

This performs a 0D constant strainrate experiment for a composite rheology structure v, and a given, constant, strainrate tensor ε and rheology arguments args. The initial deviatoric stress tensor τ, the time range t and the number of timesteps nt can be modified

source

GeoParams.MaterialParameters.MaterialParams Type

MaterialParams

Structure that holds all material parameters for a given phase

source

GeoParams.MaterialParameters.MaterialParamsInfo Type

MaterialParamsInfo

Structure that holds information (Equation, Comment, BibTex_Reference) about a given material parameter, which can be used to create parameter tables, documentation etc.

Usually used in combination with param_info(the_parameter_of_interest)

source

GeoParams.MaterialParameters.SetMaterialParams Method

SetMaterialParams(; Name::String="", Phase::Int64=1,
                    Density             =   nothing,
                    Gravity             =   nothing,
                    CreepLaws           =   nothing,
                    Elasticity          =   nothing,
                    Plasticity          =   nothing,
                    CompositeRheology   =   nothing,
                    ChemDiffusion       =   nothing,
                    Conductivity        =   nothing,
                    HeatCapacity        =   nothing,
                    RadioactiveHeat     =   nothing,
                    LatentHeat          =   nothing,
                    ShearHeat           =   nothing,
                    Permeability        =   nothing,
                    Melting             =   nothing,
                    SeismicVelocity     =   nothing,
                    CharDim::GeoUnits   =   nothing)

Sets material parameters for a given phase.

If CharDim is specified the input parameters are non-dimensionalized. Note that if Density is specified, we also set Gravity even if not explicitly listed

Examples

Define two viscous creep laws & constant density:

julia> Phase = SetMaterialParams(Name="Viscous Matrix",
                       Density   = ConstantDensity(),
                       CreepLaws = (PowerlawViscous(), LinearViscous(η=1e21Pa*s)))
Phase 1 : Viscous Matrix
        | [dimensional units]
        |
        |-- Density           : Constant density: ρ=2900 kg m⁻³
        |-- Gravity           : Gravitational acceleration: g=9.81 m s⁻²
        |-- CreepLaws         : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹
        |                       Linear viscosity: η=1.0e21 Pa s

Define two viscous creep laws & P/T dependent density and nondimensionalize

julia> CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);
julia> Phase = SetMaterialParams(Name="Viscous Matrix", Phase=33,
                              Density   = PT_Density(),
                              CreepLaws = (PowerlawViscous(n=3), LinearViscous(η=1e23Pa*s)),
                              CharDim   = CharUnits_GEO)
Phase 33: Viscous Matrix
        | [non-dimensional units]
        |
        |-- Density           : P/T-dependent density: ρ0=2.9e-16, α=0.038194500000000006, β=0.01, T0=0.21454659702313156, P0=0.0
        |-- Gravity           : Gravitational acceleration: g=9.810000000000002e18
        |-- CreepLaws         : Powerlaw viscosity: η0=0.1, n=3, ε0=0.001
        |                       Linear viscosity: η=10000.0

You can also create an array that holds several parameters:

julia> MatParam        =   Array{MaterialParams, 1}(undef, 2);
julia> Phase           =   1;
julia> MatParam[Phase] =   SetMaterialParams(Name="Upper Crust", Phase=Phase,
                            CreepLaws= (PowerlawViscous(), LinearViscous(η=1e23Pa*s)),
                            Density  = ConstantDensity(ρ=2900kg/m^3));
julia> Phase           =   2;
julia> MatParam[Phase] =   SetMaterialParams(Name="Lower Crust", Phase=Phase,
                            CreepLaws= (PowerlawViscous(n=5), LinearViscous(η=1e21Pa*s)),
                            Density  = PT_Density(ρ0=3000kg/m^3));
julia> MatParam
2-element Vector{MaterialParams}:
 Phase 1 : Upper Crust
    | [dimensional units]
    |
    |-- Density           : Constant density: ρ=2900 kg m⁻³
    |-- Gravity           : Gravitational acceleration: g=9.81 m s⁻²
    |-- CreepLaws         : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹
    |                       Linear viscosity: η=1.0e23 Pa s

 Phase 2 : Lower Crust
    | [dimensional units]
    |
    |-- Density           : P/T-dependent density: ρ0=3000 kg m⁻³, α=3.0e-5 K⁻¹, β=1.0e-9 Pa⁻¹, T0=0 °C, P0=0 MPa
    |-- Gravity           : Gravitational acceleration: g=9.81 m s⁻²
    |-- CreepLaws         : Powerlaw viscosity: η0=1.0e18 Pa s, n=5, ε0=1.0e-15 s⁻¹
    |                       Linear viscosity: η=1.0e21 Pa s

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.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.ConstantDensity Type

ConstantDensity(ρ=2900kg/m^3)

Set a constant density:

$$\rho =cst$$

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

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.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.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

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.Vector_Density Type

Vector_Density(_T)

Stores a vector with density data that can be retrieved by providing an index

source

GeoParams.MaterialParameters.Density.Vector_Density Method

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

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

source

GeoParams.MaterialParameters.Density.compute_density! Method

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

GeoParams.MaterialParameters.Density.compute_density Method

compute_density(P,T, s::PhaseDiagram_LookupTable)

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

source

GeoParams.MaterialParameters.HeatCapacity.ConstantHeatCapacity Type

ConstantHeatCapacity(Cp=1050J/mol/kg)

Set a constant heat capacity:

$$Cp=cst$$

where $Cp$ is the thermal heat capacity [$J/kg/K$].

source

GeoParams.MaterialParameters.HeatCapacity.Latent_HeatCapacity Type

Latent_HeatCapacity(Cp=ConstantHeatCapacity(), Q_L=400kJ/kg)

This takes the effects of latent heat into account by modifying the heat capacity in the temperature equation:

$$\rho C_p^{\text{new}}\frac{\partial T}{\partial t}=\frac{\partial }{\partial x_i}\left(k\frac{\partial T}{\partial x_i}\right)+H_s$$

with

$$C_p^{\text{new}}=C_p+\frac{\partial \varphi }{\partial T}{Q}_L$$

where $Q_L$ is the latent heat [$J/kg$], and $\frac{\partial \varphi }{\partial T}$ is the derivative of the melt fraction with respect to temperature

source

GeoParams.MaterialParameters.HeatCapacity.T_HeatCapacity_Whittington Type

T_HeatCapacity_Whittington()

Sets a temperature-dependent heat capacity following the parameterization of Whittington et al. (2009), Nature:

$$Cp=(a+bT-c/T^2)/m$$

where $Cp$ is the heat capacity [$J/kg/K$], and $a,b,c$ are parameters that dependent on the temperature $T$ [$K$]:

• a = 199.50 J/mol/K if T<= 846 K

• a = 199.50 J/mol/K if T> 846 K

• b = 0.0857J/mol/K^2 if T<= 846 K

• b = 0.0323J/mol/K^2 if T> 846 K

• c = 5e6J/mol*K if T<= 846 K

• c = 47.9e-6J/mol*K if T> 846 K

• molmass = 0.22178kg/mol

Note that this is slightly different than the equation in the manuscript, as Cp is in J/kg/K (rather than $J/mol/K$ as in eq.3/4 of the paper)

source

GeoParams.MaterialParameters.HeatCapacity.Vector_HeatCapacity Type

Vector_HeatCapacity(_T)

Stores a vector with heat capacity data that can be retrieved by providing an index

source

GeoParams.MaterialParameters.HeatCapacity.compute_heatcapacity! Method

compute_heatcapacity!(Cp::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat})

In-place computation of heat capacity Cp 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.

source

GeoParams.MaterialParameters.HeatCapacity.compute_heatcapacity Method

compute_heatcapacity(a::Vector_HeatCapacity; index::Int64, kwargs...)

Pointwise calculation of heat capacity from a vector where index is the index of the point

source

GeoParams.MaterialParameters.HeatCapacity.compute_heatcapacity Method

Cp = compute_heatcapacity(s:<AbstractHeatCapacity, P, T)

Returns the heat capacity Cp at any temperature T and pressure P using any of the heat capacity laws implemented.

Currently available:

• ConstantHeatCapacity

• T_HeatCapacity_Whittington

• Latent_HeatCapacity

Example

Using dimensional units

julia> T = 250.0:100:1250
julia> Cp2 = T_HeatCapacity_Whittington()
julia> Cp = similar(T)
julia> args = (; T=T)
julia> Cp =compute_heatcapacity!(Cp, Cp2, args)
11-element Vector nitful.Quantity{Float64, 𝐋² 𝚯⁻¹ 𝐓⁻², Unitful.FreeUnits{(kg⁻¹, J, K⁻¹), 𝐋² 𝚯⁻¹ 𝐓⁻², nothing}}}:
635.4269997294616 J kg⁻¹ K⁻¹
850.7470171764261 J kg⁻¹ K⁻¹
962.0959598489883 J kg⁻¹ K⁻¹
1037.542043377064 J kg⁻¹ K⁻¹
1097.351792196648 J kg⁻¹ K⁻¹
1149.274556367170 J kg⁻¹ K⁻¹
1157.791505094840 J kg⁻¹ K⁻¹
1172.355487419726 J kg⁻¹ K⁻¹
1186.919469744596 J kg⁻¹ K⁻¹
1201.483452069455 J kg⁻¹ K⁻¹
1216.0474343943067 J kg⁻¹ K⁻¹

source

GeoParams.MaterialParameters.RadioactiveHeat.ConstantRadioactiveHeat Type

ConstantRadioactiveHeat(H_r=1e-6Watt/m^3)

Set a constant radioactive heat:

$$H_r=cst$$

where $H_r$ is the radioactive heat source [$Watt/m^3$].

source

GeoParams.MaterialParameters.RadioactiveHeat.ExpDepthDependentRadioactiveHeat Type

ExpDepthDependent(H_0=1e-6Watt/m^3, h_r=10e3m, z_0=0m)

Sets an exponential depth-dependent radioactive

$$H_r=H_0\exp (-\frac{z-z_0}{h_r})$$

where $H_0$ is the radioactive heat source [$Watt/m^3$] at $z=z_0$ which decays with depth over a characteristic distance $h_r$.

source

GeoParams.MaterialParameters.Shearheating.ConstantShearheating Type

ConstantShearheating(Χ=0.0NoUnits)

Set the shear heating efficiency [0-1] parameter

$${\rm X}=cst$$

where $\text{Chi}$ is the shear heating efficiency [NoUnits]

Shear heating is computed as

$$H_s=\text{Chi}\cdot {\tau }_{ij}({\dot{\epsilon }}_{ij}-{\dot{\epsilon }}_{ij}^{el})$$source

GeoParams.MaterialParameters.Shearheating.compute_shearheating! Method

compute_shearheating!(H_s, s:<AbstractShearheating,  τ, ε, ε_el)

Computes the shear heating source term in-place

$$H_s=\text{Chi}\cdot {\tau }_{ij}({\dot{\epsilon }}_{ij}-{\dot{\epsilon }}_{ij}^{el})$$

Parameters

• $\text{Chi}$ : The efficiency of shear heating (between 0-1)

• ${\tau }_{ij}$ : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]

• ${\dot{\epsilon }}_{ij}$ : The full deviatoric strainrate tensor

• ${\dot{\epsilon }}_{ij}^{el}$ : The full elastic deviatoric strainrate tensor

NOTE: The shear heating terms require the full deviatoric stress & strain rate tensors, i.e.:

$$2D:{\tau }_{ij}=\left(\begin{array}{cc}{\tau }_{xx}& {\tau }_{xz}\\ {\tau }_{zx}& {\tau }_{zz}\end{array}\right)$$

Since ${\tau }_{zx}={\tau }_{xz}$, most geodynamic codes only take one of the terms into account; shear heating requires all components to be used!

source

GeoParams.MaterialParameters.Shearheating.compute_shearheating! Method

compute_shearheating!(H_s, s:<AbstractShearheating, τ, ε)

Computes the shear heating source term H_s in-place when there is no elasticity

$$H_s=\text{Chi}\cdot {\tau }_{ij}{\dot{\epsilon }}_{ij}$$

Parameters

• $\text{Chi}$ : The efficiency of shear heating (between 0-1)

• ${\tau }_{ij}$ : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]

• ${\dot{\epsilon }}_{ij}$ : The full deviatoric strainrate tensor

source

GeoParams.MaterialParameters.Shearheating.compute_shearheating Method

H_s = compute_shearheating(s:<AbstractShearheating, τ, ε, ε_el)

Computes the shear heating source term

$$H_s=\text{Chi}\cdot {\tau }_{ij}({\dot{\epsilon }}_{ij}-{\dot{\epsilon }}_{ij}^{el})$$

Parameters

• $\text{Chi}$ : The efficiency of shear heating (between 0-1)

• ${\tau }_{ij}$ : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]

• ${\dot{\epsilon }}_{ij}$ : The full deviatoric strainrate tensor

• ${\dot{\epsilon }}_{ij}^{el}$ : The full elastic deviatoric strainrate tensor

source

GeoParams.MaterialParameters.Shearheating.compute_shearheating Method

H_s = ComputeShearheating(s:<AbstractShearheating, τ, ε)

Computes the shear heating source term when there is no elasticity

$$H_s=\text{Chi}\cdot {\tau }_{ij}{\dot{\epsilon }}_{ij}$$

Parameters

• $\text{Chi}$ : The efficiency of shear heating (between 0-1)

• ${\tau }_{ij}$ : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]

• ${\dot{\epsilon }}_{ij}$ : The full deviatoric strainrate tensor

source