# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Isostatic deflection in 2D

# Source: Hodgetts et al. (1998). Flexural modelling of continental lithosphere deformation: a comparison of 2D and 3D techniques, Tectonophysics, 294, 1-2, p.1-2

# These are the equations being solved:
#
# $\left(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}\right) D \left( \frac{\partial^{2} w_{(x,y)}}{\partial y^{2}} + \frac{2 \partial^{2} w_{(x,y)}}{\partial x \partial y} + \frac{\partial^{2} w_{(x,y)}}{\partial x^{2}}\right) + \left( \rho_{m} - \rho \right) g w_{(x,y)} = l_{(x,y)}$
#
# This is solved using a Fourier transform solution:
#
# $W_{(u,v)} = R_{(u,v)} \cdot L_{(u,v)}$
#
# Where $W(u,v)$ is the Fourier transform of the deflections, $L(u,v)$ is the Fourier transform of the surface loads (equal to $\rho g h$), and $R(u,v)$ is a response function, defined as:
#
# $R_{(u,v)} = \frac{1}{\left( \rho_{m} - \rho \right) g + D\left(u^{2} + v^{2}\right)^{2}}$
#
# In the particular case of Curtis' Santa Cruz Mountains problem, we are interested in knowing the rock uplift that is associated with a given amount of crustal thickening.  Noting that $h = t - w$, where $t$ is the thickening, we can rewrite Equation (10) as:
#
# $\left(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}\right) D \left( \frac{\partial^{2} w_{(x,y)}}{\partial y^{2}} + \frac{2 \partial^{2} w_{(x,y)}}{\partial x \partial y} + \frac{\partial^{2} w_{(x,y)}}{\partial x^{2}}\right) + \rho_{m} g w_{(x,y)} = \rho g t_{(x,y)}$
#
# And Equation 12 as:
#
# $\frac{1}{\rho_{m} g + D\left(u^{2} + v^{2}\right)^{2}}$
#
# In this case, $l(x,y)$ now becomes the crustal thickening ($t(x,y)$).
#
# Once deflections are computed, rock uplift can be computed as $u(x,y) = t(x,y) - w(x,y)$

import numpy as np
# %matplotlib widget
import matplotlib.pylab as plt

# We will use a model with E=50 GPa, h = 20 km as a demonstration.  Values of E=10 GPa, h = 5 km produce way too much deflection.  This is probably due to the fact that this model uses thickening rate (and so $(\rho_{m} - \rho)$ in the original form is smaller than $\rho_{m}$), the fact that the model is not periodic, and that the model does not extend forever in the out-of-plane direction.  Note however, that the fraction of Airy isostacy approaches about what we would like ($f\approx 0.38$) for the wavelengths of the SCM.

# +
# Define constants:

E = 10E9
g = 9.8
rho_m = 3200
rho = 2700
h = 5E3
v = 0.25

D = E*np.power(h,3) / (12*(1-np.power(v,2)))

print('D = ', D)

# +
# Define extent of plots:

bounding_box = np.array([[5.149823487603397807e+05, 4.162743473135999404e+06],
                         [5.592764889708703849e+05, 4.195161883133378811e+06],
                         [6.377426705260890303e+05, 4.087951441662845202e+06],
                         [5.934485303155583097e+05, 4.055533031665465795e+06]])

extent = [np.min(bounding_box[:,0]), np.max(bounding_box[:,0]), np.min(bounding_box[:,1]), np.max(bounding_box[:,1])]

# +
# Read thickening grid and define dimensions:

import pickle as p
(X, Y, UX, UY, UZ) = p.load(open('data/dislocation_safonly_nolock.p','rb'))

dx = np.mean(np.diff(X)[1,:])*1000
thickening_disloc = UZ*1E6*4*1000 # This will give us units of meters for a 4Myr model

disloc_extent = np.array([np.min(X[:]), np.max(X[:]), np.min(Y[:]), np.max(Y[:])])*1E3

# +
# Calculate wavenumbers:

wx_disloc = np.fft.fftfreq(thickening_disloc.shape[1],d=dx)*2.0*np.pi
wy_disloc = np.fft.fftfreq(thickening_disloc.shape[0],d=dx)*2.0*np.pi

[WX_disloc, WY_disloc] = np.meshgrid(wx_disloc,wy_disloc)

# +
# Build response function:

R_disloc = np.power(rho_m*g + D*np.power(np.power(WX_disloc,2)+np.power(WY_disloc,2),2),-1)

# +
# Transform thickening grid:

T_disloc = np.fft.fft2(thickening_disloc*rho*g)

# +
# Convolve and back-transform:

