LITHO1.0 for gravity assessment

Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 31 Jan 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.

I interpolated the LITHO1.0 model onto a 0.5 degree by 0.5 degree grid in longitude and latitude.

I provide ASCII files which give the elevations of interfaces in km and densities in layers in kg/m3 at these 720(lon) by 360(lat) gridpoints:

  • litho1_0.5deg_top1.txt
  • litho1_0.5deg_top2.txt
  • litho1_0.5deg_top3.txt
  • litho1_0.5deg_top4.txt
  • litho1_0.5deg_top5.txt
  • litho1_0.5deg_top6.txt
  • litho1_0.5deg_top7.txt
  • litho1_0.5deg_top8.txt
  • litho1_0.5deg_top9.txt
  • litho1_0.5deg_top10.txt
  • litho1_0.5deg_rho1.txt
  • litho1_0.5deg_rho2.txt
  • litho1_0.5deg_rho3.txt
  • litho1_0.5deg_rho4.txt
  • litho1_0.5deg_rho5.txt
  • litho1_0.5deg_rho6.txt
  • litho1_0.5deg_rho7.txt
  • litho1_0.5deg_rho8.txt
  • litho1_0.5deg_rho9.txt

Elevations refer to a reference Earth shape, positive elevations above, negative below.

Density is set to zero where layer thickness is zero.

These 19 files can be downloaded in a single zip file at: http://geo.mff.cuni.cz/~sramek/ipynb/LITHO1_for_gravity_assessment.zip

These two images clarify the file names and grid point ordering:

title

title

In the following, I just 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
In [2]:
top1 = np.loadtxt('litho1_0.5deg_top1.txt')
top2 = np.loadtxt('litho1_0.5deg_top2.txt')
top3 = np.loadtxt('litho1_0.5deg_top3.txt')
top4 = np.loadtxt('litho1_0.5deg_top4.txt')
top5 = np.loadtxt('litho1_0.5deg_top5.txt')
top6 = np.loadtxt('litho1_0.5deg_top6.txt')
top7 = np.loadtxt('litho1_0.5deg_top7.txt')
top8 = np.loadtxt('litho1_0.5deg_top8.txt')
top9 = np.loadtxt('litho1_0.5deg_top9.txt')
top10 = np.loadtxt('litho1_0.5deg_top10.txt')

rho1 = np.loadtxt('litho1_0.5deg_rho1.txt')
rho2 = np.loadtxt('litho1_0.5deg_rho2.txt')
rho3 = np.loadtxt('litho1_0.5deg_rho3.txt')
rho4 = np.loadtxt('litho1_0.5deg_rho4.txt')
rho5 = np.loadtxt('litho1_0.5deg_rho5.txt')
rho6 = np.loadtxt('litho1_0.5deg_rho6.txt')
rho7 = np.loadtxt('litho1_0.5deg_rho7.txt')
rho8 = np.loadtxt('litho1_0.5deg_rho8.txt')
rho9 = np.loadtxt('litho1_0.5deg_rho9.txt')
In [3]:
def fig_top(x, y, z, title):
    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('km')
    plt.title(title)
    return fig, ax

def fig_rho(x, y, z, title):
    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('kg/m3')
    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)
In [4]:
plt.rcParams['figure.figsize'] = (18,9)

Plots of elevation (km) of tops of layers

In [5]:
fig, ax = fig_top(x05, y05, top1, 'top1 (ice)')
plt.show()
In [6]:
fig, ax = fig_top(x05, y05, top2, 'top2 (water)')
plt.show()
In [7]:
fig, ax = fig_top(x05, y05, top3, 'top3 (sediment1)')
plt.show()
In [8]:
fig, ax = fig_top(x05, y05, top4, 'top4 (sediment2)')
plt.show()
In [9]:
fig, ax = fig_top(x05, y05, top5, 'top5 (sediment3)')
plt.show()
In [10]:
fig, ax = fig_top(x05, y05, top6, 'top6 (upper crust)')
plt.show()
In [11]:
fig, ax = fig_top(x05, y05, top7, 'top7 (middle crust)')
plt.show()
In [12]:
fig, ax = fig_top(x05, y05, top8, 'top8 (lower crust)')
plt.show()
In [13]:
fig, ax = fig_top(x05, y05, top9, 'top9 (lithospheric mantle)')
plt.show()
In [14]:
fig, ax = fig_top(x05, y05, top10, 'top10 (asthenospheric mantle)')
plt.show()

Plots of density (kg/m3) in layers

In [15]:
fig, ax = fig_rho(x05, y05, rho1, 'rho1 (ice)')
plt.show()
In [16]:
fig, ax = fig_rho(x05, y05, rho2, 'rho2 (water)')
plt.show()
In [17]:
fig, ax = fig_rho(x05, y05, rho3, 'rho3 (sediment1)')
plt.show()
In [18]:
fig, ax = fig_rho(x05, y05, rho4, 'rho4 (sediment2)')
plt.show()
In [19]:
fig, ax = fig_rho(x05, y05, rho5, 'rho5 (sediment3)')
plt.show()
In [20]:
fig, ax = fig_rho(x05, y05, rho6, 'rho6 (upper crust)')
plt.show()
In [21]:
fig, ax = fig_rho(x05, y05, rho7, 'rho7 (middle crust)')
plt.show()