NOTE: This notebook contains an interactive figure with sliders. it relies on python modules ipympl and mpl_interactions.
If you want to run this notebook with static figures, restart the kernel, clear any history and set interactive = False in the first code block.
In the mantle, it is common to assume that convecting material is at chemical equilibrium; all of the reactions between phases keep pace with the changes in pressure and temperature. Because of this relaxation, physical properties such as heat capacity $C_P$, thermal expansion $\alpha$ and compressibility $\beta$ must be computed by numerical differentiation of the entropy $\mathcal{S}$ and volume $\mathcal{V}$. It is these values, rather than the unrelaxed values output as standard by BurnMan and PerpleX which should be used in geodynamic simulations.
Relaxed properties can sometimes be very different from their unrelaxed counterparts. Take, for example, the univariant reaction forsterite -> Mg-wadsleyite. These transformation involves a step change in volume, and thus the relaxed compressibility at the transition is infinite. Obviously, if geodynamics software uses compressibility as an input parameter, then whichever meshing is chosen, it will completely miss the transition. There are two solutions to this problem:
The second method is used here to create 1D material property profiles which can be directly used by $ASPECT$. The user of this notebook can vary important mineral physics parameters (rock type, potential temperature, surface gravity) and smoothing parameters (Gaussian widths).
First, let's install some modules that we need to run interactive figures.
interactive = True
if interactive:
%matplotlib ipympl
import ipywidgets as widgets
import mpl_interactions.ipyplot as iplt
Now a few imports
import matplotlib.pyplot as plt
import numpy as np
import burnman
from burnman import Layer
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, brentq
from scipy.integrate import odeint
from scipy.interpolate import UnivariateSpline
plt.style.use('bmh')
The following code block sets up some parameters for our problem, including the PerpleX model we will use to determine physical properties of the planet.
perplex_filename = '../../burnman/data/input_perplex/in23_1.tab' # 'example23_hires.tab' # '../../burnman/data/input_perplex/in23_1.tab'
potential_temperature = 1550.
outer_radius = 6371.e3
thickness = 550.e3
n_points = 251
pressure_top = 1.e5
gravity_bottom = 10.
depths = np.linspace(thickness, 0., n_points)
rock = burnman.PerplexMaterial(perplex_filename)
layer = Layer(name='Mantle', radii=outer_radius-depths)
layer.set_material(rock)
layer.set_temperature_mode(temperature_mode='adiabatic',
temperature_top=1550.)
layer.set_pressure_mode(pressure_mode='self-consistent',
pressure_top=1.e5,
gravity_bottom=gravity_bottom)
layer.make()
truncate = 4. # truncates the convolution Gaussian at 4 sigma
The following code block reinterpolates the entropy and volume onto a grid for smoothing purposes. As we're using a PerpleX table, we could just have read the original grid in directly, but this way, we could more easily adapt the code to use other BurnMan materials. This step takes a few seconds (10--30 s on a fast ultrabook), but we only do it once.
n_gridpoints = (501, 101)
min_grid_pressure = rock.bounds[0][0]
max_grid_pressure = rock.bounds[0][1]
min_grid_temperature = rock.bounds[1][0]
max_grid_temperature = rock.bounds[1][1]
grid_pressures = np.linspace(min_grid_pressure, max_grid_pressure, n_gridpoints[0])
grid_temperatures = np.linspace(min_grid_temperature, max_grid_temperature, n_gridpoints[1])
pp, TT = np.meshgrid(grid_pressures, grid_temperatures)
mesh_shape = pp.shape
pp = np.ndarray.flatten(pp)
TT = np.ndarray.flatten(TT)
grid_entropies = np.zeros_like(pp)
grid_volumes = np.zeros_like(pp)
grid_entropies, grid_volumes = layer.material.evaluate(['S', 'V'], pp, TT)
grid_entropies = grid_entropies.reshape(mesh_shape)
grid_volumes = grid_volumes.reshape(mesh_shape)
This long code block sets up three functions, which:
# Define function to find an isentrope given a
# 2D entropy interpolation function
# Here we use fsolve, because we'll normally have a good starting guess
# from the previous pressure
def interp_isentrope(interp, pressures, entropies, T_guess):
def _deltaS(args, S, P):
T = args[0]
return interp(P, T)[0] - S
sol = [T_guess]
temperatures = np.empty_like(pressures)
for i in range(len(pressures)):
sol = fsolve(_deltaS, sol, args=(entropies[i], pressures[i]))
temperatures[i] = sol[0]
return temperatures
def relaxed_profile(layer, pressure_stdev, temperature_stdev,
truncate):
unsmoothed_T_spline = UnivariateSpline(layer.pressure[::-1], layer.temperature[::-1])
unsmoothed_grid_isentrope_temperatures = unsmoothed_T_spline(grid_pressures)
# Having defined the grid and calculated unsmoothed properties,
# we now calculate the smoothed entropy and volume and derivatives with
# respect to pressure and temperature.
S_interps = burnman.tools.math.interp_smoothed_array_and_derivatives(array=grid_entropies,
x_values=grid_pressures,
y_values=grid_temperatures,
x_stdev=pressure_stdev,
y_stdev=temperature_stdev,
truncate=truncate)
interp_smoothed_S, interp_smoothed_dSdP, interp_smoothed_dSdT = S_interps
V_interps = burnman.tools.math.interp_smoothed_array_and_derivatives(array=grid_volumes,
x_values=grid_pressures,
y_values=grid_temperatures,
x_stdev=pressure_stdev,
y_stdev=temperature_stdev,
truncate=truncate)
interp_smoothed_V, interp_smoothed_dVdP, interp_smoothed_dVdT = V_interps
# Now we can calculate and plot the relaxed and smoothed properties along the isentrope
smoothed_temperatures = interp_isentrope(interp_smoothed_S, layer.pressure[::-1], layer.S[::-1], layer.temperature[-1])[::-1]
densities = layer.material.evaluate(['rho'], layer.pressure, smoothed_temperatures)[0]
volumes = np.array([interp_smoothed_V(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])
dSdT = np.array([interp_smoothed_dSdT(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])
dVdT = np.array([interp_smoothed_dVdT(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])
dVdP = np.array([interp_smoothed_dVdP(p, T)[0] for (p, T) in zip(*[layer.pressure, smoothed_temperatures])])
alphas_relaxed = dVdT / volumes
compressibilities_relaxed = -dVdP / volumes
specific_heats_relaxed = smoothed_temperatures * dSdT / (densities[0]*volumes[0])
dT = 0.1
Vpsub, Vssub = layer.material.evaluate(['p_wave_velocity', 'shear_wave_velocity'],
layer.pressure, smoothed_temperatures-dT/2.)
Vpadd, Vsadd = layer.material.evaluate(['p_wave_velocity', 'shear_wave_velocity'],
layer.pressure, smoothed_temperatures+dT/2.)
Vps = (Vpadd + Vpsub)/2.
Vss = (Vsadd + Vssub)/2.
dVpdT = (Vpadd - Vpsub)/dT
dVsdT = (Vsadd - Vssub)/dT
depths = layer.outer_radius - layer.radii
return (smoothed_temperatures, layer.pressure, depths, layer.gravity, densities,
alphas_relaxed, compressibilities_relaxed, specific_heats_relaxed,
Vss, Vps, dVsdT, dVpdT)
flag_index = {'T': 0, 'P': 1, 'z': 2, 'g': 3, 'rho': 4,
'alpha': 5, 'beta_T': 6, 'Cp': 7,
'Vs': 8, 'Vp': 9, 'dVsdT': 10,
'dVpdT': 11}
def save_relaxed_properties(layer, P_GPa_gaussian, T_K_gaussian, outfile='isentrope_properties.txt'):
"""
A function to output smoothed, relaxed properties for use in ASPECT
depth, pressure, temperature, density, gravity, Cp (per kilo), thermal expansivity
"""
d = relaxed_profile(layer, P_GPa_gaussian*1.e9, T_K_gaussian, truncate)[::-1]
np.savetxt(outfile, X=np.array([d[flag_index['z']],
d[flag_index['P']],
d[flag_index['T']],
d[flag_index['rho']],
d[flag_index['g']],
d[flag_index['alpha']],
d[flag_index['Cp']],
d[flag_index['beta_T']],
d[flag_index['Vs']],
d[flag_index['Vp']],
d[flag_index['dVsdT']],
d[flag_index['dVpdT']]]).T,
header=('# This ASPECT-compatible file contains material '
'properties calculated along an isentrope by the '
f'BurnMan software.\n# POINTS: {n_points}\n'
'# depth (m), pressure (Pa), temperature (K), '
'density (kg/m^3), gravity (m/s^2), '
'thermal expansivity (1/K), specific heat (J/K/kg), '
'compressibility (1/Pa), seismic Vs (m/s), '
'seismic Vp (m/s), seismic dVs/dT (m/s/K), '
'seismic dVp/dT (m/s/K)\n'
'depth pressure '
'temperature density '
'gravity thermal_expansivity '
'specific_heat compressibility '
'seismic_Vs seismic_Vp '
'seismic_dVs_dT seismic_dVp_dT'),
fmt='%.10e', delimiter='\t', comments='')
print('File saved to {0}'.format(outfile))
Save file
save_relaxed_properties(layer, P_GPa_gaussian=0.25, T_K_gaussian=0.25)
The following code block attempts to make updating the multipart interactive figure as efficient as possible. It does this by calculating all of the properties in an inner function, and storing them in a global parameter (global_stored_properties). That function is then wrapped in an outer function that dictates which properties to return. If the inner function has been previously called with the same input parameters, the previously stored properties are returned, otherwise, all properties are calculated from the new input parameters.
global global_PT_smooth
global_PT_smooth = [None, None]
global_stored_properties = [None]
def plot_y(flag):
index = flag_index[flag]
def f(x, P_GPa_gaussian, T_K_gaussian):
if P_GPa_gaussian == global_PT_smooth[0] and T_K_gaussian == global_PT_smooth[1]:
pass
else:
global_PT_smooth[0] = P_GPa_gaussian
global_PT_smooth[1] = T_K_gaussian
f = relaxed_profile(layer, P_GPa_gaussian*1.e9, T_K_gaussian,
truncate)
global_stored_properties[0] = f
return global_stored_properties[0][index]
return f
plt.rcParams['figure.figsize'] = 8, 5 # inches
fig = plt.figure()
px, py = [2, 3]
depths = layer.outer_radius - layer.radii
gravity = layer.gravity
x = depths/1.e3
xlabel = 'Depths (km)'
ax_T = fig.add_subplot(px, py, 1)
ax_T.plot(x, layer.temperatures, label='unrelaxed')
ax_T.set_ylabel('Temperature (K)')
ax_T.set_xlabel(xlabel)
ax_g = fig.add_subplot(px, py, 2)
ax_g.plot(x, layer.gravity)
ax_g.set_ylabel('Gravity (m/s^2)')
ax_g.set_xlabel(xlabel)
ax_rho = fig.add_subplot(px, py, 3)
ax_rho.plot(x, layer.rho, label='$\rho$ (kg/m$^3$)')
ax_rho.plot(x, layer.v_p, label='P (km/s)')
ax_rho.plot(x, layer.v_s, label='S (km/s)')
ax_rho.set_ylabel('Densities/Velocities')
ax_rho.set_xlabel(xlabel)
ax_alpha = fig.add_subplot(px, py, 4)
ax_alpha.plot(x, layer.alpha)
ax_alpha.set_ylabel('alpha (/K)')
ax_alpha.set_xlabel(xlabel)
ax_beta = fig.add_subplot(px, py, 5)
ax_beta.plot(x, layer.beta_T)
ax_beta.set_ylabel('compressibilities (/Pa)')
ax_beta.set_xlabel(xlabel)
ax_cp = fig.add_subplot(px, py, 6)
ax_cp.plot(x, layer.C_p/layer.molar_mass)
ax_cp.set_ylabel('Cp (J/K/kg)')
ax_cp.set_xlabel(xlabel)
# Relaxed, unsmoothed properties
ax_T.plot(x, plot_y('T')(x, 0., 0.), label='relaxed, unsmoothed')
ax_g.plot(x, plot_y('g')(x, 0., 0.))
ax_alpha.plot(x, plot_y('alpha')(x, 0., 0.))
ax_beta.plot(x, plot_y('beta_T')(x, 0., 0.))
ax_cp.plot(x, plot_y('Cp')(x, 0., 0.))
if interactive:
# Interactive smoothing
P_GPa_gaussian = np.linspace(0., 3., 51)
T_K_gaussian = np.linspace(0., 30, 41)
controls = iplt.plot(x, plot_y('T'), P_GPa_gaussian=P_GPa_gaussian, T_K_gaussian=T_K_gaussian, ax=ax_T, label='relaxed, smoothed')
_ = iplt.plot(x, plot_y('g'), controls=controls, ax=ax_g)
_ = iplt.plot(x, plot_y('alpha'), controls=controls, ax=ax_alpha)
_ = iplt.plot(x, plot_y('beta_T'), controls=controls, ax=ax_beta)
_ = iplt.plot(x, plot_y('Cp'), controls=controls, ax=ax_cp)
else:
# Non-interactive smoothing
P_GPa_gaussian = 0.5
T_K_gaussian = 20.
ax_T.plot(x, plot_y('T')(x, P_GPa_gaussian, T_K_gaussian), label='relaxed, smoothed')
ax_g.plot(x, plot_y('g')(x, P_GPa_gaussian, T_K_gaussian))
ax_alpha.plot(x, plot_y('alpha')(x, P_GPa_gaussian, T_K_gaussian))
ax_beta.plot(x, plot_y('beta_T')(x, P_GPa_gaussian, T_K_gaussian))
ax_cp.plot(x, plot_y('Cp')(x, P_GPa_gaussian, T_K_gaussian))
ax_T.legend(loc='lower right',prop={'size':8})
fig.set_tight_layout(True)