Ctrl+K

EMSC Week 11 - A detailed look at the depth-age relationship for the seafloor.

EMSC Week 11 - A detailed look at the depth-age relationship for the seafloor.#

Extra - what happens if we try to fit the data ?

We have a theoretical relationship that we have been trying to see in the data

\[ d = d_0 + \alpha \sqrt{t \,} \]

We can always try to see what happens if we fit the data to find \(d_0\) and \(\alpha\)

First a grid at fine resolution#

import stripy
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

even_mesh = stripy.spherical_meshes.icosahedral_mesh(include_face_points=True, tree=True, refinement_levels=6)
number_of_mesh_points = even_mesh.npoints

latitudes_in_radians = even_mesh.lats
latitudes_in_degrees = np.degrees(latitudes_in_radians) 

longitudes_in_radians = even_mesh.lons 
longitudes_in_degrees = np.degrees(longitudes_in_radians)%360.0 - 180.0

Find the age and depth values on these points#

As before, we interpolate each of our datasets to the same set of grid points. First we need to define the interpolation routine again.

def map_raster_to_mesh(mesh, latlongrid):
    
    raster = latlongrid.T

    latitudes_in_radians  = mesh.lats
    longitudes_in_radians = mesh.lons 
    latitudes_in_degrees  = np.degrees(latitudes_in_radians) 
    longitudes_in_degrees = np.degrees(longitudes_in_radians)%360.0 - 180.0

    dlons = np.mod(longitudes_in_degrees+180.0, 360.0)
    dlats = np.mod(latitudes_in_degrees+90, 180.0)

    ilons = raster.shape[0] * dlons / 360.0
    ilats = raster.shape[1] * dlats / 180.0

    icoords = np.array((ilons, ilats))

    from scipy import ndimage

    mvals = ndimage.map_coordinates(raster, icoords , order=3, mode='nearest').astype(float)
    
    return mvals

Interpolate age data to fine, triangular grid.#

(You can plot the results as before to see that you have not made a mistake)

plt.figure(figsize=(6, 6))
ax = plt.subplot(111, projection=ccrs.Orthographic(central_longitude=0.1))

ax.add_feature(coastline, edgecolor="black", linewidth=0.5, zorder=3)

plt.scatter(longitudes_in_degrees, latitudes_in_degrees, c=meshages, cmap="RdYlBu",
            vmin=0, vmax=250, s=5,
            transform=ccrs.Geodetic())

import xarray

age_dataset = "data/age.3.6.nc"
age_data = xarray.open_dataset(age_dataset)
subs_data = age_data.sel(x=slice(-180,180, 1), y=slice(-90, 90, 1))

lons = subs_data.coords.get('x')
lats = subs_data.coords.get('y')
vals = subs_data['z']

x,y = np.meshgrid(lons.data, lats.data)
age = vals.data / 100.0
age[np.isnan(age)] = -1.0

meshages = map_raster_to_mesh(even_mesh, age)

(left, bottom, right, top) = (-180, -90, 180, 90)
map_extent = ( left, right, bottom, top)

etopo_dataset = "http://thredds.socib.es/thredds/dodsC/ancillary_data/bathymetry/ETOPO1_Bed_g_gmt4.nc"
etopo_data = xarray.open_dataset(etopo_dataset)
subs_data = etopo_data.sel(x=slice(left,right, 30), y=slice(bottom, top, 30))

lons = subs_data.coords.get('x')
lats = subs_data.coords.get('y')
vals = subs_data['z']

x,y = np.meshgrid(lons.data, lats.data)
height = vals.data

meshheights = map_raster_to_mesh(even_mesh, height)

Interpolate height data to fine, triangular grid.#

(You can plot the results as before to see that you have not made a mistake)

You also should make a decision about the resolution of the data you want to download.

plt.figure(figsize=(6, 6))
ax = plt.subplot(111, projection=ccrs.Orthographic(central_longitude=0.1))

ax.add_feature(coastline, edgecolor="black", linewidth=0.5, zorder=3)

plt.scatter(longitudes_in_degrees, latitudes_in_degrees, c=meshheights, cmap="terrain",
            vmin=-5000, vmax=5000, s=2,
            transform=ccrs.Geodetic())

plt.figure(figsize=(6, 6))
ax = plt.subplot(111)
ax.set_xlim(0,150)
ax.set_ylim(-7000,-2000)

plt.scatter( meshages[meshheights<-2000],
             meshheights[meshheights<-2000], 
             alpha=0.2, marker="+")

Now use scipy.optimise#

The scipy bundle of useful tools has some curve-fitting routines.

You can read the instructions for scipy.optimise.curve_fit() at the scipy website

from scipy import optimize
import numpy as np

## First define the function that we are going to fit

def depth_age_fn(age, d0, alpha):
    return d0 + alpha * np.sqrt(age)

## Now we need to check the data (interpolation can do weird things)

print(meshages.min(), meshages.max())
print(meshheights.min(), meshheights.max())

## Only fit data for valid ages

valid_age_locations = (meshages > 0)

meshages_1 = meshages[valid_age_locations]
meshheights_1 = meshheights[valid_age_locations]

## Now limit the data to the deep ocean

valid_depth_locations = (meshheights_1 < -2000)

meshages_2 = meshages_1[valid_depth_locations]
meshheights_2 = meshheights_1[valid_depth_locations]

## Something like this:
## popt, pcov = optimize.curve_fit(f, ages, heights )
## d0, alpha = popt
## fitted_heights = depth_age_fn(ages, d0, alpha)

Plot the results of the fitting exercise#

Is this what you would expect from looking at this by eye ?

Now add the curves from Stein & Stein (1992)

\[\begin{split} h(t) = \begin{cases} -2600 - 365 t^{1/2} & \textrm{when} \quad t < 20 \textrm{Myr}\\ -5651 + 2473 \exp \left(-0.0278 t \right) & \textrm{when} \quad t \ge 20 \textrm{Myr}\\ \end{cases} \end{split}\]

Note their values for \(d_0\) and \(\alpha\) !

Contents