Chamber.jl

IC_Finder(composition::Mafic, M_h2o::Float64, M_co2::Float64, M_tot::Float64, P::Float64, T::Float64, V::Float64, rho_m::Float64, param_IC::ParamICFinder{Float64})::NamedTuple{(:eps_g0, :X_co20, :mco2_diss, :phase),NTuple{4,Float64}}

An iterative function that uses thermodynamic modeling to determine the gas phase composition of a mafic magma chamber.

Arguments

• composition: Mafic()
• M_h2o: total mass of water in magma (kg)
• M_co2: total mass of CO2 in magma (kg)
• M_tot: total mass of magma (kg)
• P: Pressure (Pa)
• T: Temperature (K)
• V: chamber volume (m³)
• rho_m: density of the melt (kg/m³)
• param_IC: An instance of the ParamICFinder composite type containing the parameters for IC_Finder

Returns

A NamedTuple with the following fields:

• eps_g0: Initial gas fraction of the magma chamber
• X_co20: Initial mole fraction of CO2 in the gas phase
• mco2_diss: Total mass of CO2 dissolved in the melt (kg)
• phase: Integer representing the state of the magma

Details

The function uses a thermodynamic model to find the initial gas fraction (eps_g0) and initial mole fraction of CO2 in the gas phase (X_co20) for a mafic magma chamber. The model assumes that the chamber contains two phases: a liquid melt and a gas phase. The model calculates the saturation state of the magma with respect to CO2, and determines whether the gas phase contains only CO2 or a mixture of CO2 and H2O. If the chamber is in two-phase state, the model returns eps_g0 = 0.0, X_co20 = 0.0, mco2_diss = Total mass of CO2 dissolved in the melt, and phase = 2. If the chamber is in gas or liquid state, the function iteratively solves for eps_g0 and X_co20 using the solve_X_co2 and get_eps_g helper functions until the relative difference between eps_g0 and X_co20 in two consecutive iterations is less than the tolerance Tol, or until the maximum number of iterations max_count is reached.

IC_Finder(composition::Silicic, M_h2o::Float64, M_co2::Float64, M_tot::Float64, P::Float64, T::Float64, V::Float64, rho_m::Float64, param_IC::ParamICFinder{Float64})::NamedTuple{(:eps_g0, :X_co20, :mco2_diss, :phase),NTuple{4,Float64}}

An iterative function that uses thermodynamic modeling to determine the gas phase composition of a silicic magma chamber.

Arguments

• composition: Silicic()
• M_h2o: total mass of water in magma (kg)
• M_co2: total mass of CO2 in magma (kg)
• M_tot: total mass of magma (kg)
• P: Pressure (Pa)
• T: Temperature (K)
• V: chamber volume (m³)
• rho_m: density of the melt (kg/m³)
• param_IC: An instance of the ParamICFinder composite type containing the parameters for IC_Finder

Returns

A NamedTuple with the following fields:

• eps_g0: Initial gas fraction of the magma chamber
• X_co20: Initial mole fraction of CO2 in the gas phase
• mco2_diss: Total mass of CO2 dissolved in the melt (kg)
• phase: Integer representing the state of the magma

Details

The function uses a thermodynamic model to find the initial gas fraction (eps_g0) and initial mole fraction of CO2 in the gas phase (X_co20) for a silicic magma chamber. The model assumes that the chamber contains two phases: a liquid melt and a gas phase. The model calculates the saturation state of the magma with respect to CO2, and determines whether the gas phase contains only CO2 or a mixture of CO2 and H2O. If the chamber is in two-phase state, the model returns eps_g0 = 0.0, X_co20 = 0.0, mco2_diss = Total mass of CO2 dissolved in the melt, and phase = 2. If the chamber is in gas or liquid state, the function iteratively solves for eps_g0 and X_co20 using the solve_X_co2 and get_eps_g helper functions until the relative difference between eps_g0 and X_co20 in two consecutive iterations is less than the tolerance Tol, or until the maximum number of iterations max_count is reached.

This is for a11 and a12 a11: rho, drhodP, V, dVdP a12: rho, drhodT, V, dVdT

affect!(int, idx, sw::SW{Int8}, param::Param{Float64}, param_saved_var::ParamSaved{Float64}, param_IC_Finder::ParamICFinder{Float64})

Re-initialize the condition when the event happens. This function modifies the current state of the integrator (int.u) when a particular event occurs during the simulation. The function adjusts various parameters based on the current state of the integrator and the custom parameters that were passed in.

Arguments

