List of GeoParams functions

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

GeoParams.GeoParams Module

module GeoParams

Typical geodynamic simulations involve a large number of material parameters that have units that are often inconvenient to be directly used in numerical models This package has two main features that help with this:

• Create a nondimensionalization object, which can be used to transfer dimensional to non-dimensional parameters (usually better for numerical solvers)

• Create an object in which you can specify material parameters employed in the geodynamic simulations

The material parameter object is designed to be extensible and can be passed on to the solvers, such that new creep laws or features can be readily added. We also implement some typically used creep law parameters, together with tools to plot them versus and compare our results with those of published papers (to minimize mistakes).

source

GeoParams.DislocationCreep_info Constant

SetDislocationCreep["Name of Dislocation Creep"]

Sets predefined dislocation creep data from a dictionary

source

GeoParams.AbstractComposite Type

AbstractComposite <: AbstractMaterialParam

Supertype of composite rheologies that combine several constitutive laws (see CompositeRheology and Parallel).

source

GeoParams.AbstractConstitutiveLaw Type

AbstractConstitutiveLaw{T} <: AbstractMaterialParam

Supertype of the constitutive laws that relate stress and strain rate (creep laws, elasticity, plasticity). T is the numeric element type.

source

GeoParams.AbstractMaterialParam Type

AbstractMaterialParam

Supertype of every individual material property parameterization (density, elasticity, viscosity, conductivity, …). Concrete subtypes are callable and evaluate the property they describe.

source

GeoParams.AbstractMaterialParamsStruct Type

AbstractMaterialParamsStruct

Supertype of the per-phase container that bundles all AbstractMaterialParams belonging to a single material phase (see MaterialParams).

source

GeoParams.AbstractPhaseDiagramsStruct Type

AbstractPhaseDiagramsStruct <: AbstractMaterialParam

Supertype of material parameters obtained by interpolating a precomputed phase diagram lookup table as a function of pressure and temperature.

source

GeoParams.ConstTemp Type

ConstTemp(T=1000)

Sets a constant temperature inside the box

Parameters:

• T : the value

source

GeoParams.ConstantDensityTrait Type

DensityTrait marking a constant density; returned by isconstant.

source

GeoParams.DensityTrait Type

DensityTrait

Supertype of the Holy traits classifying a density parameterization as either ConstantDensityTrait or NonConstantDensityTrait, returned by isconstant.

source

GeoParams.ElasticRheologyTrait Type

RheologyTrait marking a rheology that contains elasticity; returned by isviscoelastic.

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

RheologyTrait marking a linear rheology; returned by islinear.

source

GeoParams.NonConstantDensityTrait Type

DensityTrait marking a pressure/temperature-dependent density; returned by isconstant.

source

GeoParams.NonElasticRheologyTrait Type

RheologyTrait marking a rheology with no elasticity; returned by isviscoelastic.

source

GeoParams.NonLinearRheologyTrait Type

RheologyTrait marking a non-linear rheology; returned by islinear.

source

GeoParams.NonPlasticRheologyTrait Type

RheologyTrait marking a rheology with no plasticity; returned by isplasticity.

source

GeoParams.PlasticRheologyTrait Type

RheologyTrait marking a rheology that contains plasticity; returned by isplasticity.

source

GeoParams.RheologyTrait Type

RheologyTrait

Supertype of the Holy traits used to classify a rheology. The concrete singletons come in mutually exclusive pairs — LinearRheologyTrait/NonLinearRheologyTrait, ElasticRheologyTrait/NonElasticRheologyTrait, PlasticRheologyTrait/NonPlasticRheologyTrait — and are returned by islinear, isviscoelastic, and isplasticity to enable branch-free dispatch on a rheology's properties.

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

PlotConductivity(k::AbstractConductivity; kwargs...)

Plots thermal conductivity as a function of temperature for the parameterization k.

source

GeoParams.PlotDeformationMap Function

PlotDeformationMap(v; kwargs...)

Plots a deformation-mechanism map (deformation regime as a function of temperature and stress or strain rate) for the rheology v.

source

GeoParams.PlotDiffusionCoefArrhenius Function

PlotDiffusionCoefArrhenius(x; kwargs...)

Arrhenius plot of the diffusion coefficient (log(D) versus 10⁴/T) for one or more ChemicalDiffusionData structures x.

source

GeoParams.PlotHeatCapacity Function

PlotHeatCapacity(Cp::AbstractHeatCapacity; kwargs...)

Plots heat capacity as a function of temperature for the parameterization Cp.

source

GeoParams.PlotMeltFraction Function

PlotMeltFraction(p::AbstractMeltingParam; kwargs...)

Plots melt fraction and dϕ/dT as a function of temperature for the parameterization p.

source

GeoParams.PlotPhaseDiagram Function

PlotPhaseDiagram(p::AbstractPhaseDiagramsStruct, fieldname::Symbol; kwargs...)

Plots the field fieldname of a phase diagram as a function of temperature (x-axis) and pressure (y-axis).

source

GeoParams.PlotPressureStressTime_0D Function

PlotPressureStressTime_0D(x; εII, εvol, kwargs...)

Plots the pressure and stress evolution over time of a 0-D visco-elasto-(plastic) model for the rheology x.

source

GeoParams.PlotStrainrateStress Function

PlotStrainrateStress(x; kwargs...)

Plots deviatoric stress versus deviatoric strain rate for one or more creep laws x.

source

GeoParams.PlotStrainrateViscosity Function

PlotStrainrateViscosity(x; kwargs...)

Plots effective viscosity versus deviatoric strain rate for one or more creep laws x.

source

GeoParams.PlotStressStrainrate Function

PlotStressStrainrate(x; kwargs...)

Plots deviatoric strain rate versus deviatoric stress for one or more creep laws x (the transpose of PlotStrainrateStress).

source

GeoParams.PlotStressTime_0D Function

PlotStressTime_0D(x; εII, kwargs...)

Plots the stress evolution over time of a 0-D visco-elasto-(plastic) model for the rheology x.

source

GeoParams.PlotStressViscosity Function

PlotStressViscosity(x; kwargs...)

Plots effective viscosity versus deviatoric stress for one or more creep laws x.

source

GeoParams.Plot_TAS_diagram Function

Plot_TAS_diagram(point; kwargs...)

Plots a TAS (total-alkali versus silica) classification diagram for the given composition point.

source

GeoParams.Plot_ZirconAge_PDF Function

Plot_ZirconAge_PDF(time_Ma, PDF_zircons, time_Ma_average, PDF_zircon_average)

Plots the zircon-age probability density function computed from a simulation.

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_elasticviscosity Method

compute_elasticviscosity(v, args)

Returns the viscosity contributed by the elastic elements of the rheology v alone (G·dt), combining them harmonically for a composite rheology.

source

GeoParams.compute_elastoviscosity Method

compute_elastoviscosity(v, η, dt)
compute_elastoviscosity(G, η, dt)

Returns the effective visco-elastic viscosity that combines a viscous viscosity η and elasticity (shear modulus G, or a ConstantElasticity v) over a time step dt, as the harmonic sum (1/η + 1/(G·dt))⁻¹. The time step may instead be passed via args.dt.

source

GeoParams.compute_elastoviscosity_εII Method

compute_elastoviscosity_εII(v, εII, args)

Returns the effective visco-elastic viscosity of the rheology v evaluated from the deviatoric strain rate εII, combining the non-plastic elements harmonically.

source

GeoParams.compute_elastoviscosity_τII Method

compute_elastoviscosity_τII(v, τII, args)

Returns the effective visco-elastic viscosity of the rheology v evaluated from the deviatoric stress τII, combining the non-plastic elements harmonically.

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 Method

compute_viscosity(v, args)

Returns the effective viscosity of the linear rheology v (a LinearViscous, ConstantElasticity, CompositeRheology, or MaterialParams), combining non-plastic elements harmonically. For multiple phases, v may be a tuple of MaterialParams indexed by a phase or weighted by phase ratios. Throws for non-linear rheologies.

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

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

Computes 2D deviatoric stress components (Txx,Tyy,Txy) given deviatoric strainrate components (Exx,Eyy,Exy) and old deviatoric stresses (Txx_o, Tyy_o, Txy_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_τ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.diffusion_law_list Method

Returns the list of pre-defined diffusion creep laws (entries of the Diffusion submodule).

source

GeoParams.dislocation_law_list Method

Returns the list of pre-defined dislocation creep laws (entries of the Dislocation submodule).

source

GeoParams.get_G Function

get_G(r)

Returns the elastic shear modulus G of the rheology r (a ConstantElasticity, CompositeRheology, or MaterialParams), summing the contributions of elastic elements and returning 0 when r is non-elastic. Also available as get_shearmodulus.

source

GeoParams.get_Kb Function

get_Kb(r)

Returns the elastic bulk modulus Kb of the rheology r (a ConstantElasticity, CompositeRheology, or MaterialParams), summing the contributions of elastic elements and returning 0 when r is non-elastic. Also available as get_bulkmodulus.

source

GeoParams.grainboundarysliding_law_list Method

Returns the list of pre-defined grain-boundary-sliding creep laws (entries of the GBS submodule).

source

GeoParams.isconstant Method

isconstant(v) -> DensityTrait

Returns ConstantDensityTrait if the density parameterization v (or MaterialParams, or tuple thereof) is pressure- and temperature-independent, and NonConstantDensityTrait otherwise.

source

GeoParams.islinear Method

islinear(v) -> RheologyTrait

Returns LinearRheologyTrait if the rheology v (a constitutive law, CompositeRheology, MaterialParams, or tuple thereof) is linear, and NonLinearRheologyTrait otherwise. When given two arguments, returns the linear trait only if both are linear.

source

GeoParams.isplasticity Method

isplasticity(v) -> RheologyTrait

Returns PlasticRheologyTrait if the rheology v contains a plastic element, and NonPlasticRheologyTrait otherwise.

source

GeoParams.isviscoelastic Method

isviscoelastic(v) -> RheologyTrait

Returns ElasticRheologyTrait if the rheology v contains an elastic element, and NonElasticRheologyTrait otherwise. When given two arguments, returns the elastic trait if either is elastic.

source

GeoParams.nonlinearpeierls_law_list Method

Returns the list of pre-defined non-linear Peierls creep laws (entries of the NonLinearPeierls submodule).

source

GeoParams.ntuple_idx Method

ntuple_idx(args::NamedTuple, I...)

Indexes each array-valued field of the NamedTuple args at the indices I, returning a new NamedTuple with the same keys and the scalar values at that position.

source

GeoParams.param_info Function

param_info(s::AbstractMaterialParam) -> MaterialParamsInfo

Returns a MaterialParamsInfo describing the parameterization s: its governing equation (as a LaTeXString) and, where available, a comment and BibTeX reference. Each concrete material-parameter type provides its own method.

source

GeoParams.peierls_law_list Method

Returns the list of pre-defined Peierls creep laws (entries of the Peierls submodule).

source

GeoParams.rotate_elastic_stress Method

rotate_elastic_stress(ω, τ, dt)

Applies the Jaumann co-rotation to the elastic stress tensor τ (in Voigt notation) over a time step dt, given the vorticity ω. Dispatches to the 2D or 3D rotation depending on the length of ω.

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 Method

second_invariant(A)
second_invariant(xx, yy, xy)
second_invariant(xx, yy, zz, yz, xz, xy)

Computes the second invariant $\sqrt{\frac{1}{2}{A}_{ij}{A}_{ij}}$ of a deviatoric tensor A. The tensor may be given as an NTuple, SVector, SMatrix, or Matrix, or by its independent components (Voigt order) for the 2D (xx, yy, xy) or 3D (xx, yy, zz, yz, xz, xy) case; the component forms are differentiable.

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

value_and_partial(f, x::Real) -> (f(x), f'(x))

Evaluates the scalar function f and its first derivative at x in a single forward-mode automatic-differentiation pass using ForwardDiff.

source

GeoParams.@use Macro

@use GeoParamsAliases property = alias ...

Creates a GeoParamsAliases module in the calling scope defining short-hand alias functions for the compute_* routines. Each keyword maps one of the supported properties (e.g. density, conductivity, heatcapacity, …) to the name you want to call it by. For example, @use GeoParamsAliases density = ρ makes ρ an alias for compute_density.

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

No_MaterialParam{_T}()

Placeholder AbstractMaterialParam representing the absence of a given property, used to pad the per-phase parameter fields so that every phase has the same field layout.

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,
                    Solubility          =   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.AbstractDensity Type

AbstractDensity{T} <: AbstractMaterialParam

Supertype of density parameterizations (e.g. ConstantDensity, PT_Density).

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

ConduitDensity{T} <: AbstractDensity{T}

Supertype of density parameterizations for volcanic conduit flow.

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

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

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.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::AbstractPhaseDiagramsStruct)

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

source

GeoParams.MaterialParameters.Density.compute_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_ratio Method

compute_density_ratio(phase_ratios, MatParam, args)

Returns the phase-fraction-weighted average density at a point, summing each phase's density (compute_density) times its fraction in phase_ratios.

source

GeoParams.MaterialParameters.HeatCapacity.ConstantHeatCapacity Type

ConstantHeatCapacity(Cp=1050J/kg/K)

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 = 229.32 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(P,T, s::AbstractPhaseDiagramsStruct)

Interpolates heat capacity as a function of T,P from a lookup table

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.RadioactiveHeat.compute_radioactive_heat! Method

compute_radioactive_heat!(H_r, s, args)

In-place version of compute_radioactive_heat that fills the array H_r for the whole domain.

source

GeoParams.MaterialParameters.RadioactiveHeat.compute_radioactive_heat Method

compute_radioactive_heat(s::AbstractRadioactiveHeat, args)

Returns the radioactive heat production for the parameterization s (which may depend on depth z).

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