This notebook runs the StatGeochem.jl package, which implements an interface for interacting with Perple_X from the Julia programming language, including from Jupyter notebooks such as this.
StatGeochem.jl also includes some of the codes and utilities used in Keller & Schoene 2012, Keller et al. 2015 and Keller & Schoene 2018.
If running this notebook as an online Binder notebook and the webpage times out, click the badge at left to relaunch (refreshing will not work). Note that any changes will be lost!
Hint: shift
-enter
to run a single cell, or from the Cell
menu select Run All
to run the whole file. Any code from this notebook can be copied and pasted into the Julia REPL or a .jl
script.
## --- Load (and install if neccesary) the StatGeochem package which has the resampling functions we'll want
try
using StatGeochem
catch
using Pkg
Pkg.add("StatGeochem")
using StatGeochem
end
using Statistics, StatsBase, DelimitedFiles
using Plots; gr();
┌ Info: Recompiling stale cache file /Users/cbkeller/.julia/compiled/v1.1/StatGeochem/Ht7Cf.ji for StatGeochem [df4de05a-b714-11e8-3c2a-c30fb13e804c] └ @ Base loading.jl:1184
## --- Configure
# Absolute paths to perplex resources
perplexdir = joinpath(resourcepath,"perplex-stable")
scratchdir = "./scratch/" # Location of directory to store output files
# Attempt to install perplex, if not already extant
if !isfile(joinpath(perplexdir,"vertex"))
# Make sure resourcepath exists
run(`mkdir -p $resourcepath`)
# Try to compile PerpleX from source; if that fails, try to download linux binaries
try
# Check if there is a fortran compiler
run(`gfortran -v`)
# Download Perplex v6.8.7 -- known to work with interface used here
file = download("https://storage.googleapis.com/statgeochem/perplex-6.8.7-source.zip", joinpath(resourcepath,"perplex-stable.zip"))
# # For a more updated perplex version, you might also try
# file = download("https://petrol.natur.cuni.cz/~ondro/perplex-sources-stable.zip", joinpath(resourcepath,"perplex-stable.zip"))
run(`unzip -u $file -d $resourcepath`) # Extract
system("cd $perplexdir; make") # Compile
catch
@warn "Failed to compile from source, trying precompiled linux binaries instead"
run(`mkdir -p $perplexdir`)
file = download("https://petrol.natur.cuni.cz/~ondro/Perple_X_6.8.7_Linux_64_gfortran.tar.gz","perplex-6.8.7-linux.tar.gz")
run(`tar -xzf $file -C $perplexdir`)
end
end
## --- # # # # # # # # # # # # # Initial composition # # # # # # # # # # # # # #
## McDonough Pyrolite
#elements = [ "SIO2", "TIO2", "AL2O3", "FEO", "MNO", "MGO", "CAO", "NA2O", "K2O", "H2O", "CO2",]
#composition = [45.1242, 0.2005, 4.4623, 8.0723, 0.1354, 37.9043, 3.5598, 0.3610, 0.0291, 0.1511, 0.0440,]
## Kelemen (2014) primitive continental basalt. H2O and CO2 are guesses
#elements = [ "SIO2", "TIO2", "AL2O3", "FEO", "MNO", "MGO", "CAO", "NA2O", "K2O", "H2O", "CO2",]
#composition = [50.0956, 0.9564, 15.3224, 8.5103, 0.1659, 9.2520, 9.6912, 2.5472, 0.8588, 2.0000, 0.6000,]
# Kelemen (2014) primitive continental basalt excluding Mn and Ti since most melt models can"t handle them..
elements = [ "SIO2", "AL2O3", "FEO", "MGO", "CAO", "NA2O", "K2O", "H2O", "CO2",]
composition = [50.0956, 15.3224, 8.5103, 9.2520, 9.6912, 2.5472, 0.8588, 2.0000, 0.6000,]
## Average Archean basalt (EarthChem data)
#elements = [ "SIO2", "TIO2", "AL2O3", "FEO", "MNO", "MGO", "CAO", "NA2O", "K2O", "H2O", "CO2",]
#composition = [49.2054, 0.8401, 12.0551, 11.4018, 0.2198, 12.3997, 9.3113, 1.6549, 0.4630, 1.8935, 0.5555,]
9-element Array{Float64,1}: 50.0956 15.3224 8.5103 9.252 9.6912 2.5472 0.8588 2.0 0.6
## --- # # # # # # # # # # # Some solution model options # # # # # # # # # # # #
# Emphasis on phases from Green (2016) -- developed for metabasites, includes what is probably the best (and most expensive) amphibole model. Use with hp11ver.dat
G_solution_phases = "Augite(G)\nOpx(JH)\ncAmph(G)\noAmph(DP)\nO(JH)\nSp(JH)\nGrt(JH)\nfeldspar_B\nMica(W)\nBio(TCC)\nChl(W)\nCtd(W)\nCrd(W)\nSa(WP)\nSt(W)\nIlm(WPH)\nAtg(PN)\nT\nB\nF\nDo(HP)\nScap\nChum\nNeph(FB)\n"
G_excludes ="ged\nfanth\ngl\nilm\nilm_nol\n"
# Emphasis on phases from White (2014) -- developed for metapelites. Use with hp11ver.dat (though can apparenty run with hp02ver.dat without crashing)
W_solution_phases = "Omph(HP)\nOpx(W)\ncAmph(DP)\noAmph(DP)\nO(JH)\nSp(JH)\nGt(W)\nfeldspar_B\nMica(W)\nBi(W)\nChl(W)\nCtd(W)\nCrd(W)\nSa(WP)\nSt(W) \nIlm(WPH)\nAtg(PN)\nT\nB\nF\nDo(HP)\nScap\nChum\nPu(M)\n"
W_excludes = "andr\nts\nparg\ngl\nged\nfanth\n"
# Emphasis on phases from Jennings and Holland (2015) -- developed for mantle melting. Use with hp11ver.dat
JH_solution_phases = "Cpx(JH)\nOpx(JH)\ncAmph(DP)\noAmph(DP)\nO(JH)\nSp(JH)\nGrt(JH)\nfeldspar_B\nMica(W)\nBio(TCC)\nChl(W)\nCtd(W)\nCrd(W)\nSa(WP)\nSt(W)\nIlm(WPH)\nAtg(PN)\nT\nB\nF\nDo(HP)\nScap\nChum\nNeph(FB)\n"
JH_excludes = "ts\nparg\ngl\nged\nfanth\n"
# Emphasis on phases from Holland and Powell -- all phases can be used with hp02ver.dat.
HP_solution_phases = "Omph(HP)\nOpx(HP)\nGlTrTsPg\nAnth\nO(HP)\nSp(HP)\nGt(HP)\nfeldspar_B\nMica(CF)\nBio(TCC)\nChl(HP)\nCtd(HP)\nSapp(HP)\nSt(HP)\nIlHm(A)\nDo(HP)\nT\nB\nF\n"
HP_excludes = ""
""
## --- # # # # # # # # # # # # # Isobaric example # # # # # # # # # # # # # # # #
# Input parameters
P = 10000 # Pressure, bar
T_range = [500+273.15, 1500+273.15] # Temperature range, Kelvin
melt_model = "melt(G)"
# Configure (run build and vertex)
@time perplex_configure_isobar(perplexdir, scratchdir, composition, elements,
P, T_range, dataset="hp11ver.dat", npoints=100, excludes=G_excludes,
solution_phases=melt_model*"\n"*G_solution_phases)
## --- Query all properties at a single temperature -- results returned as text
T = 1450+273.15
data_isobaric = perplex_query_1d(perplexdir, scratchdir, T) |> print
54.792925 seconds (464.08 k allocations: 23.311 MiB, 0.02% gc time) ---------------------------------------- Stable phases at: T(K) = 1723.15 P(bar) = 10000.0 Phase Compositions (molar proportions): wt % vol % mol % mol H2O CO2 FEO MGO CAO NA2O K2O SIO2 AL2O3 melt(G) 99.35 98.28 97.33 0.590 0.18379 0.00000 0.20071 0.38896 0.29283 0.06964 0.01545 1.41272 0.25463 F 0.65 1.72 2.67 0.162E-01 0.15752 0.84248 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 Phase speciation (molar proportions): melt(G) foTL: 0.11411, faTL: 0.05889, abL: 0.16344, silL: 0.02543, kspL: 0.03626, woL: 0.17010, q8L: 0.04254 h2oL: 0.21569, oanL: 0.17354 F H2O: 0.15752, CO2: 0.84248 Molar Properties and Density: N(g) G(J) S(J/K) V(J/bar) Cp(J/K) Alpha(1/K) Beta(1/bar) Cp/Cv Density(kg/m3) melt(G) 166.44 -2819391 506.13 6.3017 256.47 0.87739E-04 0.40948E-05 1.0865 2641.2 F 39.92 -647615 214.51 4.0329 60.884 0.26245E-03 0.43779E-04 1.2189 989.75 System 98.88 -1674426 302.18 3.7844 152.35 0.90752E-04 0.47792E-05 1.0796 2612.7 System - fluid 98.23 -1663946 298.71 3.7192 151.37 0.87739E-04 0.40948E-05 1.0865 2641.2 Seismic Properties: Gruneisen_T Ks(bar) Mu(bar) V0(km/s) Vp(km/s) Vs(km/s) Poisson ratio melt(G) 0.57200 0.26533E+06 88443. 3.1695 3.8093 1.8299 0.35000 F 0.48403 27842. 0.0000 1.6772 1.6772 0.0000 0.50000 System 0.57049 0.24627E+06 43459. 3.0701 3.4123 1.2897 0.41667 System - fluid 0.57200 0.26533E+06 88443. 3.1695 3.8093 1.8299 0.35000 Bulk Composition: Complete Assemblage Solid+Melt Only mol g wt % mol/kg mol g wt % mol/kg H2O 0.111 2.000 2.023 1.123 0.108 1.954 1.989 1.104 CO2 0.014 0.600 0.607 0.138 0.000 0.000 0.000 0.000 FEO 0.118 8.510 8.607 1.198 0.118 8.510 8.664 1.206 MGO 0.230 9.252 9.357 2.322 0.230 9.252 9.419 2.337 CAO 0.173 9.691 9.801 1.748 0.173 9.691 9.866 1.759 NA2O 0.041 2.547 2.576 0.416 0.041 2.547 2.593 0.418 K2O 0.009 0.859 0.869 0.092 0.009 0.859 0.874 0.093 SIO2 0.834 50.096 50.664 8.432 0.834 50.096 50.997 8.488 AL2O3 0.150 15.322 15.496 1.520 0.150 15.322 15.598 1.530 Other Bulk Properties: Enthalpy (J/kg) = -.116682E+08 Specific Enthalpy (J/m3) = -.304862E+11 Entropy (J/K/kg) = 3056.09 Specific Entropy (J/K/m3) = 0.798480E+07 Heat Capacity (J/K/kg) = 1540.80 Specific Heat Capacity (J/K/m3) = 0.402573E+07 Solid Enthalpy (J/kg) = -.116992E+08 Solid Secific Enthalpy (J/m3) (2) = -.309002E+11 Solid Entropy (J/K/kg) = 3040.85 Solid Specific Entropy (J/K/m3) = 0.803158E+07 Solid Heat Capacity (J/K/kg) (1) = 1557.70 Solid Specific Heat Capacity (J/K/m3) (1) = 0.406988E+07 N.B.: Aggregate properties represent the entire stable assemblage. Solid aggregate properties represent solid and melt properties, but do not include molecular fluid properties. Chemical Potentials (J/mol): H2O CO2 FEO MGO CAO NA2O K2O SIO2 AL2O3 -509890. -673366. -478356. -738298. -864198. -0.102378E+07 -0.111639E+07 -0.106581E+07 -0.194481E+07 Variance (c-p+2) = 9 ----------------------------------------
## --- Query the full isobar -- results returned as dict
bulk = perplex_query_system(perplexdir, scratchdir) # Get system data for all temperatures. Set include_fluid = "n" to get solid+melt only
modes = perplex_query_modes(perplexdir, scratchdir) # || phase modes
melt = perplex_query_phase(perplexdir, scratchdir, melt_model) # || melt data
# Melt wt.% seems to be slightly inaccurate; use values from modes instead
melt["wt_pct"] = modes[melt_model]
# Create dictionary to hold solid composition and fill it using what we know from system and melt
solid = Dict()
solid["wt_pct"] = 100 .- melt["wt_pct"]
for e in ["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"]
solid[e] = (bulk[e] - (melt[e] .* melt["wt_pct"]/100)) ./ (solid["wt_pct"]/100)
end
renormalize!(solid,["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"],total=100)
┌ Warning: Perplex seems to be reporting mole fractions instead of weight percentages, attempting to correct └ @ StatGeochem /Users/cbkeller/Documents/julia/packages/StatGeochem.jl/src/utilities/Geochemistry.jl:649
## --- Plot melt composition as a function of melt percent
h = plot(xlabel="Percent melt", ylabel="Wt. % in melt", title="$melt_model + G_solution_phases, $P bar")
i = 0
for e in ["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"]
plot!(h, melt["wt_pct"], melt[e], label=e, color=lines[global i += 1])
plot!(h, melt["wt_pct"], bulk[e], label="", color=lines[i], linestyle=:dot)
end
plot!(h,fg_color_legend=:white, framestyle=:box)
# savefig(h,"MeltComposition.pdf")
display(h)
## --- Plot solid composition as a function of melt percent
h = plot(xlabel="Percent melt", ylabel="Wt. % in solid", title="$melt_model + G_solution_phases, $P bar")
i = 0
for e in ["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"]
plot!(h, melt["wt_pct"], solid[e], label=e, color=lines[global i +=1])
end
plot!(h,fg_color_legend=:white, framestyle=:box, legend=:topleft)
# savefig(h,"SolidComposition.pdf")
display(h)
## --- Plot modes of all phases as a function of temperature
h = plot(xlabel="T (C)", ylabel="Weight percent", title="$melt_model + G_solution_phases, $P bar")
for m in modes["elements"][3:end]
plot!(h, modes["T(K)"] .- 273.15, modes[m], label=m)
end
plot!(h,fg_color_legend=:white, framestyle=:box)
# savefig(h,"PhaseModes.pdf")
display(h)
## --- # # # # # # # # # # # Geothermal gradient example # # # # # # # # # # # #
# Input parameters
P_range = [280, 28000] # Pressure range to explore, bar (roughly 1-100 km depth)
T_surf = 273.15 # Temperature of surface (K)
geotherm = 0.1 # Geothermal gradient of 0.1 K/bar == about 28.4 K/km
melt_model = ""
# Configure (run build and vertex)
@time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
P_range, T_surf, geotherm, dataset="hp02ver.dat", excludes=HP_excludes,
solution_phases=HP_solution_phases, npoints=200, index=2)
# # Alternative configuration, using hpha02ver.dat
# @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
# P_range, T_surf, geotherm, dataset="hpha02ver.dat", excludes="qGL\n"*HP_excludes,
# solution_phases=HP_solution_phases, npoints=200, index=2)
# # Alternative configuration, using hpha02ver.dat and new phases for metapelites
# @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
# P_range, T_surf, geotherm, dataset="hpha02ver.dat", excludes="qGL\n"*W_excludes,
# solution_phases=W_solution_phases, npoints=200, index=2)
MethodError: no method matching perplex_configure_geotherm(::String, ::String, ::Array{Float64,1}, ::Array{String,1}, ::Array{Int64,1}, ::Float64, ::Float64; dataset="hp02ver.dat", excludes="", solution_phases="Omph(HP)\nOpx(HP)\nGlTrTsPg\nAnth\nO(HP)\nSp(HP)\nGt(HP)\nfeldspar_B\nMica(CF)\nBio(TCC)\nChl(HP)\nCtd(HP)\nSapp(HP)\nSt(HP)\nIlHm(A)\nDo(HP)\nT\nB\nF\n", npoints=200, index=2) Closest candidates are: perplex_configure_geotherm(::String, ::String, ::Array{#s47,N} where N where #s47<:Number; elements, P_range, T_surf, geotherm, dataset, solution_phases, excludes, index, npoints) at /Users/cbkeller/Documents/julia/packages/StatGeochem.jl/src/utilities/Geochemistry.jl:427 Stacktrace: [1] top-level scope at util.jl:156 [2] top-level scope at In[10]:8
## --- Plot modes of all phases as a function of temperature
# Get phase modes
modes = perplex_query_modes(perplexdir, scratchdir, index=2)
h = plot(xlabel="T (C)", ylabel="Weight percent")
for m in modes["elements"][3:end]
plot!(h, modes["T(K)"] .- 273.15, modes[m], label=m)
end
plot!(h,fg_color_legend=:white, framestyle=:box)
# savefig(h,"GeothermPhaseModes.pdf")
display(h)
## --- Plot seismic properties
# Query seismic properties along the whole profile
seismic = perplex_query_seismic(perplexdir, scratchdir, index=2)
seismic["vp/vs"][seismic["vp/vs"] .> 100] .= NaN # Exclude cases where vs drops to zero
h = plot(xlabel="Pressure", ylabel="Property")
plot!(h,seismic["P(bar)"],seismic["vp,km/s"], label="vp,km/s")
plot!(h,seismic["P(bar)"],seismic["vp/vs"], label="vp/vs")
plot!(h,seismic["P(bar)"],seismic["rho,kg/m3"]/1000, label="rho, g/cc")
plot!(h,seismic["P(bar)"],seismic["T(K)"]/1000, label="T(K)/1000")
# savefig(h,"GeothermSeismicProperties.pdf")
display(h)