• int: The current state of the integrator. It's format is from the DifferentialEquations.jl package
• idx: The index of the event that caused the function to be called.
• sw: A custom parameter used to control simulation behavior.
• param: A custom parameter containing physical constants and other model parameters.
• param_saved_var: A custom parameter used to store values from the previous time step.
• param_IC_Finder: A custom parameter used to control the behavior of the IC_Finder function.

The arguments int and idx are from the DifferentialEquations.jl package. These argument formats are specific to the DifferentialEquations.jl package.

boundary_conditions_new(P::Float64, T::Float64, V::Float64, rho_m::Float64, rho_x::Float64, c::Float64, sw::SW{Int8}, T_in::Float64, M_h2o::Float64, M_co2::Float64, total_Mass::Float64, param::Param{Float64}, param_saved_var::ParamSaved{Float64})

Arguments

• P: Pressure (Pa)
• T: Temperature (K)
• V: chamber volume (m³)
• rho_m: density of melt
• rho_x: density of crystal of magma
• c: heat of magma
• sw: eruption/coolingmodule/viscousrelaxation control
• T_in: Temperature
• M_h2o: total mess of H2O in the magma
• M_co2: total mess of CO2 in the magma
• total_Mass: total mess of magma chamber

build_co2(Pw::T, Pc::T, Temp::T, dPwdP::T, dPcdP::T, dPwdXco2::T, dPcdXco2::T)::NamedTuple{(:C_co2, :dC_co2dT, :dC_co2dP, :dC_co2dXco2),NTuple{4,T}} where {T<:Float64}

• This function is used within the exsolve function.

Returns

A NamedTuple with the following fields:

• C_co2
• dC_co2dT
• dC_co2dP
• dC_co2dXco2

build_mdot_in(fluxing::Bool, rho_m0::Float64, log_vfr::Float64, P_0::Float64, T_in::Float64)::Float64

• This function is used within the chamber function.

Returns

• mdot_in: mass inflow rate

build_meq(composition::Mafic, P::T, Temp::T, X_co2::T)::NamedTuple{(:meq, :dmeqdT, :dmeqdP, :dmeqdXco2),NTuple{4,T}} where {T<:Float64}

• This function is used within the exsolve function.

Returns

A NamedTuple with the following fields:

• meq
• dmeqdT
• dmeqdP
• dmeqdXco2

build_meq(composition::Silicic, Pw::T, Pc::T, Temp::T, dPwdP::T, dPcdP::T, dPwdXco2::T, dPcdXco2::T)::NamedTuple{(:meq, :dmeqdT, :dmeqdP, :dmeqdXco2),NTuple{4,T}} where {T<:Float64}

• This function is used within the exsolve function.

Returns

A NamedTuple with the following fields:

• meq
• dmeqdT
• dmeqdP
• dmeqdXco2

build_rho_rc(eps_m::T, eps_g::T, eps_x::T, rho_m::T, rho_g::T, rho_x::T, drho_m_dP::T, drho_g_dP::T, drho_x_dP::T, drho_m_dT::T, drho_g_dT::T, drho_x_dT::T, c_x::T, c_m::T, c_g::T, deps_x_dP::T, deps_x_dT::T)::Vector{T} where {T<:Float64}

• This function is used within the odeChamber function.

Returns

[rho, drho_dP, drho_dT, drho_deps_g, rc, drc_dP, drc_dT]

chamber(composition::Union{Silicic,Mafic}, end_time::Float64, log_volume_km3_vector::Union{Float64,Vector{Float64}}, InitialConc_H2O_vector::Union{Float64,Vector{Float64}}, InitialConc_CO2_vector::Union{Float64,Vector{Float64}}, log_vfr_vector::Union{Float64,Vector{Float64}}, depth_vector::Union{Float64,Vector{Float64}}, output_dirname::String; kwargs...)

Simulate the eruption of a volcano using a model for the frequency of eruptions of upper crustal magma chambers based on Degruyter and Huber (2014).

Arguments

• composition: The magma composition. Use Silicic() for rhyolite composition (arc magma) or Mafic() for basalt composition (rift).
• end_time: Maximum magma chamber evolution duration in seconds.
• log_volume_km3: The initial volume of the chamber in logarithmic scale. The actual initial chamber volume is calculated as 10^(log_volume_km3) in km³.
• InitialConc_H2O: The initial weight fraction of water in the magma (exsolved + dissolved).
• InitialConc_CO2: The initial weight fraction of CO₂ in the magma (exsolved + dissolved).
• log_vfr: Magma recharge rate in km³/yr calculated as 10^(log_vfr).
• depth: Depth of the magma chamber in meters.
• output_dirname(optional): Name of the output directory. Defaults to current timestamp.

Keyword Arguments

• plotfig(optional): (default: true). Generate and plot figures for each result if true.

Returns

A DataFrame containing the solution with columns:

• time: Simulation timestamps in seconds.
• P+dP: Pressure in Pa.
• T: Temperature in K.
• eps_g: Gas volume fraction.
• V: Volume of the magma chamber in m³.
• rho_m: Density of the melt in kg/m³.
• rho_x: Density of magma crystal in kg/m³.
• X_CO2: Mole fraction of CO2 in the gas.
• total_mass: Total mass of magma chamber in kg.
• total_mass_H2O: Total mass of water in the magma in kg.
• total_mass_CO2: Total mass of CO₂ in the magma in kg.
• eps_x: Crystal volume fraction.

Outputs

A directory named after output_dirname or the default value, containing the following files:

• out.csv: a CSV file containing the solution columns listed above.
• eruptions.csv, A CSV file containing the datas of eruptions with the following columns: timeoferuption (sec), durationoferuption (sec), masserupted (kg) and volumeerupted (km³).
• Figures for P+dP(t), T(t), epsg(t), V(t), XCO2(t), total_mass(t).

References

• W. Degruyter and C. Huber (2014). A model for eruption frequency of upper crustal silicic magma chambers. Earth Planet. Sci. Lett. https://doi.org/10.1016/j.epsl.2014.06.047

Examples

# Run a simulation with silicic magma chamber
julia> composition = Silicic()
julia> end_time = 3e9
julia> log_volume_km3 = 0.2
julia> InitialConc_H2O = [0.04, 0.045]
julia> InitialConc_CO2 = [0.001, 0.0012]
julia> log_vfr = -3.3
julia> depth = 8e3

julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth)

# Run a simulation with mafic magma chamber, with custom directory name "MyDirname"
julia> composition = Mafic()
julia> end_time = 3e9
julia> log_volume_km3 = 0.2
julia> InitialConc_H2O = [0.01, 0.012]
julia> InitialConc_CO2 = 0.001
julia> log_vfr = -3.3
julia> depth = [7e3, 8e3]
julia> output_dirname = "MyDirname"

julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth, output_dirname)

chamber(composition::Union{Silicic,Mafic}, end_time::Float64, log_volume_km3::Float64, InitialConc_H2O::Float64, InitialConc_CO2::Float64, log_vfr::Float64, depth::Float64, output_dirname::String; kwargs...)

Simulate the eruption of a volcano using a model for the frequency of eruptions of upper crustal magma chambers based on Degruyter and Huber (2014).

Arguments

• composition: The magma composition. Use Silicic() for rhyolite composition (arc magma) or Mafic() for basalt composition (rift).
• end_time: Maximum magma chamber evolution duration in seconds.
• log_volume_km3: The initial volume of the chamber in logarithmic scale. The actual initial chamber volume is calculated as 10^(log_volume_km3) in km³.
• InitialConc_H2O: The initial weight fraction of water in the magma (exsolved + dissolved).
• InitialConc_CO2: The initial weight fraction of CO₂ in the magma (exsolved + dissolved).
• log_vfr: Magma recharge rate in km³/yr calculated as 10^(log_vfr).
• depth: Depth of the magma chamber in meters.
• output_dirname(optional): Name of the output directory. Defaults to current timestamp.

Keyword Arguments

• plotfig(optional): (default: true). Generate and plot figures for each result if true.

Returns

A DataFrame containing the solution with columns:

• time: Simulation timestamps in seconds.
• P+dP: Pressure in Pa.
• T: Temperature in K.
• eps_g: Gas volume fraction.
• V: Volume of the magma chamber in m³.
• rho_m: Density of the melt in kg/m³.
• rho_x: Density of magma crystal in kg/m³.
• X_CO2: Mole fraction of CO2 in the gas.
• total_mass: Total mass of magma chamber in kg.
• total_mass_H2O: Total mass of water in the magma in kg.
• total_mass_CO2: Total mass of CO₂ in the magma in kg.
• eps_x: Crystal volume fraction.

Outputs

A directory named after output_dirname or the default value, containing the following files:

• out.csv: a CSV file containing the solution columns listed above.
• eruptions.csv, A CSV file containing the datas of eruptions with the following columns: timeoferuption (sec), durationoferuption (sec), masserupted (kg) and volumeerupted (km³).
• Figures for P+dP(t), T(t), epsg(t), V(t), XCO2(t), total_mass(t).

References

• W. Degruyter and C. Huber (2014). A model for eruption frequency of upper crustal silicic magma chambers. Earth Planet. Sci. Lett. https://doi.org/10.1016/j.epsl.2014.06.047

Examples

# Run a simulation with silicic magma chamber
julia> composition = Silicic()
julia> end_time = 3e9
julia> log_volume_km3 = 0.2
julia> InitialConc_H2O = 0.04
julia> InitialConc_CO2 = 0.001
julia> log_vfr = -3.3
julia> depth = 8e3

julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth)

# Run a simulation with mafic magma chamber, with custom directory name "MyDirname"
julia> composition = Mafic()
julia> end_time = 3e9
julia> log_volume_km3 = 0.2
julia> InitialConc_H2O = 0.01
julia> InitialConc_CO2 = 0.001
julia> log_vfr = -3.3
julia> depth = 8e3
julia> output_dirname = "MyDirname"

julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth, output_dirname)

check_for_duplicates(log_volume_km3_vector::Union{Float64,Vector{Float64}}, InitialConc_H2O_vector::Union{Float64,Vector{Float64}}, InitialConc_CO2_vector::Union{Float64,Vector{Float64}}, log_vfr_vector::Union{Float64,Vector{Float64}}, depth_vector::Union{Float64,Vector{Float64}})::Nothing

This function checks if any of the input arrays contain duplicate elements. If duplicates are found, it raises an error indicating which arguments have duplicates.

• This function is used within the chamber function

compute_dXdP_dXdT(u::Float64, param::Param, var::String)::Tuple{Float64,Float64,Float64}

• This function is used within the odeChamber function.

crystal_fraction(composition::Mafic, T::Float64, P::Float64, mH2O::Float64, mCO2::Float64)::NamedTuple{(:eps_x, :deps_x_dP, :deps_x_dT, :deps_x_deps_g, :deps_x_dmco2_t, :deps_x_dmh2o_t), NTuple{6, Float64}}

Calculate crystal volume fraction and its derivatives in magma at given temperature, pressure, and water and CO2 contents.

Arguments

• composition: Mafic()
• T: Temperature (K)
• P: Pressure (Pa)
• mH2O: Weight fration of the H2O in magma.
• mCO2: Weight fration of the CO2 in magma.

Returns

A NamedTuple with the following fields:

• eps_x: Crystal volume fraction.
• deps_x_dP: Derivative of eps_x with respect to pressure.
• deps_x_dT: Derivative of eps_x with respect to temperature.
• deps_x_deps_g: Derivative of eps_x with respect to gas volume fraction.
• deps_x_dmco2_t: Derivative of eps_x with respect to CO2 content.
• deps_x_dmh2o_t: Derivative of eps_x with respect to H2O content.

crystal_fraction(composition::Silicic, T::Float64, P::Float64, mH2O::Float64, mCO2::Float64)::NamedTuple{(:eps_x, :deps_x_dP, :deps_x_dT, :deps_x_deps_g, :deps_x_dmco2_t, :deps_x_dmh2o_t), NTuple{6, Float64}}

Calculate crystal volume fraction and its derivatives in magma at given temperature, pressure, and water and CO2 contents.

Arguments

• composition: Silicic()
• T: Temperature (K)
• P: Pressure (Pa)
• mH2O: Weight fration of the H2O in magma.
• mCO2: Weight fration of the CO2 in magma.

Returns

A NamedTuple with the following fields:

• eps_x: Crystal volume fraction.
• deps_x_dP: Derivative of eps_x with respect to pressure.
• deps_x_dT: Derivative of eps_x with respect to temperature.
• deps_x_deps_g: Derivative of eps_x with respect to gas volume fraction.
• deps_x_dmco2_t: Derivative of eps_x with respect to CO2 content.
• deps_x_dmh2o_t: Derivative of eps_x with respect to H2O content.

crystal_fraction_eps_x(composition::Mafic, T::Float64, P::Float64, mH2O::Float64, mCO2::Float64)::Float64

