earthsciencekit.volcanology

Domain hub · History · Examples · API reference

Inspectable volcanology and geochemistry models; SI units unless stated otherwise.

earthsciencekit.volcanology.arrhenius_viscosity(temperature, reference_viscosity, activation_energy, reference_temperature)

Arrhenius temperature dependence of melt viscosity.

Silicate melts stiffen by orders of magnitude as they cool, which controls how lava flows and how magma erupts.

Parameters:

• temperature (float or array_like) – K; constraint: positive.

• reference_viscosity (float or array_like) – Pa s, viscosity at the reference temperature; constraint: positive.

• activation_energy (float or array_like) – J/mol; constraint: nonnegative.

• reference_temperature (float or array_like) – K; constraint: positive.

Returns:

Model value [Pa s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Constant activation energy over the range of interest; crystal- and bubble-free Newtonian melt. Equation: reference_viscosity * np.exp(activation_energy / R_GAS * (1 / temperature - 1 / reference_temperature)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.column_height(heat_flux)

Height of a Plinian eruption column from its heat flux.

The thermal energy released at the vent, not the eruption velocity, controls how high a volcanic column rises.

Parameters:

heat_flux (float or array_like) – W, rate of heat release mass_eruption_rate * specific_heat * temperature_excess; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Wilson, Sparks, Huang and Watkins (1978), JGR 83: 1829-1836, for a standard atmosphere and a strong steady column; wind bends weak plumes lower. Equation: 8.2 * heat_flux**0.25. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.concordia_ratios(time)

Point on the Wetherill U-Pb concordia.

Both uranium-lead clocks agree only on this curve, so a zircon off it has lost lead and its chord dates the loss.

Parameters:

time (float or array_like) – s; constraint: nonnegative.

Returns:

Model value [1]: [207Pb*/235U, 206Pb*/238U], the horizontal and vertical concordia coordinates.

Return type:

ndarray, shape (…, 2)

Notes

Wetherill (1956): closed system with the Jaffey et al. (1971) decay constants. Episodic lead loss puts samples on a straight chord (discordia) between the concordia points of formation and loss. Equation: [np.expm1(LAMBDA_U235 * time), np.expm1(LAMBDA_U238 * time)]. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.decay_constant(half_life)

Decay constant from a half-life.

Every radioactive nuclide has a fixed probability per unit time of decaying, set by its half-life.

Parameters:

half_life (float or array_like) – s; constraint: positive.

Returns:

Model value [1/s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Rutherford and Soddy (1902): first-order decay, independent of temperature, pressure and chemistry. Equation: np.log(2) / half_life. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.delta_value(ratio, standard_ratio)

Isotope delta value relative to a standard.

Isotope ratios vary in the fourth decimal place, so geochemists report deviations from a standard in per mil.

Parameters:

• ratio (float or array_like) – 1, heavy-to-light isotope ratio of the sample; constraint: positive.

• standard_ratio (float or array_like) – 1, the same ratio in the standard (e.g. VSMOW, VPDB); constraint: positive.

Returns:

Model value [per mil]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Linear in the ratio; delta values of different standards must be converted before they are compared. Equation: (ratio / standard_ratio - 1) * 1000. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.dike_ascent_velocity(width, density_contrast, viscosity, gravity)

Mean laminar velocity of buoyant magma in a dike.

Buoyancy drives magma up a planar crack, and the viscous drag on its walls sets how fast it rises.

Parameters:

• width (float or array_like) – m, full dike width; constraint: positive.

• density_contrast (float or array_like) – kg/m^3, host rock minus magma; constraint: real.

• viscosity (float or array_like) – Pa s; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m/s], positive upward. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Wilson and Head (1981): plane Poiseuille flow between parallel walls driven by buoyancy alone; check that the Reynolds number magma_density * velocity * width / viscosity stays below about 2000. Equation: width**2 * density_contrast * gravity / (12 * viscosity). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.equilibrium_crystallization_melt(initial_concentration, melt_fraction, partition_coefficient)

Trace-element concentration in a melt during equilibrium crystallization.

Crystals that keep reacting with the melt, as in Bowen’s continuous series, limit how far the melt can evolve.

Parameters:

• initial_concentration (float or array_like) – Concentration in the parent melt (any unit, e.g. ppm); constraint: nonnegative.

• melt_fraction (float or array_like) – 1, fraction of melt remaining; constraint: fraction.

• partition_coefficient (float or array_like) – 1, bulk solid/melt concentration ratio; constraint: positive.

Returns:

Model value [unit of initial_concentration]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Batch mass balance (Shaw 1970): solid and melt stay in equilibrium, so an incompatible element is enriched at most by 1 / partition_coefficient. Equation: initial_concentration / (partition_coefficient + melt_fraction * (1 - partition_coefficient)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.fractional_crystallization_melt(initial_concentration, melt_fraction, partition_coefficient)

Trace-element concentration in a melt during fractional crystallization.

Crystals removed as soon as they form leave the residual melt enriched in the elements they reject.

Parameters:

• initial_concentration (float or array_like) – Concentration in the parent melt (any unit, e.g. ppm); constraint: nonnegative.

• melt_fraction (float or array_like) – 1, fraction of melt remaining; constraint: fraction above zero.

• partition_coefficient (float or array_like) – 1, bulk solid/melt concentration ratio; constraint: nonnegative.

Returns:

Model value [unit of initial_concentration]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Rayleigh fractional crystallization (Bowen 1928; Gast 1968): crystals do not re-equilibrate and the bulk partition coefficient is constant. Equation: initial_concentration * melt_fraction**(partition_coefficient - 1). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.isochron_fit(parent_ratio, daughter_ratio, decay_constant)

Age and initial ratio from a whole-rock or mineral isochron.

Cogenetic minerals with different parent contents line up on a straight line whose slope gives the age without knowing the initial daughter.

Parameters:

• parent_ratio (array_like, shape (n,), n >= 2) – 1, parent over a stable reference isotope (e.g. 87Rb/86Sr); constraint: nonnegative, not all equal.

• daughter_ratio (array_like, shape (n,)) – 1, daughter over the same reference isotope (e.g. 87Sr/86Sr); constraint: nonnegative.

• decay_constant (float) – 1/s; constraint: positive scalar.

Returns:

[age, initial_ratio]: the age [s] from the least-squares slope and the intercept [1].

Return type:

numpy.ndarray, shape (2,)

Notes

Nicolaysen (1961): every sample formed at the same time with the same initial ratio and stayed closed; ordinary least squares without measurement-error weighting (a York fit would weight them). The slope must be positive. Equation: daughter_ratio = initial_ratio + parent_ratio * np.expm1(decay_constant * age). Invalid, empty, nonreal or nonfinite inputs raise ValueError.

earthsciencekit.volcanology.isotope_temperature(delta_calcite, delta_water)

Calcite-water oxygen-isotope paleotemperature.

The oxygen-isotope fractionation between a shell and the water it grew in falls with temperature, so fossil shells record ancient ocean temperatures.

Parameters:

• delta_calcite (float or array_like) – per mil, delta18O of the carbonate (CO2 from the calcite, PDB scale); constraint: real.

• delta_water (float or array_like) – per mil, delta18O of the water on the same scale; constraint: real.

Returns:

Model value [K]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Urey (1947) predicted the effect; this is the empirical calibration of Epstein, Buchsbaum, Lowenstam and Urey (1953), valid for about 7-30 degC and equilibrium calcite. Equation: 273.15 + 16.5 - 4.3 * d + 0.14 * d**2 with d = delta_calcite - delta_water. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.lattice_strain_partition(ionic_radius, optimal_radius, peak_partition, youngs_modulus, temperature)

Mineral-melt partition coefficient from the strain of the crystal lattice.

Ions close in size to the site they replace enter a crystal easily, which is Goldschmidt’s first rule of substitution.

Parameters:

• ionic_radius (float or array_like) – m; constraint: positive.

• optimal_radius (float or array_like) – m, radius that fits the site without strain; constraint: positive.

• peak_partition (float or array_like) – 1, partition coefficient of an ion of the optimal radius; constraint: positive.

• youngs_modulus (float or array_like) – Pa, apparent Young’s modulus of the site; constraint: positive.

• temperature (float or array_like) – K; constraint: positive.

Returns:

Model value [1]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Blundy and Wood (1994) quantification of Goldschmidt (1937) for isovalent ions on one site; charge differences need a separate parabola. Equation: peak_partition * np.exp(-4 * np.pi * youngs_modulus * AVOGADRO * (optimal_radius / 2 * (ionic_radius - optimal_radius)**2 + (ionic_radius - optimal_radius)**3 / 3) / (R_GAS * temperature)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.lead_lead_age(lead_ratio)

Age from a radiogenic 207Pb/206Pb ratio.

Patterson dated the Earth from the slope of a meteorite lead isochron, which equals this ratio.

Parameters:

lead_ratio (float or array_like) – 1, radiogenic 207Pb*/206Pb* or the slope of a Pb-Pb isochron; constraint: between the ratios at one year (about 0.0461) and at 10 Ga.

Returns:

Model value [s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Inverts radiogenic_lead_ratio by Brent’s method between one year and 10 Ga; same closed-system assumptions. Equation: radiogenic_lead_ratio(age) = lead_ratio. Invalid, empty, nonreal or nonfinite inputs raise ValueError.

earthsciencekit.volcanology.mogi_displacement(radial_distance, depth, volume_change, poisson_ratio)

Surface displacement above an inflating point source.

Magma entering a chamber lifts the ground in a bell-shaped dome whose width reveals the chamber depth.

Parameters:

• radial_distance (float or array_like) – m, horizontal distance from the point above the source; constraint: nonnegative.

• depth (float or array_like) – m; constraint: positive.

• volume_change (float or array_like) – m^3, signed change of the source volume; constraint: real.

• poisson_ratio (float or array_like) – 1; constraint: between 0 and 0.5.

Returns:

Model value [m]: [radial, vertical] displacement, outward and upward positive.

Return type:

ndarray, shape (…, 2)

Notes

Mogi (1958) point source in an elastic half-space; the chamber radius must be small compared with its depth. The uplift halves at radial_distance = 0.766 * depth. Equation: (1 - poisson_ratio) * volume_change / np.pi * [radial_distance, depth] / (radial_distance**2 + depth**2)**1.5. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.parent_fraction(decay_constant, time)

Fraction of radioactive parent atoms remaining.

Exponential decay is the clock behind every radiometric age.

Parameters:

• decay_constant (float or array_like) – 1/s; constraint: positive.

• time (float or array_like) – s; constraint: nonnegative.

Returns:

Model value [1]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Rutherford and Soddy (1902): a closed system with no parent added after time zero. Equation: np.exp(-decay_constant * time). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.plume_rise_height(buoyancy_flux, buoyancy_frequency)

Maximum height of a buoyant plume in a stratified fluid.

Turbulent entrainment dilutes a plume until the stratification removes its buoyancy, which caps its height at a value set by only the source buoyancy flux and the stratification.

Parameters:

• buoyancy_flux (float or array_like) – m^4/s^3, source reduced_gravity * velocity * radius**2 (the factor pi omitted); constraint: positive.

• buoyancy_frequency (float or array_like) – 1/s; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Morton, Taylor and Turner (1956) for a steady point source of buoyancy in uniform stratification, with Briggs’s constant 5.0 (entrainment about 0.1); source momentum, wind and moisture are neglected. Equation: 5.0 * buoyancy_flux**0.25 * buoyancy_frequency**-0.75. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.radiocarbon_age(fraction_modern)

Conventional radiocarbon age.

Living matter exchanges carbon-14 with the atmosphere; after death the carbon-14 decays, dating organic remains.

Parameters:

fraction_modern (float or array_like) – 1, sample 14C/12C relative to the modern standard; constraint: positive.

Returns:

Model value [s]; negative for bomb carbon above the standard. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Libby et al. (1949) with the conventional 5568-year (Libby) half-life, mean life LIBBY_MEAN_LIFE; it is a radiocarbon age, not a calendar age, until calibrated. Equation: -LIBBY_MEAN_LIFE * np.log(fraction_modern). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.radiogenic_daughter(parent, decay_constant, time)

Radiogenic daughter accumulated beside the surviving parent.

Forward accumulation shows how the daughter-to-parent ratio of a rock records its age.

Parameters:

• parent (float or array_like) – Surviving parent amount today (atoms, moles or any consistent unit); constraint: nonnegative.

• decay_constant (float or array_like) – 1/s; constraint: positive.

• time (float or array_like) – s; constraint: nonnegative.

Returns:

Model value [unit of parent]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Closed system with no initial daughter; the inverse of radiometric_age. Holmes (1913) used it to build the first radiometric time scale. Equation: parent * np.expm1(decay_constant * time). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.radiogenic_lead_ratio(time)

Radiogenic 207Pb/206Pb ratio produced over a given age.

Two uranium isotopes decaying at different rates make the ratio of their lead daughters a clock that needs no uranium measurement.

Parameters:

time (float or array_like) – s; constraint: positive.

Returns:

Model value [1]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Patterson (1956): closed system with today’s 238U/235U = U238_U235 and the Jaffey et al. (1971) decay constants. The ratio tends to LAMBDA_U235 / (U238_U235 * LAMBDA_U238) at zero age. Equation: np.expm1(LAMBDA_U235 * time) / (U238_U235 * np.expm1(LAMBDA_U238 * time)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.radiometric_age(daughter_parent_ratio, decay_constant)

Age from a radiogenic daughter-to-parent ratio.

The ratio of accumulated daughter atoms to surviving parent atoms grows with time, so measuring it dates the rock.

Parameters:

• daughter_parent_ratio (float or array_like) – 1, radiogenic daughter atoms per parent atom; constraint: nonnegative.

• decay_constant (float or array_like) – 1/s; constraint: positive.

Returns:

Model value [s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Closed system with no initial daughter; Boltwood (1907) applied it to lead in uranium minerals. Equation: np.log1p(daughter_parent_ratio) / decay_constant. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.rayleigh_fractionation(initial_ratio, remaining_fraction, fractionation_factor)

Isotope ratio of a reservoir depleted by Rayleigh distillation.

Removing a product that is always slightly enriched in one isotope drives the remaining reservoir further and further from its starting composition.

Parameters:

• initial_ratio (float or array_like) – 1, heavy-to-light isotope ratio of the reservoir at the start; constraint: positive.

• remaining_fraction (float or array_like) – 1, fraction of the reservoir left; constraint: fraction above zero.

• fractionation_factor (float or array_like) – 1, product ratio over reservoir ratio at each instant; constraint: positive.

Returns:

Model value [1]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Rayleigh (1896): the product leaves as soon as it forms, the reservoir stays well mixed and the factor is constant. Equation: initial_ratio * remaining_fraction**(fractionation_factor - 1). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.shaw_viscosity(temperature, slope)

Shaw’s empirical viscosity of an anhydrous or hydrous silicate melt.

One composition-dependent slope predicts the viscosity of melts from basalt to rhyolite.

Parameters:

• temperature (float or array_like) – K; constraint: positive.

• slope (float or array_like) – 1, Shaw’s s from the molar composition (about 2-3 for basalt, 3.5-4 for rhyolite); constraint: positive.

Returns:

Model value [Pa s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Shaw (1972), Am. J. Sci. 272: 870-893; all lines pass near ln(eta [poise]) = -6.4 at 1e4/T = 1.5. Crystal- and bubble-free Newtonian melt; about 0.2 log units of accuracy. Equation: 0.1 * np.exp(slope * (1e4 / temperature) - 1.5 * slope - 6.4). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.simulate_eruption_column(vent_radius, vent_velocity, reduced_gravity, buoyancy_frequency, time, *, entrainment=0.12)

Follow a parcel up the axis of a steady buoyant eruption column.

Parameters:

• vent_radius (float) – Positive scalar source radius [m].

• vent_velocity (float) – Positive scalar source velocity [m/s].

• reduced_gravity (float) – Positive scalar source buoyancy g (rho_air - rho_plume) / rho_air [m/s^2].

• buoyancy_frequency (float) – Nonnegative scalar Brunt-Vaisala frequency of the atmosphere [1/s].

• time (array_like) – Nonnegative strictly increasing travel times since the parcel left the vent [s].

• entrainment (float) – Positive scalar top-hat entrainment coefficient [1]; 0.12 reproduces the constant 5.0 of plume_rise_height.

Returns:

Parcel height [m], vertical velocity [m/s] and reduced gravity [m/s^2].

Return type:

TimeSeries

Notes

The Morton, Taylor and Turner (1956) top-hat plume equations for the fluxes Q = b**2 w, M = b**2 w**2 and F = b**2 w g’ (pi omitted), carried by a parcel moving at w = M / Q: dz/dt = M / Q, dQ/dt = 2 alpha abs(M)**1.5 / Q, dM/dt = F and dF/dt = -N**2 M. In a uniform atmosphere the parcel stops near plume_rise_height. Past that height the steady theory ends; the parcel is continued as a fountain that falls back with the same entrainment law and oscillates about the level of neutral buoyancy, near three quarters of the maximum height, where the umbrella cloud spreads. This illustrative closure omits wind, the dense gas-thrust region, particle fallout and moisture.

earthsciencekit.volcanology.volcanic_explosivity_index(erupted_volume)

Volcanic Explosivity Index from the erupted tephra volume.

A single logarithmic number lets eruptions of very different style and size be compared and catalogued.

Parameters:

erupted_volume (float or array_like) – m^3, bulk tephra volume; constraint: nonnegative.

Returns:

Model value [1], an integer from 0 to 8. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Newhall and Self (1982), volume criterion only: VEI 0 below 1e4 m^3, 1 below 1e6 m^3, then one step per factor of ten up to VEI 8 at 1e12 m^3 and above. Column height and other criteria are not used. Equation: np.where(erupted_volume < 1e4, 0, np.where(erupted_volume < 1e6, 1, np.clip(np.floor(np.log10(erupted_volume)) - 4, 2, 8))). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.volcanology.water_solubility(pressure, solubility_constant)

Dissolved water in a melt saturated at a given pressure.

Water dissolves as hydroxyl, so its solubility grows as the square root of pressure and ascending magma exsolves bubbles.

Parameters:

• pressure (float or array_like) – Pa; constraint: nonnegative.

• solubility_constant (float or array_like) – 1/sqrt(Pa), mass fraction per root pascal (about 4.1e-6 for rhyolite); constraint: positive.

Returns:

Model value [1], mass fraction. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Square-root (Henry-type) law for pure water vapour, valid up to a few hundred MPa; ignores CO2 and temperature. Equation: solubility_constant * np.sqrt(pressure). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.