|
| 1 | +""" |
| 2 | +Modeled CPO evolution for a SSA parcel over Pine Island Glacier, Antarctica |
| 3 | +""" |
| 4 | + |
| 5 | +import numpy as np |
| 6 | +from scipy import interpolate |
| 7 | +import pandas as pd |
| 8 | + |
| 9 | +import matplotlib.pyplot as plt |
| 10 | + |
| 11 | +from matplotlib import rc |
| 12 | +#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) |
| 13 | +rc('font',**{'family':'serif','serif':['Palatino']}) |
| 14 | +rc('text', usetex=True) |
| 15 | + |
| 16 | +from specfabpy import specfab as sf |
| 17 | +from specfabpy import common as sfcom |
| 18 | +from specfabpy import plotting as sfplt |
| 19 | + |
| 20 | +### Init |
| 21 | + |
| 22 | +L = 12 # expansion series truncation |
| 23 | +lm, nlm_len = sf.init(L) |
| 24 | + |
| 25 | +### Velocity gradient experienced by parcel |
| 26 | + |
| 27 | +H = 3027 # ice thickness (Montagnat et al., 2014) |
| 28 | +a = 0.24 # meter ice equiv. per yr (Montagnat et al., 2014) |
| 29 | + |
| 30 | +tau = H/a # e-folding time scale |
| 31 | +ugrad = -1/tau * np.diag([-0.5, -0.5, 1]) # uniaxial compression along z-axis |
| 32 | +D = (ugrad+np.transpose(ugrad))/2 # symmetric part (strain rate tensor) |
| 33 | +W = (ugrad-np.transpose(ugrad))/2 # anti-symmetric part (spin tensor) |
| 34 | +S = D # stress tensor (assume coaxiality with strain-rate tensor; magnitude does not matter for our purpose) |
| 35 | + |
| 36 | +### Fabric dynamics |
| 37 | + |
| 38 | +# Lattice rotation |
| 39 | +iota, zeta = 1, 0 # "deck of cards" behavior |
| 40 | + |
| 41 | +# DDRX |
| 42 | +A = 1.1e7 # rate prefactor (tunable parameter) |
| 43 | +Q = 3.36e4 # activation energy (see Richards et al. (2021) and Lilien et al. (2023)) |
| 44 | +R = 8.314 # gas constant |
| 45 | +Gamma0 = lambda D, T: A*np.sqrt(np.einsum('ij,ji',D,D)/2)*np.exp(-Q/(R*(T+273.15))) # DDRX rate factor |
| 46 | + |
| 47 | +### Numerics |
| 48 | + |
| 49 | +Nt = 500 # number of time steps |
| 50 | +dt = 100 # time step size (yr) |
| 51 | +ti = np.arange(0,Nt) * dt # time vector |
| 52 | +zi = np.exp(-ti/tau) # relative height above bed at each point in time |
| 53 | + |
| 54 | +### Temperature profile |
| 55 | + |
| 56 | +df = pd.read_csv('../../../data/icecores/GRIP/temperature.csv') # fetch from github |
| 57 | +f = interpolate.interp1d(df['zrel'].to_numpy(), df['T'].to_numpy(), kind='nearest', fill_value='extrapolate') |
| 58 | +Ti = f(zi) # temperature vector |
| 59 | +#Ti[:] = -60 # no DDRX if very cold |
| 60 | + |
| 61 | +### Initial fabric state |
| 62 | + |
| 63 | +nlm = np.zeros((Nt,nlm_len), dtype=np.complex64) # state vector |
| 64 | +lami = np.zeros((Nt,3)) # a2 eigenvalues |
| 65 | + |
| 66 | +lxy = 0.25 # initial horizontal eigenvalues |
| 67 | +a2_0 = np.diag([lxy, lxy, 1-2*lxy]) # initial a2 surface state |
| 68 | +nlm[0,:sf.L2len] = sf.a2_to_nlm(a2_0) # initial state vector |
| 69 | +lami[0] = sfcom.eigenframe(nlm[0])[1] # eigenvalues of initial state |
| 70 | + |
| 71 | +### Euler integration |
| 72 | + |
| 73 | +for tt in np.arange(1,Nt): |
| 74 | + nlm_0 = nlm[tt-1,:] # previous solution |
| 75 | + T = Ti[tt] # temperature from borehole measurements |
| 76 | + M_LROT = sf.M_LROT(nlm_0, D, W, iota, zeta) # lattice rotation operator |
| 77 | + M_DDRX = Gamma0(D,T)*sf.M_DDRX(nlm_0, S) # DDRX operator |
| 78 | + M_REG = sf.M_REG(nlm_0, D) # regularization operator |
| 79 | + M = M_LROT + M_DDRX + M_REG |
| 80 | + nlm[tt] = nlm_0 + dt*np.matmul(M, nlm_0) # Euler step |
| 81 | + lami[tt] = sfcom.eigenframe(nlm[tt])[1] |
| 82 | + |
| 83 | +### Plot modeled eigenvalues |
| 84 | + |
| 85 | +fig = plt.figure(figsize=(3,4)) |
| 86 | +ax = plt.subplot(111) |
| 87 | + |
| 88 | +c1,c2,c3 = 'tab:green', 'tab:red', 'k' |
| 89 | + |
| 90 | +ax.plot(lami[:,0], zi, '-', c=c1, label=r'$\lambda_1$') |
| 91 | +ax.plot(lami[:,1], zi, '-', c=c2, label=r'$\lambda_2$') |
| 92 | +ax.plot(lami[:,2], zi, '--', c=c3, label=r'$\lambda_3$') |
| 93 | + |
| 94 | +ax.legend(loc=1, fancybox=False, frameon=False) |
| 95 | +ax.set_title(r'GRIP ice core') |
| 96 | + |
| 97 | +ax.set_xlabel(r'$\lambda_i$') |
| 98 | +ax.set_xticks(np.arange(0,1+.01,0.2)) |
| 99 | +ax.set_xlim([0,1]) |
| 100 | + |
| 101 | +ax.set_ylabel(r'$z/H$') |
| 102 | +ax.set_yticks(np.arange(0,1+.01,0.1)) |
| 103 | +ax.set_ylim([0,1]) |
| 104 | + |
| 105 | +### Plot modeled CPOs |
| 106 | + |
| 107 | +geo, prj = sfplt.getprojection(rotation=45, inclination=50) |
| 108 | + |
| 109 | +def plotCPO(ax, nlm, p0, HW=0.2, cmap='Greys'): |
| 110 | + axtrans = ax.transData.transform(p0) |
| 111 | + trans = fig.transFigure.inverted().transform(axtrans) |
| 112 | + axin = plt.axes([trans[0]-HW/2, trans[1]-HW/2, HW,HW], projection=prj) |
| 113 | + axin.set_global() |
| 114 | + lvlset = [np.linspace(0.05, 0.45, 8), lambda x,p:'%.1f'%x] |
| 115 | + sfplt.plotODF(nlm, lm, axin, lvlset=lvlset, cmap=cmap, showcb=False, nchunk=None) |
| 116 | + sfplt.plotcoordaxes(axin, geo, negaxes=False, color=sfplt.c_dred, axislabels='xi') |
| 117 | + return axin |
| 118 | + |
| 119 | +for _ in np.linspace(0.1, 0.9, 4): |
| 120 | + I = np.argmin(np.abs(zi-_)) |
| 121 | + plotCPO(ax, nlm[I], (1.2,_)) |
| 122 | + |
| 123 | +### Plot observations |
| 124 | + |
| 125 | +df = pd.read_csv('../../../data/icecores/GRIP/orientations.csv') # fetch from github |
| 126 | +zi = df['zrel'].to_numpy() |
| 127 | + |
| 128 | +kw = dict(marker='o', facecolor='none', zorder=1) |
| 129 | +ax.scatter(df['lam1'].to_numpy(), zi, edgecolor=c1, **kw) |
| 130 | +ax.scatter(df['lam2'].to_numpy(), zi, edgecolor=c2, **kw) |
| 131 | +ax.scatter(df['lam3'].to_numpy(), zi, edgecolor=c3, **kw) |
| 132 | + |
| 133 | +### Save plot |
| 134 | + |
| 135 | +plt.savefig('GRIP.png', dpi=175, pad_inches=0.1, bbox_inches='tight') |
0 commit comments