LITHO1.0¶
Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 6 March 2018
LITHO1.0 model: http://igppweb.ucsd.edu/~gabi/litho1.0.html
Pasyanos, Masters, Laske, Ma (2014): "LITHO1.0: An updated crust and lithospheric model of the Earth." J. Geophys. Res., 119(3): 2153–2173, doi:10.1002/2013JB010626
LITHO1.0 describes the lithosphere as a 9-layer sandwich (see image below).
At each (latitude,longitude) point, the density is constant within each layer.
Elevations refer to a reference Earth shape, positive elevations above, negative below.
Density is set to zero where layer thickness is zero.
These two images clarify the file names and grid point ordering:
In the following, I plot those files as maps.
In [1]:
#%reset -f
from mpl_toolkits.basemap import Basemap, cm
import matplotlib.pyplot as plt
import numpy as np
from math import pi
In [2]:
# read-in CRUST1.0 data
c1_iscont = np.loadtxt('crust1_1deg_iscont.txt')
In [3]:
# read-in LITHO1.0 data
l1_top1 = np.loadtxt('litho1_1deg_top1.txt')
l1_top2 = np.loadtxt('litho1_1deg_top2.txt')
l1_top3 = np.loadtxt('litho1_1deg_top3.txt')
l1_top4 = np.loadtxt('litho1_1deg_top4.txt')
l1_top5 = np.loadtxt('litho1_1deg_top5.txt')
l1_top6 = np.loadtxt('litho1_1deg_top6.txt')
l1_top7 = np.loadtxt('litho1_1deg_top7.txt')
l1_top8 = np.loadtxt('litho1_1deg_top8.txt')
l1_top9 = np.loadtxt('litho1_1deg_top9.txt')
l1_top10 = np.loadtxt('litho1_1deg_top10.txt')
l1_rho1 = np.loadtxt('litho1_1deg_rho1.txt')
l1_rho2 = np.loadtxt('litho1_1deg_rho2.txt')
l1_rho3 = np.loadtxt('litho1_1deg_rho3.txt')
l1_rho4 = np.loadtxt('litho1_1deg_rho4.txt')
l1_rho5 = np.loadtxt('litho1_1deg_rho5.txt')
l1_rho6 = np.loadtxt('litho1_1deg_rho6.txt')
l1_rho7 = np.loadtxt('litho1_1deg_rho7.txt')
l1_rho8 = np.loadtxt('litho1_1deg_rho8.txt')
l1_rho9 = np.loadtxt('litho1_1deg_rho9.txt')
l1_vp1 = np.loadtxt('litho1_1deg_vp1.txt') * 1e-3 # in km/s
l1_vp2 = np.loadtxt('litho1_1deg_vp2.txt') * 1e-3
l1_vp3 = np.loadtxt('litho1_1deg_vp3.txt') * 1e-3
l1_vp4 = np.loadtxt('litho1_1deg_vp4.txt') * 1e-3
l1_vp5 = np.loadtxt('litho1_1deg_vp5.txt') * 1e-3
l1_vp6 = np.loadtxt('litho1_1deg_vp6.txt') * 1e-3
l1_vp7 = np.loadtxt('litho1_1deg_vp7.txt') * 1e-3
l1_vp8 = np.loadtxt('litho1_1deg_vp8.txt') * 1e-3
l1_vp9 = np.loadtxt('litho1_1deg_vp9.txt') * 1e-3
l1_vs1 = np.loadtxt('litho1_1deg_vs1.txt') * 1e-3
l1_vs2 = np.loadtxt('litho1_1deg_vs2.txt') * 1e-3
l1_vs3 = np.loadtxt('litho1_1deg_vs3.txt') * 1e-3
l1_vs4 = np.loadtxt('litho1_1deg_vs4.txt') * 1e-3
l1_vs5 = np.loadtxt('litho1_1deg_vs5.txt') * 1e-3
l1_vs6 = np.loadtxt('litho1_1deg_vs6.txt') * 1e-3
l1_vs7 = np.loadtxt('litho1_1deg_vs7.txt') * 1e-3
l1_vs8 = np.loadtxt('litho1_1deg_vs8.txt') * 1e-3
l1_vs9 = np.loadtxt('litho1_1deg_vs9.txt') * 1e-3
l1_d1 = l1_top1 - l1_top2
l1_d2 = l1_top2 - l1_top3
l1_d3 = l1_top3 - l1_top4
l1_d4 = l1_top4 - l1_top5
l1_d5 = l1_top5 - l1_top6
l1_d6 = l1_top6 - l1_top7
l1_d7 = l1_top7 - l1_top8
l1_d8 = l1_top8 - l1_top9
l1_d9 = l1_top9 - l1_top10
l1_d678 = l1_top6 - l1_top9
l1_d6789 = l1_top6 - l1_top10
In [4]:
def fig1(x, y, z, title, unit):
fig, ax = plt.subplots()
m = Basemap(projection='cyl', llcrnrlat=-90, urcrnrlat=90, llcrnrlon=-180, urcrnrlon=180, resolution='c')
m.drawcoastlines()
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
x, y = m(x, y)
clevs = np.linspace(z.min(), z.max(), 19)
cs = m.contourf(x, y, z, clevs, cmap='rainbow')
#cs = m.pcolormesh(x, y, z, vmin=z.min(), vmax=z.max(), cmap='rainbow')
cbar = m.colorbar(cs)
cbar.set_label(unit)
plt.title(title)
return fig, ax
x05 = np.arange(-179.75, 180, 0.5)
y05 = np.arange(89.75, -90, -0.5)
x05, y05 = np.meshgrid(x05, y05)
x1 = np.arange(-179.5, 180, 1.0)
y1 = np.arange(89.5, -90, -1.0)
x1, y1 = np.meshgrid(x1, y1)
Continental mask from CRUST1.0¶
In [5]:
plt.rcParams['figure.figsize'] = (18,9)
In [6]:
z = c1_iscont
fig, ax = fig1(x1, y1, z, 'Continent', '1/0')
plt.show()
Plots of LITHO1.0 elevation (km) of tops of layers¶
In [7]:
z = l1_top1
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top1 (water)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [8]:
z = l1_top2
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top2 (ice)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [9]:
z = l1_top3
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top3 (sediment1)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [10]:
z = l1_top4
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top4 (sediment2)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [11]:
z = l1_top5
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top5 (sediment3)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [12]:
z = l1_top6
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top6 (upper crust)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [13]:
z = l1_top7
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top7 (middle crust)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [14]:
z = l1_top8
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top8 (lower crust)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [15]:
z = l1_top9
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top9 (lithospheric mantle)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [16]:
z = l1_top10
fig, ax = fig1(x1, y1, z, 'LITHO1.0 top10 (asthenospheric mantle)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
Plots of LITHO1.0 layer thicknesses (km)¶
In [17]:
z = l1_d1
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Water thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [18]:
z = l1_d2
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Ice thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [19]:
z = l1_d3
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Sediment1 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [20]:
z = l1_d4
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Sediment2 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [21]:
z = l1_d5
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Sediment3 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [22]:
z = l1_d6
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Upper crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [23]:
z = l1_d7
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Middle crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [24]:
z = l1_d8
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Lower crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [25]:
z = l1_d9
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Lithospheric mantle thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
LITHO1.0 crystalline crust thicknesses, layers 6–8 (km)¶
In [26]:
z = l1_d678
fig, ax = fig1(x1, y1, z, 'LITHO1.0 crustal thickness (layers 6–8)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
Plots of LITHO1.0 density (kg/m3) in layers¶
In [27]:
z = l1_rho1
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho1 (water)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [28]:
z = l1_rho2
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho2 (ice)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [29]:
z = l1_rho3
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho3 (sediment1)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [30]:
z = l1_rho4
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho4 (sediment2)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [31]:
z = l1_rho5
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho5 (sediment3)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [32]:
z = l1_rho6
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho6 (upper crust)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [33]:
z = l1_rho7
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho7 (middle crust)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [34]:
z = l1_rho8
fig, ax = fig1(x1, y1, z, 'LITHO1.0 rho8 (lower crust)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
Layer 9 (lithospheric mantle) has uniform density of 3300 km/m3 at all lat,lon.
In [35]:
z = l1_rho9
print('min =', z.min())
print('max =', z.max())
Plots of LITHO1.0 Vp (km/s)¶
In [36]:
z = l1_vp1
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp1 (water)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [37]:
z = l1_vp2
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp2 (ice)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [38]:
z = l1_vp3
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp3 (sediment1)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [39]:
z = l1_vp4
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp4 (sediment2)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [40]:
z = l1_vp5
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp5 (sediment3)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [41]:
z = l1_vp6
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp6 (upper crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [42]:
z = l1_vp7
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp7 (middle crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [43]:
z = l1_vp8
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp8 (lower crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [44]:
z = l1_vp9
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vp9 (lithospheric mantle)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
Plots of LITHO1.0 Vs (km/s)¶
In [45]:
z = l1_vs1
print('min =', z.min())
print('max =', z.max())
In [46]:
z = l1_vs2
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs2 (ice)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [47]:
z = l1_vs3
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs3 (sediment1)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [48]:
z = l1_vs4
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs4 (sediment2)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [49]:
z = l1_vs5
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs5 (sediment3)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [50]:
z = l1_vs6
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs6 (upper crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [51]:
z = l1_vs7
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs7 (middle crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [52]:
z = l1_vs8
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs8 (lower crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [53]:
z = l1_vs9
fig, ax = fig1(x1, y1, z, 'LITHO1.0 Vs9 (lithospheric mantle)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
Histograms Vp of continental crust¶
In [54]:
plt.rcParams['figure.figsize'] = (18,12)
In [55]:
# weights to account for latitude effect on surface area
nlat = l1_vp6.shape[0]
sintheta = np.sin( np.arange(0.5/nlat, 1, 1/nlat) * pi )
l1_wt = np.zeros(l1_vp6.shape)
for i in range(0,l1_vp6.shape[1]):
l1_wt[:,i] = sintheta[:]
In [56]:
dz = 0.025
zmin = min(l1_vp6.min(), l1_vp7.min(), l1_vp8.min())
zmax = max(l1_vp6.max(), l1_vp7.max(), l1_vp8.max())
zrange = ( (int(zmin/dz)+1)*dz, int(zmax/dz)*dz )
nbins = int((zrange[1]-zrange[0])/dz)
fig = plt.figure()
wt1 = l1_wt.flatten()*c1_iscont.flatten()
#
plt.subplot(3, 1, 1)
plt.hist(l1_vp6.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 Vp6 (upper crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 2)
plt.hist(l1_vp7.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 Vp7 (middle crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 3)
plt.hist(l1_vp8.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 Vp8 (lower crust) – CONTINENTS ONLY')
plt.xlabel('km/s')
plt.ylabel('#')
#
plt.show()
In [57]:
zmin = min(l1_vp6.min(), l1_vp7.min(), l1_vp8.min(), l1_vp9.min())
zmax = max(l1_vp6.max(), l1_vp7.max(), l1_vp8.max(), l1_vp9.max())
zrange = ( (int(zmin/dz)+1)*dz, int(zmax/dz)*dz )
nbins = int((zrange[1]-zrange[0])/dz)
fig = plt.figure()
wt1 = l1_wt.flatten()*c1_iscont.flatten()
#
plt.subplot(4, 1, 1)
plt.hist(l1_vp6.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 Vp6 (upper crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(4, 1, 2)
plt.hist(l1_vp7.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 Vp7 (middle crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(4, 1, 3)
plt.hist(l1_vp8.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 Vp8 (lower crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(4, 1, 4)
plt.hist(l1_vp9.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 Vp9 (lithospheric mantle) – CONTINENTS ONLY')
plt.xlabel('km/s')
plt.ylabel('#')
#
plt.show()
Histograms fractional thickness of continental crustal layers¶
In [58]:
# thickness of crystalline crustal layers
l1_d678 = l1_top6 - l1_top9
l1_fd6 = l1_d6/l1_d678
l1_fd7 = l1_d7/l1_d678
l1_fd8 = l1_d8/l1_d678
dz = 0.01
zmin = min(l1_fd6.min(), l1_fd7.min(), l1_fd8.min())
zmax = max(l1_fd6.max(), l1_fd7.max(), l1_fd8.max())
zrange = ( (int(zmin/dz)+1)*dz, int(zmax/dz)*dz )
nbins = int((zrange[1]-zrange[0])/dz)
fig = plt.figure()
wt1 = l1_wt.flatten()*c1_iscont.flatten()
#
plt.subplot(3, 1, 1)
plt.hist(l1_fd6.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 6/678 thickness (upper crust / crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 2)
plt.hist(l1_fd7.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 7/678 thickness (middle crust / crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 3)
plt.hist(l1_fd8.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('LITHO1.0 8/678 thickness (lower crust / crust) – CONTINENTS ONLY')
plt.xlabel('m/s')
plt.ylabel('#')
#
plt.show()
In [59]:
z = l1_fd6
fig, ax = fig1(x1, y1, z, 'LITHO1.0 lay6/lays678 thickness (upper crust / crust)', 'fraction')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_fd7
fig, ax = fig1(x1, y1, z, 'LITHO1.0 lay7/lays678 thickness (upper crust / crust)', 'fraction')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_fd8
fig, ax = fig1(x1, y1, z, 'LITHO1.0 lay8/lays678 thickness (upper crust / crust)', 'fraction')
plt.show()
print('min =', z.min())
print('max =', z.max())
