What function can we reasonably use in order to parameterise the upper mantle viscosity (and test Theia)
In [3]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
In [4]:
def sinha_butler(nu_j, lam, z_j, z):
"""Smoothed jump
This is the function used by Sinha and Butler
(2007; JGR 112:B10406 http://dx.doi.org/10.1029/2006JB004850)
"""
z_j = z_j / np.max(z)
z_norm = z / np.max(z)
nu = ((nu_j-1)/2)*np.tanh(lam*(z_j-z_norm))+((nu_j+1)/2)
return nu
This can be used to give something that looks quite like the function used by Sinha and Butler
In [5]:
z = np.linspace(1.0, 0.0, 1000)
lam = 50.0
z_j = 0.77
nu_1 = sinha_butler(1, lam, z_j, z)
nu_10 = sinha_butler(10, lam, z_j, z)
nu_100 = sinha_butler(100, lam, z_j, z)
plt.plot(nu_1, z, nu_10, z, nu_100, z)
Out[5]:
And, because we are normalising inside the fucntion, we can use real numbers. (However, the lambda parameter is a bit unclear)
In [6]:
z = np.linspace(6346.0, 3480.0, 1000)
lam = 50.0
z_j = 6346.0 - 670.0
nu_1 = sinha_butler(1, lam, z_j, z)
nu_10 = sinha_butler(10, lam, z_j, z)
nu_100 = sinha_butler(100, lam, z_j, z)
plt.plot(nu_1, z, nu_10, z, nu_100, z)
Out[6]:
Play around with hypobolic tangent to see what it does
In [7]:
plt.plot(0.5*np.tanh(np.linspace(-2*np.pi, 2*np.pi)), np.linspace(-2*np.pi, 2*np.pi),
0.5*np.tanh(10*(np.linspace(-2*np.pi, 2*np.pi)-2)), np.linspace(-2*np.pi, 2*np.pi))
Out[7]:
In [8]:
plt.plot(sinha_butler(0.1, lam, z_j, z), z)
Out[8]:
As long as the tapers fall off quickly, we can create a low viscosity zone
In [9]:
def sinh_notch(z, z_j1, nu_j1, lam_j1, z_j2, nu_j2, lam_j2):
z_j1 = z_j1 / np.max(z)
z_j2 = z_j2 / np.max(z)
z_norm = z / np.max(z)
nu_j2 = (nu_j2/nu_j1)
nu = ((((nu_j1-1)/2)*np.tanh(lam_j1*(z_j1-z_norm))+((nu_j1+1)/2)) *
(((nu_j2-1)/2)*np.tanh(lam_j2*(z_j2-z_norm))+((nu_j2+1)/2)))
return nu
z = np.linspace(6346.0, 3480.0, 1000)
lam_j1 = 1000
z_j1 = 6346.0 - 100.0
nu_j1 = 0.05
z = np.linspace(6346.0, 3480.0, 1000)
lam_j2 = 1000
z_j2 = 6346.0 - 670.0
nu_j2 = 11
plt.plot(sinh_notch(z, z_j1, nu_j1, lam_j1, z_j2, nu_j2, lam_j2), z,
sinh_notch(z, z_j1, 0.05, 1000, z_j2, 12, 100),z,
sinh_notch(z, z_j1, 0.05, 1000, z_j2, 13, 50),z,)
plt.xscale('log')
But this isn't great if the lower jump is too smooth
In [10]:
plt.plot(sinh_notch(z, z_j1, 0.05, 1000, z_j2, 20, 20),z)
Out[10]:
How about just doing a smooth lower boundary, and making the top arbitary:
In [69]:
def visc_profile(nu_lid, z_lid, nu_j, lam, z_j, z):
"""Smoothed lower viscosity jump for base of the asthenosphere
The upper boundary is defined by a lid nu_lid times more
viscous than the asthenosphere with a base at depth z_lid
For the lower boundary we use the function of Sinha and Butler
(2007; JGR 112:B10406 http://dx.doi.org/10.1029/2006JB004850)
which is parameterised by:
* z_j: the depth (in km) of the viscosty jump at the "base"
of the asthernosphere
* lam: a parameter giving the width of the jump (numbers
between 20 and 200 give good values), larger is smoother
* nu_j: the radial viscosity factor: how much more viscous is the
mantle below the jump, compared to the manlte above the jump.
"""
z_j = z_j / np.max(z)
z_norm = z / np.max(z)
nu = ((1-nu_j)/2)*np.tanh(lam*(z_j-z_norm))+((nu_j+1)/2)
for i in range(len(nu)):
if z[i] < z_lid:
nu[i] = nu_lid
nu[-1] = nu_j
return nu
In [70]:
z = np.linspace(0.0, 3480.0, 1000)
z = np.linspace(0.0, 1000.0, 1000)
nu_1 = visc_profile(20.0, 100, 50.0, 100.0, 670.0, z)
nu_2 = visc_profile(20.0, 100, 50.0, 25.0, 670.0, z)
nu_3 = visc_profile(20.0, 100, 50.0, 10.0, 670.0, z)
nu_4 = visc_profile(20.0, 100, 50.0, 100.0, 400.0, z)
nu_5 = visc_profile(20.0, 100, 50.0, 25.0, 400.0, z)
nu_6 = visc_profile(20.0, 100, 50.0, 10.0, 400.0, z)
nu_7 = visc_profile(20.0, 100, 50.0, 5.0, 400.0, z)
plt.plot(nu_1, z, nu_2, z, nu_3, z, nu_4, z, nu_5, z, nu_6, z, nu_7, z)
plt.gca().invert_yaxis()
In [76]:
z = np.linspace(0.0, 3480.0, 1000)
z = np.linspace(0.0, 1000.0, 100)
nu_4 = visc_profile(10.0, 100, 10.0, 100.0, 400.0, z)
nu_6 = visc_profile(10.0, 100, 10.0, 10.0, 400.0, z)
plt.plot(nu_4, z, nu_6, z)
plt.gca().invert_yaxis()
In [72]:
def write_viscosity_file(filename, header, z, nu):
f = open(filename, 'w')
f.write(header)
for zi, nui in zip(z, nu):
f.write("{:10.3f} {:10.3f}\n".format(zi, nui))
f.close()
In [77]:
write_viscosity_file("square_viscosity_profile.dat", "Profile with 10 100 10 100 400\n", z, nu_4)
In [78]:
write_viscosity_file("smooth_viscosity_profile.dat", "Profile with 10 100 10 10 400\n", z, nu_6)
In [ ]: