Final crustal thickness

Q&A
3

Hi Thomas,

One way to create a file containing present-day crustal thickness that could be used within GPlates would be via some python code (fairly sure this is not possible within GPlates itself).

Suppose you have three arrays (numpy or python lists), containing the longitude, latitude and crustal thickness at some arbitrary points. These could be used to create a gpml ‘scalar coverage’ like this:

multi_point = pygplates.MultiPointOnSphere(zip(latitude_array,longitude_array))

scalar_coverages = {
    pygplates.ScalarType.create_gpml('CrustalThickness'): thickness_array}

ct_feature = pygplates.Feature()
ct_feature.set_geometry((multi_point,scalar_coverages))
ct_feature.set_name('Crustal Thickness')

output_feature_collection = pygplates.FeatureCollection(ct_feature)

If you save this feature collection to a file, you can then load it into GPlates together with a deforming model (presumably one with meshes that overlap to some degree with the point locations), then set the point layer to ‘reconstruct by topologies’. Now the thickness values should change with reconstruction time, but with these changes computed with respect to their present-day values.

Incidentally, to get the crustal thickness values in the first place, a fairly convenient way would be to use the litho1.0 model via this python module:

Then do something like this:

import litho1pt0 as litho
top_depth = litho.layer_depth(latitude_array, longitude_array, 'CRUST1-TOP')
bottom_depth = litho.layer_depth(latitude_array, longitude_array, 'CRUST3-BOTTOM')
thickness_array = bottom_depth-top_depth

hope that helps,
Simon

1 Like