BurnMan - Compute seismic velocities - advanced case

Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S. Released under GPL v2 or later.

This example script is an extension the simple case. It gives suggestions on how to change the geotherm and composition, and how to compare to the seismic model with depth instead of pressures. One can choose to run the different cells for the various options (press shift+enter on the cell of choice).

requires:

  • geotherms

  • seismic models

  • compute seismic velocities

  • composite mineral helpers

teaches:

  • changing composition

  • changing geotherm

  • changing averaging scheme

  • computing self-consistent depth to compare to seismology

1. Import modules

Import of relevant modules.


In [1]:
import os, sys, numpy as np, matplotlib.pyplot as plt
sys.path.insert(1,os.path.abspath('..'))
import burnman
from burnman import minerals

2. Import seismic model

Here we load the PREM seismic velocity model for comparison.


In [2]:
seismic_model = burnman.seismic.PREM()

We create an array of 20 depths at which we want to evaluate PREM, and then query the seismic model for the pressure, density, P wave speed, S wave speed, and bulk sound velocity at those depths


In [3]:
depths = np.linspace(750e3, 2800e3, 20)
pressure, seis_rho, seis_vp, seis_vs, seis_vphi = seismic_model.evaluate_all_at(depths)
plt.plot(pressure/1.e9,seis_vs/1.e3,'k',label='Vs')
plt.plot(pressure/1.e9,seis_vp/1.e3,'b',label='Vp')
plt.plot(pressure/1.e9,seis_rho/1.e3,'r',label='Vphi')
plt.plot(pressure/1.e9,seis_vphi/1.e3,'g',label='rho')
plt.xlabel('pressure (GPa)')
plt.ylabel('velocity (km/s) density (kg/m^3)')
plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)
plt.title('PREM')
plt.legend()
plt.show();


3. Input composition

Beneath are three different methods to input a composition. Choose and run one of the cells. Here we show how to print the minerals in a library (needs improvement) if one wants to choose a different one.

Available mineral libraries:

  • SLB_2005
  • SLB_2011
  • Matas_etal_2007
  • Murakami_etal_2012

In [4]:
minlib=minerals.SLB_2011
print dir(minlib)


['AsymmetricRegularSolution', 'Fraction', 'IdealSolution', 'Mineral', 'SolidSolution', 'SolutionModel', 'SymmetricRegularSolution', '__builtins__', '__doc__', '__file__', '__name__', '__package__', 'ab', 'akimotoite', 'al', 'al_perovskite', 'al_post_perovskite', 'albite', 'almandine', 'alpv', 'an', 'anorthite', 'appv', 'atomic_masses', 'burnman', 'c2c', 'c2c_pyroxene', 'ca_ferrite_structured_phase', 'ca_perovskite', 'ca_tschermaks', 'capv', 'cats', 'cen', 'cf', 'clinoenstatite', 'clinopyroxene', 'co', 'coes', 'coesite', 'compositional_array', 'constants', 'corundum', 'cpx', 'di', 'dictionarize_formula', 'dictionarize_site_formula', 'diopside', 'en', 'enstatite', 'fa', 'fayalite', 'fe_akimotoite', 'fe_bridgmanite', 'fe_ca_ferrite', 'fe_perovskite', 'fe_post_perovskite', 'fe_ringwoodite', 'fe_wadsleyite', 'fec2', 'fecf', 'feil', 'fepv', 'feri', 'ferropericlase', 'ferrosilite', 'fewa', 'fo', 'formula_mass', 'forsterite', 'fppv', 'fs', 'garnet', 'gr', 'grossular', 'gt', 'hc', 'he', 'hedenbergite', 'hercynite', 'hp_clinoenstatite', 'hp_clinoferrosilite', 'hpcen', 'hpcfs', 'il', 'ilmenite_group', 'jadeite', 'jd', 'jd_majorite', 'jdmj', 'kd', 'ky', 'kyanite', 'magnesiowuestite', 'mg_akimotoite', 'mg_bridgmanite', 'mg_ca_ferrite', 'mg_fe_aluminous_spinel', 'mg_fe_bridgmanite', 'mg_fe_olivine', 'mg_fe_perovskite', 'mg_fe_ringwoodite', 'mg_fe_silicate_perovskite', 'mg_fe_wadsleyite', 'mg_majorite', 'mg_perovskite', 'mg_post_perovskite', 'mg_ringwoodite', 'mg_tschermaks', 'mg_wadsleyite', 'mgc2', 'mgcf', 'mgil', 'mgmj', 'mgpv', 'mgri', 'mgts', 'mgwa', 'mppv', 'mw', 'na_ca_ferrite', 'nacf', 'neph', 'nepheline', 'np', 'odi', 'ol', 'opx', 'ordered_compositional_array', 'ortho_diopside', 'orthopyroxene', 'pe', 'periclase', 'pkgutil', 'plag', 'plagioclase', 'post_perovskite', 'ppv', 'process_solution_chemistry', 'pv', 'py', 'pyrope', 're', 'read_masses', 'ri', 'seif', 'seifertite', 'sp', 'spinel', 'spinel_group', 'spinelloid_III', 'st', 'stishovite', 'wa', 'warnings', 'wu', 'wuestite']
Composite I. Fixed minerals

