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