W_disloc = R_disloc*T_disloc
w_disloc = np.real(np.fft.ifft2(W_disloc))

# +
# Calculate rock uplift and plot:

u_disloc = thickening_disloc - w_disloc

plt.figure()
plt.title('Dislocation - Thickening (m)')
plt.imshow(thickening_disloc, extent=disloc_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

plt.figure()
plt.title('Dislocation - Deflection (m)')
plt.imshow(w_disloc, extent=disloc_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

plt.figure()
plt.title('Dislocation - Rock / Surface Uplift (m)')
plt.imshow(u_disloc, extent=disloc_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

# +
# Read irregular points for EP model and create regular grid:

xyz_ep = np.loadtxt('data/EP_UTM_surface_nodes.txt')

xy_ep = xyz_ep[:,0:2]
z_ep = xyz_ep[:,2] - 20000.0

ep_extent = [np.min(xy_ep[:,0]), np.max(xy_ep[:,0]), np.min(xy_ep[:,1]), np.max(xy_ep[:,1])]

[Xi, Yi] = np.meshgrid(np.arange(ep_extent[0],ep_extent[1],dx), np.arange(ep_extent[2],ep_extent[3],dx))

from scipy.interpolate import griddata
thickening_ep = griddata(xy_ep, z_ep, (Xi, Yi), method='cubic', fill_value=0.0)


# +
# Transform thickening grid and calculate deflections and rock uplift:

T_ep = np.fft.fft2(thickening_ep*rho*g)

wx_ep = np.fft.fftfreq(thickening_ep.shape[1],d=dx)*2.0*np.pi
wy_ep = np.fft.fftfreq(thickening_ep.shape[0],d=dx)*2.0*np.pi
[WX_ep, WY_ep] = np.meshgrid(wx_ep,wy_ep)

R_ep = np.power(rho_m*g + D*np.power(np.power(WX_ep,2)+np.power(WY_ep,2),2),-1)

W_ep = R_ep*T_ep
w_ep = np.real(np.fft.ifft2(W_ep))
u_ep = thickening_ep - w_ep

# +
# Calculate rock uplift and plot:

# %matplotlib widget
import matplotlib.pylab as plt

u_ep = thickening_ep - w_ep

plt.figure()
plt.title('Elastoplastic - Thickening (m)')
plt.imshow(thickening_ep, extent=ep_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

plt.figure()
plt.title('Elastoplastic - Deflection (m)')
plt.imshow(w_ep, extent=ep_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

plt.figure()
plt.title('Elastoplastic - Rock / Surface Uplift (m)')
plt.imshow(u_ep, extent=ep_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

# +
# Read irregular points for E model and create regular grid:

xyz_e = np.loadtxt('data/E_UTM_surface_nodes.txt')

xy_e = xyz_e[:,0:2]
z_e = xyz_e[:,2] - 20000.0

e_extent = [np.min(xy_e[:,0]), np.max(xy_e[:,0]), np.min(xy_e[:,1]), np.max(xy_e[:,1])]

[Xi, Yi] = np.meshgrid(np.arange(e_extent[0],e_extent[1],dx), np.arange(e_extent[2],e_extent[3],dx))

from scipy.interpolate import griddata
thickening_e = griddata(xy_e, z_e, (Xi, Yi), method='cubic', fill_value=0.0)


# +
# Transform thickening grid and calculate deflections and rock uplift:

T_e = np.fft.fft2(thickening_e*rho*g)

wx_e = np.fft.fftfreq(thickening_e.shape[1],d=dx)*2.0*np.pi
wy_e = np.fft.fftfreq(thickening_e.shape[0],d=dx)*2.0*np.pi
[WX_e, WY_e] = np.meshgrid(wx_e,wy_e)

R_e = np.power(rho_m*g + D*np.power(np.power(WX_e,2)+np.power(WY_e,2),2),-1)

W_e = R_e*T_e
w_e = np.real(np.fft.ifft2(W_e))
u_e = thickening_e - w_e

# +
# Calculate rock uplift and plot:

# %matplotlib widget
import matplotlib.pylab as plt

u_e = thickening_e - w_e

plt.figure()
plt.title('Elastic - Thickening (m)')
plt.imshow(thickening_e, extent=e_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

plt.figure()
plt.title('Elastic - Deflection (m)')
plt.imshow(w_e, extent=e_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

plt.figure()
plt.title('Elastic - Rock / Surface Uplift (m)')
plt.imshow(u_e, extent=e_extent, origin='lower', vmin = 0, vmax = 2500)
plt.axis(extent)
plt.colorbar()
plt.show()

# + jupyter={"source_hidden": true}

# -



StarcoderdataJupyterScriptsDedupFiltered