We define composite object and name it "rock". A composite is made by giving burnman.composite a list of minerals and their molar fractions. Here "rock" has two constituent minerals: it is 80% Mg perovskite and 20% periclase. More than two minerals can simply be added to the composite.


In [5]:
amount_of_perovskite=0.8
rock = burnman.Composite([amount_of_perovskite, 1.0-amount_of_perovskite],
    [minerals.SLB_2011.mg_perovskite(), minerals.SLB_2011.periclase()])
Composite II. Including iron

In [6]:
amount_perovskite = 0.95
fe_num_perovskite = 0.08
fe_num_periclase  = 0.2
perovskite=minerals.SLB_2011.mg_fe_perovskite()
ferropericlase=minerals.SLB_2011.ferropericlase()
perovskite.set_composition([1.0 - fe_num_perovskite, fe_num_perovskite, 0.0])
ferropericlase.set_composition([1.0 - fe_num_periclase, fe_num_periclase])
rock = burnman.Composite([amount_perovskite, 1.0-amount_perovskite],
    [perovskite, ferropericlase])
Composite III. Weight percentages

Here is the implementations for weight percentages, where the phase fractions are calculated. The partition coefficient of iron is allowed to vary with depth.


In [7]:
weight_percents = {'Mg':0.38, 'Fe': 0.16, 'Si':0.46, 'Ca':0., 'Al':0.}
Kd_0 = .59 #Fig 5 Nakajima et al 2012

phase_fractions,relative_molar_percent = burnman.\
        					calculate_phase_percents(weight_percents)
iron_content = lambda p,t: burnman.calculate_partition_coefficient\
        					(p,t,relative_molar_percent,Kd_0)

perovskite.set_state(1.e9, 1000.)
def pv_fper_pt_dependent(iron_content):
    pv_fe, fper_fe = iron_content(perovskite.pressure,perovskite.temperature)
    perovskite.set_composition([1.0-pv_fe, pv_fe, 0.0])
    ferropericlase.set_composition([1.0-fper_fe, fper_fe])
    return perovskite, ferropericlase

rock = burnman.Composite([phase_fractions['pv'], phase_fractions['fp']],
    [pv_fper_pt_dependent(iron_content)[0], pv_fper_pt_dependent(iron_content)[1]])
Set method of composite

At this point we want to tell the rock which equation of state to use for its thermoelastic calculations. Methods available are:

  • bm2 -- second order method (Stixrude & Lithgow-Berteloni 2005)- no temperature effects
  • bm3 -- third order method (Stixrude & Lithgow-Berteloni 2005) - no temperature effects
  • slb2 -- second order method (Stixrude & Lithgow-Berteloni 2005)
  • slb3 -- third order method (Stixrude & Lithgow-Berteloni 2005)
  • mgd2 -- second order method (Matas et al. 2007) non-benchmarked
  • mgd3 -- third order method (Matas et al. 2007) non-benchmarked

In general, we recommend the 'slb3' equation of state as the most self-consistent model. The parameters from the SLB_2011 mineral library are fit using this model and we suggest not changing this method when using this library.


In [8]:
rock.set_method('slb3')

4. Input temperature

Now we get an array of temperatures at which will be used for computing the seismic properties of the rock. Here one can choose between the Brown and Shankland (1981) model or an adiabat.

Geotherm I. Brown & Shankland


In [9]:
temperature = burnman.geotherm.brown_shankland(pressure)
title='Brown&Shankland'
Geotherm II. Adiabat

In [10]:
T0=1900 #K anchor temperature
temperature = burnman.geotherm.adiabatic(pressure, T0, rock)
title='adiabat'

In [11]:
# plot temperature
plt.plot(pressure/1.e9,temperature,'r')
plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)
plt.title(title)
plt.xlabel('pressure (GPa)')
plt.ylabel('temperature (K)')
plt.show();


4. Calculate velocities