• Calculates crystal volume fraction at given temperature, pressure, and water and CO2 contents.
• Spetialized version of crystal_fractionthat computeseps_x` only.

Returns

• eps_x: Crystal volume fraction.

crystal_fraction_eps_x(composition::Silicic, T::Float64, P::Float64, mH2O::Float64, mCO2::Float64)::Float64

• Calculates crystal volume fraction at given temperature, pressure, and water and CO2 contents.
• Spetialized version of crystal_fraction that computes eps_x only.

Returns

• eps_x: Crystal volume fraction.

dM_X_t_dt_f(rho, V, Mdot_v_in, Mdot_v_out, m_h2o, Mdot_in, Mdot_out)

For dM_h2o_t_dt & dM_co2_t_dt

• MdotXin: Mdot_v_in or Mdot_c_in
• MdotXout: Mdot_v_out or Mdot_c_out
• mX: `mh2oorm_co2`

dX_dxdydz(composition::Union{Silicic,Mafic}, s::String, x::T, y::T, z::T)::NamedTuple{(:X, :dXdx, :dXdy, :dXdz),NTuple{4,T}} where {T<:Float64}

• This function is used within the parameters_melting_curve function.
• Calculate value from different parameter sets (# of parameter sets: Silicic: 3, Mafic: 2)

Arguments

• composition: Silicic() or Mafic()
• s: String that represents the parameter set to use. For Silicic, s can be "a", "b", or "c". For Mafic, s can be "a" or "b".
• x: Water (H2O)
• y: Gas (CO2)
• z: Pressure (P)

Returns

A NamedTuple with the following fields:

• X: the calculated value for the given parameter set
• dXdx: the derivative of X with respect to x
• dXdy: the derivative of X with respect to y
• dXdz: the derivative of X with respect to z

drc_dX_f(;eps_x::T, c_x::T, drho_x_dX::T, eps_g::T, c_g::T, drho_g_dX::T, eps_m::T, c_m::T, drho_m_dX::T, rho_x::T, rho_m::T, deps_x_dX::T)::T where {T<:Float64}

Computing the derivatives of the product of density and specific heat for the mixture

• X: represent P or T
• This function is used within the build_rho_rc function.

Returns

• drc_dX: represent drc_dP or drc_dT

drho_dX_f(;eps_m::T, eps_g::T, eps_x::T, drho_m_dX::T, drho_g_dX::T, drho_x_dX::T, rho_x::T, rho_m::T, deps_x_dX::T)::T where {T<:Float64}

• X: represent P or T
• This function is used within the build_rho_rc function.

Returns

• drho_dX: represent drho_dP or drho_dT

dwater_dx(composition::Mafic, p::T, t::T, x::T)::T where {T<:Float64}

• Derivative of Water Partitioning Function with respect to X_CO2
• This function is used within the exsolve3 function.

Arguments

• p: pressure (Pa)
• t: temperature (K)
• x: The previous mole fraction of CO2 (X_CO2)

dwater_dx(composition::Silicic, p::T, t::T, x::T)::T where {T<:Float64}

• Function for the derivative of water with reference to X_CO2
• This function is used within the exsolve3 function.

Arguments

• p: pressure (Pa)
• t: temperature (K)
• x: The previous mole fraction of CO2 (X_CO2)

eos_g(P::Float64, T::Float64)::EosG{Float64}

parametrization of redlich kwong taken from Huber et al. 2010

Arguments

• P: Pressure (Pa)
• T: Temperature (K)

Returns

An instance of the EosG struct, containing the following fields:

• rho_g: density of the gas (kg/m³).
• drho_g_dP: partial derivative of the density with respect to pressure (kg/m³/Pa).
• drho_g_dT: partial derivative of the density with respect to temperature (kg/m³/K).

eos_g_rho_g(P::Float64, T::Float64)::Float64

Spetialized version of eos_g that computes rho_g only.

Arguments

• P: Pressure (Pa)
• T: Temperature (K)

Returns

• rho_g: density of the gas in kg/m³.

exsolve(composition::Mafic, P::Float64, T::Float64, X_co2::Float64)::NamedTuple{(:meq, :dmeqdP, :dmeqdT, :dmeqdXco2, :C_co2, :dC_co2dP, :dC_co2dT, :dC_co2dXco2), NTuple{8, Float64}}

This script uses Liu et al. (2006) to calculate the solubility of water

Arguments

• composition: Mafic()
• P: Pressure (Pa)
• T: Temperature (K)
• X_co2: mole fraction of CO2 in gas.

exsolve(composition::Silicic, P::Float64, T::Float64, X_co2::Float64)::NamedTuple{(:meq, :dmeqdP, :dmeqdT, :dmeqdXco2, :C_co2, :dC_co2dP, :dC_co2dT, :dC_co2dXco2), NTuple{8, Float64}}

This script uses Liu et al. (2006) to calculate the solubility of water

Arguments

• composition: Silicic()
• P: Pressure (Pa)
• T: Temperature (K)
• X_co2: mole fraction of CO2 in gas.

exsolve3(composition::Mafic, P::Float64, T::Float64, m_eq::Float64)::NamedTuple{(:C_co2, :X_co2), NTuple{2, Float64}}

Takes pressure, temperature, and amount of water to solve for the concentration of CO2 and X_CO2 using a Newton Raphson scheme

Arguments

• composition: Mafic()
• P: pressure (Pa)
• T: temperature (K)
• m_eq: amount of water

Returns

A NamedTuple with the following fields:

• C_co2: The concentration of CO2 in the melt.
• X_co2: The mole fraction of CO2 in the gas.

exsolve3(composition::Silicic, P::Float64, T::Float64, m_eq::Float64)::NamedTuple{(:C_co2, :X_co2), NTuple{2, Float64}}

Takes pressure, temperature, and amount of water to solve for the concentration of CO2 and X_CO2 using a Newton Raphson scheme

Arguments

• composition: Silicic()
• P: pressure (Pa)
• T: temperature (K)
• m_eq: amount of water

Returns

A NamedTuple with the following fields:

• C_co2: The concentration of CO2 in the melt.
• X_co2: The mole fraction of CO2 in the gas.

exsolve_meq(composition::Mafic, P::Float64, T::Float64, X_co2::Float64)::Float64

Spetialized version of exsolve that computes meq only.

Arguments

• composition: Mafic()
• P: Pressure (Pa)
• T: Temperature (K)
• X_co2: mole fraction of CO2 in gas.

Returns

• meq: amount of water

exsolve_meq(composition::Silicic, P::Float64, T::Float64, X_co2::Float64)::Float64

Spetialized version of exsolve that computes meq only.

Arguments

• composition: Silicic()
• P: Pressure (Pa)
• T: Temperature (K)
• X_co2: mole fraction of CO2 in gas.

Returns

• meq: amount of water

find_liq(composition::Mafic, water::Float64, co2::Float64, P::Float64, ini_eps_x::Float64)::Float64

Given the weight fractions of H2O and CO2 in magma, and the pressure, finds the temperature at which the magma will crystallize and reach the desired volume fraction of crystals. The function uses the parameters of the melting curve and an error function to estimate the liquidus temperature.

Arguments

• composition: Mafic()
• water: Weight fration of the H2O in magma.
• co2: Weight fration of the CO2 in magma.
• P: Pressure (Pa)
• ini_eps_x: The starting volumn fraction of crystal.

Returns

• Tl: the estimated liquidus temperature (K).

find_liq(composition::Silicic, water::Float64, co2::Float64, P::Float64, ini_eps_x::Float64)::Float64

Given the weight fractions of H2O and CO2 in magma, and the pressure, finds the temperature at which the magma will crystallize and reach the desired volume fraction of crystals. The function uses the parameters of the melting curve and an error function to estimate the liquidus temperature.

Arguments

• composition: Silicic()
• water: Weight fration of the H2O in magma.
• co2: Weight fration of the CO2 in magma.
• P: Pressure (Pa)
• ini_eps_x: The starting volumn fraction of crystal.

Returns

• Tl: the estimated liquidus temperature (K).

fun(x::Float64, P::Float64, T::Float64, eps_g0::Float64, eps_x0::Float64, V::Float64, rho_m::Float64, rho_g::Float64, mm_co2::Float64, mm_h2o::Float64, M_co2::Float64)::Float64

• This function is used within the solve_X_co2 function to solve X_co20.

gas_heat_capacity(X_co2::Float64)

Arguments

• X_co2: mole fraction of CO2 in the gas.

get_eps_g(composition::Union{Silicic,Mafic}, eps_g_prev::Float64, X_co20::Float64, P::Float64, T::Float64, eps_x0::Float64, V::Float64, rho_m::Float64, rho_g::Float64, M_h2o::Float64, M_co2::Float64)::NamedTuple{(:eps_g0, :mco2_diss),NTuple{2,Float64}}

• This function is used within the IC_Finder function to get epsg0 and mco2diss.

Returns

A NamedTuple with the following fields:

• eps_g0
• mco2_diss

get_phase(composition::Union{Silicic,Mafic}, P::Float64, T::Float64, V::Float64, rho_m::Float64, M_h2o::Float64, M_co2::Float64, eps_x0::Float64)::NamedTuple{(:phase, :m_co2_melt),NTuple{2,Float64}}

• This function is used within the IC_Finder function to determine the initial phase.

Returns

A NamedTuple with the following fields:

• phase: 2 or 3
• m_co2_melt

heat_conduction_chamberCH(maxn::Int64, a::Float64, c::Float64, dr::Float64, kappa::Float64, rho::Float64, cp::Float64, Tb::Float64, param_sv::ParamSaved{Float64})::Float64

The function calculates the heat conduction within a magma chamber by solving the heat equation in cylindrical coordinates.

Arguments

• maxn: number of terms
• a: radius of magma chamber (m)
• c: radius of outer shell (m)
• dr: grid spacing for calculating the heat transfer
• kappa: thermal diffusivity of host rocks (m^2/s)
• rho: density of the magma (kg/m³)
• cp: specific heat of magma (J/(kg K))
• Tb: boundary temperature of the outer shell (K)
• param_sv: a struct ParamSaved stores the previous results of the function

Returns

• Q: heat flow rate through the magma chamber (J/s)

heat_conduction_chamber_profileCH(maxn::Int64, a::Float64, c::Float64, r::Float64, kappa::Float64, Tb::Float64, param_sv::ParamSaved{Float64})::Float64

Calculates the heat transfer in a magma chamber based on the thermal diffusivity of the host rocks and the temperature profile of the outer shell.

Arguments

• maxn: number of terms
• a: radius of magma chamber (m)
• c: radius of outer shell (m)
• r: grid spacing of calculate the heat transfer
• kappa: thermal diffusivity of host rocks
• Tb: boundary temperature of the outer shell (K)
• param_sv: a struct ParamSaved stores the previous results of the function

Returns

• Trt

ic_phase_conversion(phase_here::T, composition::Union{Silicic,Mafic}, M_h2o::T, M_co2::T, M_tot::T, P::T, Temp::T, V::T, rho_m::T, param_IC::ParamICFinder{T})::NamedTuple{(:eps_g_temp, :X_co2_temp, :C_co2_temp, :phase),NTuple{4,T}} where {T<:Float64}

This function is used to handle phase conversion by iteratively modifying the max_count or Tol parameters of function IC_Finder until the correct phase is obtained.

• This function is used within the affect! function.

mco2_dissolved_sat(X::Float64, P::Float64, T::Float64)::Float64

• This function is used within the solve_X_co2 function

Arguments

• X: mole fraction of CO2 in gas
• P: Pressure (Pa)
• T: Temperature (K)

meq_water(composition::Mafic, X::Float64, P::Float64, T::Float64)::Float64

• This function is used within the get_eps_g function

Arguments

• X: mole fration of H2O in gas
• P: Pressure (Pa)
• T: Temperature (K)

meq_water(composition::Silicic, X::Float64, P::Float64, T::Float64)::Float64

• This function is used within the get_eps_g function

Arguments

• X: mole fration of H2O in gas
• P: Pressure (Pa)
• T: Temperature (K)

odeChamber(du::Vector{Float64}, u::Vector{Float64}, params::Tuple{Param{Float64}, ParamSaved{Float64}, SW{Int8}}, t::Float64)

• Define the ODE equation.
• Solve the model for eruption frequency of upper crustal magma chambers using an ODE solver.

Arguments

• du: An array that stores the output of the ODE solver, i.e., the values of the derivatives of the solution u with respect to time t.
• u: An array that stores the values of the solution at each time step t.
• p: A tuple that stores the model parameters and some saved variables, which are described in more detail below.
• t: The time points corresponding to the saved values of the ODE solution.

The arguments du, u, p, and t are from the DifferentialEquations.jl package. These argument formats are specific to the DifferentialEquations.jl package.

Parameters

• param: A custom parameter containing physical constants and other model parameters.
• param_saved_var: A custom parameter used to store values from the previous time step.
• sw: A custom parameter used to control simulation behavior.

Returns

The function modifies du in place to store the values of the derivatives of the solution u with respect to time t.

parameters_melting_curve(composition::Mafic, mH2O::Float64, mCO2::Float64, P::Float64)::NamedTuple{(:a, :dadx, :dady, :dadz, :b, :dbdx, :dbdy, :dbdz), NTuple{8, Float64}}

• This function is used within the find_liq, crystal_fraction, crystal_fraction_eps_x function.

Arguments

• composition: Mafic()
• mH2O: Weight fration of the H2O in magma.
• mCO2: Weight fration of the CO2 in magma.
• P: Pressure (Pa)

parameters_melting_curve(composition::Silicic, mH2O::Float64, mCO2::Float64, P::Float64)::NamedTuple{(:a, :dadx, :dady, :dadz, :b, :dbdx, :dbdy, :dbdz, :c, :dcdx, :dcdy, :dcdz), NTuple{12, Float64}}

• This function is used within the find_liq, crystal_fraction, crystal_fraction_eps_x function.

Arguments

• composition: Silicic()
• mH2O: Weight fration of the H2O in magma.
• mCO2: Weight fration of the CO2 in magma.
• P: Pressure (Pa)

rc_f(;rho_x::T, eps_x::T, c_x::T, rho_m::T, eps_m::T, c_m::T, rho_g::T, eps_g::T, c_g::T)::T where {T<:Float64}

Computing the product of density and specific heat

• This function is used within the build_rho_rc function.

Returns

• rc: the product of density and specific heat

record_erupt_end(time::T, erupt_saved::EruptSaved{T}, param::Param{T}) where {T<:Float64}

Record duration, mass and volume of eruptions to EruptSaved

• This function is used within the affect! function.

record_erupt_start(time::T, eps_g::T, eps_x::T, rho_m::T, rho_x::T, rho_g::T, erupt_saved::EruptSaved{T}) where {T<:Float64}

Record time, density of eruptions to EruptSaved

• This function is used within the affect! function.

rho_0_f(eps_g0::Float64, eps_x0::Float64, rho_g0::Float64, rho_m0::Float64, rho_x0::Float64)::Float64

• This function is used within the chamber function.

Returns

• rho_0: initial bulk density (kg/m³)

rho_f(;eps_m::T, eps_g::T, eps_x::T, rho_m::T, rho_g::T, rho_x::T)::T where {T<:Float64}

• This function is used within the build_rho_rc function.

Returns

• rho: bulk density (kg/m³)

solve_NR(f, f_prime, errorTol::T, count_max::T, Xc_initial::T)::T where {T<:Float64}

• Finding mole fraction of CO2 in gas (X_CO2).
• This function is used within the exsolve3 function.

Arguments

• f: Water Paritioning Function
• f_prime: Function for the derivative of water with reference to X_CO2
• errorTol: error tolerance, the default value is 1e-10
• count_max: Maximum loop count, the default value is 1e2
• Xc_initial: initial guess of X_CO2, the dedault value is 1e-2

Returns

• X_CO2: mole fraction of CO2 in gas

solve_X_co2(composition::Mafic, eps_g0::Float64, X_co2_prev::Float64, P::Float64, T::Float64, eps_x0::Float64, V::Float64, rho_m::Float64, rho_g::Float64, M_co2::Float64, Tol::Float64)::Float64

• This function is used within the IC_Finder function to solve X_co20.

Returns

• X_co20

solve_X_co2(composition::Silicic, eps_g0::Float64, X_co2_prev::Float64, P::Float64, T::Float64, eps_x0::Float64, V::Float64, rho_m::Float64, rho_g::Float64, M_co2::Float64, Tol::Float64)::Float64

• This function is used within the IC_Finder function to solve X_co20.

Returns

• X_co20

stopChamber_MT(out, u::Vector{Float64}, t::Float64, int, sw::SW{Int8}, param::Param{Float64})

Define the stopping criteria for an ODE solver that simulates a magma chamber.

Arguments

• out: An array where the function should save the condition value at the right index. The maximum index of out should be specified in the len property of callback, which allows for a chain of len events, triggering the ith event when out[i] = 0. The function returns the value of out[8] as the last condition. Checking Event Handling and Callback Functions page of DifferentialEquations.jl for more details.
• u: A vector containing the state of the system at time t.
• t: The current time of the ODE solver.
• int: The current state of the integrator. It's format is from the DifferentialEquations.jl package
• sw: A custom parameter used to control simulation behavior.
• param: A custom parameter containing physical constants and other model parameters.

The arguments out, u, t, and int are from the DifferentialEquations.jl package. These argument formats are specific to the DifferentialEquations.jl package.

Returns

The out array is modified in-place to contain the condition values at the current state of the system. The function computes several quantities based on the current state of the system and the model parameters, such as the crystal fraction, the maximum amount of exsolved fluid, and the amount of CO2 in the melt. It then uses these values to calculate the condition values to be stored in the out array, as follows:

• out[1]: the crystal fraction.
• out[2]: the ratio of crystal fraction to liquid fraction, minus 0.8.
• out[3]: if an eruption is not occurring, the pressure difference between the lithostatic pressure and the current pressure, minus the critical pressure difference. Otherwise, negative the critical pressure difference.
• out[4]: if an eruption is occurring, the lithostatic pressure minus the current pressure. Otherwise, negative the critical pressure difference.
• out[5]: the crystal fraction minus 0.5.
• out[6]: the difference between the amount of water in the melt and the maximum amount of exsolved fluid.
• out[7]: negative the difference between the initial lithostatic pressure and the current pressure plus the critical pressure difference.
• out[8]: the difference between the amount of CO2 in the melt and the saturation concentration.

Note that the out and u arguments are in the format expected by the DifferentialEquations.jl package, and the function is intended to be used as a condition for a callback function.

water(composition::Silicic, p::T, t::T, x::T, c::T)::T where {T<:Float64}

• Water Paritioning Function
• This function is used within the exsolve3 function.

Arguments

• p: pressure (Pa)
• t: temperature (K)
• x: The previous mole fraction of CO2 (X_CO2)
• c: amount of water

water(composition::Silicic, p::T, t::T, x::T, c::T)::T where {T<:Float64}

• Water Paritioning Function
• This function is used within the exsolve3 function.

Arguments

• p: pressure (Pa)
• t: temperature (K)
• x: The previous mole fraction of CO2 (X_CO2)
• c: amount of water

Settings for ODE Solver Calculation

Parameters

Settings for IC_Finder

Rheology of the crust

eruption/coolingmodule/viscousrelaxation control