Help with grdflexure

Hello everyone,

I’m ultimately trying to run a gmt script that outputs a flexure profile for a line load on an elastic plate over a viscous half space, but first I’m testing the simpler problem of a line load on an elastic plate and I’ve pasted the code that I’m using to do this below. The code runs and I’m able to get the plots of the flexure, however the issues that I’m having are when setting the densities. I’ve been trying to compare the elastic solutions from grdflexure to an elastic solution from Anthony Watts, and when the infill and load densities are set to be some non-zero value the two solutions agree very well. However, if I try to set the infill density to be 0 in the GMT code with a nonzero load density, then there is no flexure where the Watts solution predicts flexure. Setting the load density to be 0 in the GMT code with a nonzero infill density does produce flexure, but this flexure is overpredicted compared to the Watt’s solution. I was wondering if there was something that I was missing that might help me understand this command better.

I appreciate the help!
Daniel

#!/usr/bin/env bash

# Testing firmoviscous flexure code

gmt set GMT_FFT kiss

# Set densities

rhol=0

rhom=3300

rhos=2700

rhow=0

# Create a truncated seamount load or 5km height

gmt grdseamount -R-5000000/5000000/-5000000/5000000 -I5000 -Gsmt.nc -Cd -Dk -E -Z0 << EOF

0 0 0 100000 100 10000 75 73

EOF
​
# Select times for calculation

cat << EOF > times.txt

3k red

10k orange

30k green

70k blue

1M black

EOF

rm -f flist

while read t color; do 

gmt grdflexure smt.nc -D${rhom}/${rhol}/${rhos}/${rhow} -E30k -Cp0.5 -Cy3e10 -Ge_%.12g.nc -Na -T$t -L >> flist

done < times.txt

# Get cross-sections of all solutions

gmt math -T-1000000/1000000/1000 0 = t.txt

gmt grdtrack t.txt -G+lflist > a.txt

Hi Daniel. I have never tried the two cases you are discussing. In general, the FFT solution is exact when rlol == rhos. For any other combination you have different domains that the FFT solution no longer applies. In grdflexure I handle this case with an approximation: Use rhos for calculating the flexural wavelength as well as the load, but scale load magnitude to compensate for the lower density. Without this approximation one must solve the flexure problem numerically instead (e.g., finite difference). You can read about the approximation and the justification in Wessel [2016]. See the definition of ksi in my equation (6). It you plug rhos = rhow = 0 in there then you get a zero scale and hence zero flexure. Your cases are not ones I considered when making these approximations. Ref: Wessel, P. (2016), Regional–residual separation of bathymetry and revised estimates of Hawaii plume flux, Geophysical Journal International , 204 (2), 932–947, doi:10.1093/gji/ggv472.