Here is the step which does the heavy lifting. burnman.velocities_from_rock sets the state of the rock at each of the pressures and temperatures defined,then calculates the elastic moduli and density of each individual phase. After that,it performs elastic averaging on the phases to get a single bulk and shear modulus for the rock. This averaging scheme defaults to Voigt-Reuss-Hill, but see example_averaging.py for other options. Finally, it calculates the seismic wave speeds for the whole rock. It returns a tuple of density, p-wave velocity s-wave velocity, bulk sound speed, bulk modulus, and shear modulus.


In [12]:
density, vp, vs, vphi, K, G = burnman.velocities_from_rock(rock, pressure, temperature, burnman.averaging_schemes.VoigtReussHill())

5. Plot results

All the work is done except the plotting! Here we want to plot the seismic wave speeds and the density against PREM using the matplotlib plotting tools. First we plot the results with pressure and then convert to depths to plot with depth.

Plots I. As a function of pressure

In [13]:
# First, we plot the s-wave speed verses the PREM s-wave speed
plt.figure(figsize=(16,4.5))
plt.subplot(1,3,1)
plt.plot(pressure/1.e9,vs/1.e3,color='b',linestyle='-',marker='o', markerfacecolor='b',markersize=4,label='computation')
plt.plot(pressure/1.e9,seis_vs/1.e3,color='k',linestyle='-',marker='o', markerfacecolor='k',markersize=4,label='reference')
plt.title("S wave speed (km/s)")
plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)
plt.xlabel('pressure (GPa)')
plt.legend(loc='lower right')
plt.ylim(5,8.0)
  

# Next, we plot the p-wave speed verses the PREM p-wave speed
plt.subplot(1,3,2)
plt.plot(pressure/1.e9,vp/1.e3,color='b',linestyle='-',marker='o',markerfacecolor='b',markersize=4)
plt.plot(pressure/1.e9,seis_vp/1.e3,color='k',linestyle='-',marker='o',markerfacecolor='k',markersize=4)
plt.title("P wave speed (km/s)")
plt.xlabel('pressure (GPa)')
plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)
plt.ylim(10,16)
    
# Next, we plot the density verses the PREM density
plt.subplot(1,3,3)
plt.plot(pressure/1.e9,density/1.e3,color='b',linestyle='-',marker='o', markerfacecolor='b',markersize=4)
plt.plot(pressure/1.e9,seis_rho/1.e3,color='k',linestyle='-',marker='o', markerfacecolor='k',markersize=4)
plt.xlim(min(pressure)/1.e9,max(pressure)/1.e9)
plt.xlabel("Pressure (GPa)")

plt.title("density (kg/m^3)");


Plot II. As a function of depth

In [14]:
# calculate the self-consistent depth at these pressures
depths2=burnman.depths_for_rock(rock,pressure,temperature)
print 'The maximum shift in depths is', round(max(abs( depths2-depths))/1.e3), ' km.'
print 'The closer the composition will lie to PREM, the less this shift will be.'
plt.figure(figsize=(16,4.5))
# First, we plot the s-wave speed verses the PREM s-wave speed
plt.subplot(1,3,1)
plt.plot(depths/1.e3,vs/1.e3,color='b',linestyle='-',marker='o', markerfacecolor='b',markersize=4,label='computation')
plt.plot(depths2/1.e3,seis_vs/1.e3,color='k',linestyle='-',marker='o', markerfacecolor='k',markersize=4,label='reference')
plt.title("S wave speed (km/s)")
plt.xlim(min(depths)/1.e3,max(depths)/1.e3)
plt.legend(loc='lower right')
plt.ylim(5,8.0)
    
# Next, we plot the p-wave speed verses the PREM p-wave speed
plt.subplot(1,3,2)
plt.plot(depths/1.e3,vp/1.e3,color='b',linestyle='-',marker='o',markerfacecolor='b',markersize=4)
plt.plot(depths2/1.e3,seis_vp/1.e3,color='k',linestyle='-',marker='o',markerfacecolor='k',markersize=4)
plt.title("P wave speed (km/s)")
plt.xlim(min(depths)/1.e3,max(depths)/1.e3)
plt.ylim(10,16)
    
# Next, we plot the density verses the PREM density
plt.subplot(1,3,3)
plt.plot(depths/1.e3,density/1.e3,color='b',linestyle='-',marker='o', markerfacecolor='b',markersize=4)
plt.plot(depths2/1.e3,seis_rho/1.e3,color='k',linestyle='-',marker='o', markerfacecolor='k',markersize=4)
plt.xlim(min(depths)/1.e3,max(depths)/1.e3)
plt.xlabel("Depth (km)")
plt.title("density (kg/m^3)");


The maximum shift in depths is 12.0  km.
The closer the composition will lie to PREM, the less this shift will be.

In [ ]: