Julia-Perple_X interface demo

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.

Launch Binder notebook

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 required Julia packages

In [1]:
## --- Load (and install if neccesary) the StatGeochem package which has the resampling functions we'll want

try
    using StatGeochem
catch
    using Pkg
    Pkg.add(PackageSpec(url="https://github.com/brenhinkeller/StatGeochem.jl"))
    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

Try to download and install Perple_X

In [2]:
## --- 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

Configure Perple_X options

In [3]:
## --- # # # # # # # # # # # # # 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,]
Out[3]:
9-element Array{Float64,1}:
 50.0956
 15.3224
  8.5103
  9.252 
  9.6912
  2.5472
  0.8588
  2.0   
  0.6   
In [4]:
## --- # # # # # # # # # # # 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 = ""
Out[4]:
""

Run Perple_X

In [5]:
## --- # # # # # # # # # # # # # 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


----------------------------------------

In [6]:
## --- 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 results

In [7]:
## --- 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)
0 25 50 75 100 0 10 20 30 40 50 60 melt(G) + G_solution_phases, 10000 bar Percent melt Wt. % in melt SIO2 AL2O3 FEO MGO CAO NA2O K2O
In [8]:
## --- 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)
0 25 50 75 100 0 10 20 30 40 50 melt(G) + G_solution_phases, 10000 bar Percent melt Wt. % in solid SIO2 AL2O3 FEO MGO CAO NA2O K2O
In [9]:
## --- 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)
600 800 1000 1200 1400 0 25 50 75 100 melt(G) + G_solution_phases, 10000 bar T (C) Weight percent Augite(G) cAmph(G) Chl(W) Mica(W) Bio(TCC) Do(HP) zo ab q F feldspar_B Opx(JH) melt(G)
In [10]:
## --- # # # # # # # # # # # 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
In [11]:
## --- 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)
0 500 1000 1500 2000 2500 0 10 20 30 40 50 T (C) Weight percent Bio(TCC) Ctd(HP) Gt(HP) GlTrTsPg sid acti q8L Do(HP) cz mic ky ab Mica(CF) q Anth zo feldspar_B cc Chl(HP) F Opx(HP) Omph(HP) O(HP) abL faGL kspL anL diL enL anl qGL
In [12]:
## --- 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)
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2.5×10 4 1 2 3 4 5 6 7 Pressure Property vp,km/s vp/vs rho, g/cc T(K)/1000
In [ ]: