Crustal gravity – Jaroslav Jaroš project

Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 5 January 2021

• HARMAN ... subroutines of Martinec 1991 (https://doi.org/10.1016/0010-4655(91)90057-R)
• SHTOOLS ... subroutines of Wieczorek & Meschede 2018 (https://doi.org/10.1029/2018GC007529, https://shtools.github.io/SHTOOLS/)
• VELIMSKY ... subroutines by Jakub Velímský (personal communication Nov/Dec 2020)

import numpy as np
import matplotlib.pyplot as plt

Testing HARMAN vs. SHTOOLS vs. VELIMSKY

Findings:

• HARMAN has an issue at degree 1, as the power at j=1 decreases with increasing jmax. Power at j=1 should probably be a numerical zero – although SHTOOLS gives a consistent non-zero value (but 10 orders of magnitude below degree 0).
• HARMAN, SHTOOLS, VELIMSKY are consistent up to degree ~250.
• HARMAN gives inconsistent results for j>300, where power depends on jmax. (It may be due to the character of the input dataset, Vr_XGM2019.unf, which is a model up to jmax=720, however, based on GOCE model up to jmax=300.)
• SHTOOLS gives consistent but incorrect results for j>250

Will be using VELIMSKY code.

def figure_spectrum():

    dir = 'jaros_project/compare_SHA_GOCE/'
    
    x0, y0 = np.loadtxt(dir+'powj_harman_449.jp', unpack=True)
    x1, y1 = np.loadtxt(dir+'powj_shtools_449.jp', unpack=True)
    x2, y2 = np.loadtxt(dir+'powj_velimsky_449.jp', unpack=True)

    plt.figure(figsize=(16, 14))

    plt.subplot(211)
    plt.semilogy(x0, y0, label='harman')
    plt.semilogy(x1, y1, label='shtools')
    plt.semilogy(x2, y2, label='velimsky')
    plt.xlabel('Spherical harmonic degree j')
    plt.ylabel('Power at degree j')
    plt.legend()
    plt.title('Power spectrum')

    plt.subplot(212)
    n = 9
    plt.semilogy(x0[0:n], y0[0:n], 'o-', label='harman')
    plt.semilogy(x1[0:n], y1[0:n], 'o-', label='shtools')
    plt.semilogy(x2[0:n], y2[0:n], 'o-', label='velimsky')
    plt.xlabel('Spherical harmonic degree j')
    plt.ylabel('Power at degree j')
    plt.legend()

    plt.show()

figure_spectrum()

plt.figure(figsize=(16, 6))
dir = 'jaros_project/compare_SHA_GOCE/'
x, y = np.loadtxt(dir+'powj_harman_008.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 8')
x, y = np.loadtxt(dir+'powj_harman_016.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 16')
x, y = np.loadtxt(dir+'powj_harman_032.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 232')
x, y = np.loadtxt(dir+'powj_harman_064.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 64')
x, y = np.loadtxt(dir+'powj_harman_128.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 128')
x, y = np.loadtxt(dir+'powj_harman_256.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 256')
x, y = np.loadtxt(dir+'powj_harman_449.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 449')
x, y = np.loadtxt(dir+'powj_harman_512.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 512')
x, y = np.loadtxt(dir+'powj_harman_897.jp', unpack=True) ;  plt.semilogy(x, y, label='jmax = 897')
#plt.xlim(0,8)
plt.legend()
plt.xlabel('Spherical harmonic degree j')
plt.ylabel('Power at degree j')
plt.title('HARMAN')
plt.show()

Plots of GOCE data using GMT4

def figure_powj(filename, jcut):
    x, y = np.loadtxt(filename, unpack=True)
    plt.figure(figsize=(16, 5))
    plt.subplot(121)
    plt.semilogy(x, y, label=filename)
    plt.xlabel('degree j')
    plt.legend()
    plt.subplot(122)
    plt.semilogy(x, y, 'o-')
    plt.xlabel('degree j')
    plt.xlim(0,16)
    plt.show()

V

figure_powj('jaros_project/prepareGOCEdata/V_GOCE.jp', 16)

Vr

figure_powj('jaros_project/prepareGOCEdata/Vr_GOCE.jp', 16)

Vrr

figure_powj('jaros_project/prepareGOCEdata/Vrr_GOCE.jp', 16)

T

figure_powj('jaros_project/prepareGOCEdata/T_GOCE.jp', 16)

Tr

figure_powj('jaros_project/prepareGOCEdata/Tr_GOCE.jp', 16)

Trr

figure_powj('jaros_project/prepareGOCEdata/Trr_GOCE.jp', 16)

Convergence of calculation with layered crust

Looking at max and min value of calculated V_rr on a 120 x 90 grid.

With ~30 layers the calculation converges at 4 significant digits.

def figure_convergence():
    n, t, min, max = np.loadtxt('jaros_project/figures/convergence.dat', unpack=True)

    plt.figure(figsize=(8, 11))

    plt.subplot(311)
    plt.plot(n, t, 'o-')
    plt.xlabel('Number of layers n')
    plt.ylabel('CPU time in seconds')
    plt.title('Convergence ov V_rr, computed on grid 120 x 60')

    plt.subplot(312)
    plt.semilogx(n, max, 'o-')
    plt.xlabel('Number of layers n')
    plt.ylabel('max V_rr')
    #plt.legend()

    plt.subplot(313)
    plt.semilogx(n, min, 'o-')
    plt.xlabel('Number of layers n')
    plt.ylabel('min V_rr')

    plt.show()

figure_convergence()

Comparison of CRUST1.0 prediction vs. GOCE measurement

def figure_spectrum(label):
    prefix = 'jaros_project/calculation/'+label
    x0, y0 = np.loadtxt(prefix+'_meas.jp', unpack=True)
    x1, y1 = np.loadtxt(prefix+'_pred.jp', unpack=True)
    x2, y2 = np.loadtxt(prefix+'_diff.jp', unpack=True)

    plt.figure(figsize=(10, 10))

    plt.subplot(211)
    plt.semilogy(x0, y0, label='GOCE')
    plt.semilogy(x1, y1, label='prediction')
    plt.semilogy(x2, y2, label='difference')
    plt.xlabel('Spherical harmonic degree j')
    plt.ylabel('Power at degree j')
    plt.legend()
    plt.title('Power spectrum')

    plt.subplot(212)
    n = 31
    plt.semilogy(x0[0:n], y0[0:n], 'o-', label='GOCE')
    plt.semilogy(x1[0:n], y1[0:n], 'o-', label='prediction')
    plt.semilogy(x2[0:n], y2[0:n], 'o-', label='difference')
    plt.xlabel('Spherical harmonic degree j')
    plt.ylabel('Power at degree j')
    plt.legend()

    plt.show()

def figure_correlation(label):
    filein = 'jaros_project/calculation/'+label+'_corr.jc'
    x0, y0 = np.loadtxt(filein, unpack=True)

    plt.figure(figsize=(10, 10))

    plt.subplot(211)
    plt.plot(x0, y0)
    plt.plot((x0[0],x0[-1]), (0,0), 'k--')
    plt.xlabel('Spherical harmonic degree j')
    plt.ylabel('Correlation at degree j')
    plt.title('Correlation spectrum')

    plt.subplot(212)
    n = 31
    plt.plot(x0[0:n], y0[0:n], 'o-')
    plt.plot((x0[0],x0[n-1]), (0,0), 'k--')
    plt.xlabel('Spherical harmonic degree j')
    plt.ylabel('Correlation at degree j')

    plt.show()

Trr

• RMS pred = 1135 mE -> 1156 mE
• RMS meas = 246.0 mE -> 254.1 mE
• RMS diff = 1133 mE -> 1179 mE
• Correlation pred/meas = 0.01883

figure_spectrum('Trr')

figure_correlation('Trr')

Tr

• RMS pred = 74.26 mGal -> 76.70 mGal
• RMS meas = 16.73 mGal -> 17.13 mGal
• RMS diff = 72.88 mGal -> 78.19 mGal
• Correlation pred/meas = 0.02383

figure_spectrum('Tr')

figure_correlation('Tr')

T

• RMS pred = 1146 m2/s2 -> 1224 m2/s2
• RMS meas = 240.4 m2/s2 -> 253.6 m2/s2
• RMS diff = 1099 m2/s2 -> 1204 m2/s2
• Correlation pred/meas = 0.1835

figure_spectrum('T')

figure_correlation('T')