|
| 1 | +""" |
| 2 | +Model CPO profile of the GRIP ice core, Greenland |
| 3 | +""" |
| 4 | + |
| 5 | +import numpy as np |
| 6 | +from scipy.integrate import solve_ivp |
| 7 | +from scipy import interpolate |
| 8 | +import pandas as pd |
| 9 | + |
| 10 | +from specfabpy import specfab as sf |
| 11 | +from specfabpy import common as sfcom |
| 12 | +from specfabpy import plotting as sfplt |
| 13 | + |
| 14 | +import matplotlib.pyplot as plt |
| 15 | + |
| 16 | +from matplotlib import rc |
| 17 | +#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) |
| 18 | +rc('font',**{'family':'serif','serif':['Palatino']}) |
| 19 | +rc('text', usetex=True) |
| 20 | + |
| 21 | +""" |
| 22 | +Setup |
| 23 | +""" |
| 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 | +tau = H/a # e-folding time scale |
| 30 | +ugrad = -1/tau*np.diag([-0.5, -0.5, 1]) # uniaxial compression along z-axis |
| 31 | + |
| 32 | +### Numerics |
| 33 | + |
| 34 | +tend = 50e3 # time (in years) to trace out trajectory, starting from the surface |
| 35 | +#tend = -tau*np.log(0.05) # alternatively, set tend so that simulation stops at 95% thinning |
| 36 | + |
| 37 | +Nt = 1000 # number of time steps taken (increase this until results are robust) |
| 38 | +ti = np.linspace(0, tend, Nt) # time points at which to evaluate the solution |
| 39 | +z = np.exp(-ti/tau) # relative height above bed at each point in time |
| 40 | + |
| 41 | +L = 12 # CPO expansion series truncation |
| 42 | +kw_ivp = dict(method='RK45', vectorized=False) # kwargs for solve_ivp() |
| 43 | + |
| 44 | +### CPO dynamics |
| 45 | + |
| 46 | +# Lattice rotation |
| 47 | +iota, zeta = 1, 0 # "deck of cards" behavior |
| 48 | + |
| 49 | +# DDRX |
| 50 | +A = 1.1e7 # rate prefactor (tunable parameter) |
| 51 | +Q = 3.36e4 # activation energy (see Richards et al. (2021) and Lilien et al. (2023)) |
| 52 | +R = 8.314 # gas constant |
| 53 | +Gamma0 = lambda D, T: A*np.sqrt(np.einsum('ij,ji',D,D)/2)*np.exp(-Q/(R*(T+273.15))) # DDRX rate factor |
| 54 | + |
| 55 | +### Temperature profile |
| 56 | + |
| 57 | +df = pd.read_csv('../../../data/icecores/GRIP/temperature.csv') # fetch from github |
| 58 | +fz = interpolate.interp1d(df['zrel'].to_numpy(), df['T'].to_numpy(), kind='linear', fill_value='extrapolate') |
| 59 | +T = fz(z) # temperature profile |
| 60 | + |
| 61 | +""" |
| 62 | +Solve CPO evolution |
| 63 | +""" |
| 64 | + |
| 65 | +D = (ugrad+np.transpose(ugrad))/2 # symmetric part (strain rate tensor) |
| 66 | +W = (ugrad-np.transpose(ugrad))/2 # anti-symmetric part (spin tensor) |
| 67 | +S = D # stress tensor (assume coaxiality with strain-rate tensor; magnitude does not matter for our purpose) |
| 68 | + |
| 69 | +#T[:] = -60 # no DDRX if ice is very cold |
| 70 | + |
| 71 | +lm, nlm_len = sf.init(L) # initialize specfab |
| 72 | +lh_init = 0.25 # initial horizontal eigenvalues |
| 73 | +a2_init = np.diag([lh_init, lh_init, 1-2*lh_init]) # initial a2 (at surface) |
| 74 | + |
| 75 | +nlm_init = np.zeros((nlm_len), dtype=np.complex64) # CPO expansion coefficients |
| 76 | +nlm_init[:sf.L2len] = sf.a2_to_nlm(a2_init) # initial state vector (at surface) |
| 77 | + |
| 78 | +def ODE(t, nlm): |
| 79 | + # d/dt nlm_i = M_ij . nlm_j, where nlm_i is the state vector (aka s_i) |
| 80 | + I = np.argmin(np.abs(ti-t)) # index for closest point on trajectory at time t |
| 81 | + M = sf.M_LROT(nlm, D, W, iota, zeta) # lattice rotation always present |
| 82 | + M += Gamma0(D,T[I])*sf.M_DDRX(nlm, S) # DDRX active if sufficient warm |
| 83 | + #M += Lambda0*sf.M_CDRX(nlm) # CDRX (neglected in this example) |
| 84 | + M += sf.M_REG(nlm, D) # regularization |
| 85 | + return np.matmul(M, nlm) |
| 86 | + |
| 87 | +nlm = solve_ivp(ODE, (0, tend), nlm_init, t_eval=ti, vectorized=False).y.T # CPO state along trajectory |
| 88 | +mi, lami = sfcom.eigenframe(nlm) # a2 eigenvectors and eigenvalues |
| 89 | + |
| 90 | +""" |
| 91 | +Plot results |
| 92 | +""" |
| 93 | + |
| 94 | +### Plot modeled eigenvalues |
| 95 | + |
| 96 | +fig = plt.figure(figsize=(3,4)) |
| 97 | +ax = plt.subplot(111) |
| 98 | + |
| 99 | +c1,c2,c3 = 'tab:blue', 'tab:red', 'k' |
| 100 | + |
| 101 | +ax.plot(lami[:,0], z, '-', c=c1, label=r'$\lambda_1$') |
| 102 | +ax.plot(lami[:,1], z, '-', c=c2, label=r'$\lambda_2$') |
| 103 | +ax.plot(lami[:,2], z, '--', c=c3, label=r'$\lambda_3$') |
| 104 | + |
| 105 | +ax.legend(loc=1, fancybox=False, frameon=False) |
| 106 | +ax.set_title(r'GRIP ice core') |
| 107 | + |
| 108 | +ax.set_xlabel(r'$\lambda_i$') |
| 109 | +ax.set_xticks(np.arange(0,1+.01,0.2)) |
| 110 | +ax.set_xlim([0,1]) |
| 111 | + |
| 112 | +ax.set_ylabel(r'$z/H$') |
| 113 | +ax.set_yticks(np.arange(0,1+.01,0.1)) |
| 114 | +ax.set_ylim([0,1]) |
| 115 | + |
| 116 | +### Plot CPOs |
| 117 | + |
| 118 | +geo, prj = sfplt.getprojection(rotation=45, inclination=50) |
| 119 | + |
| 120 | +def plotCPO(ax, nlm, p0, HW=0.2, cmap='Greys'): |
| 121 | + axtrans = ax.transData.transform(p0) |
| 122 | + trans = fig.transFigure.inverted().transform(axtrans) |
| 123 | + axin = plt.axes([trans[0]-HW/2, trans[1]-HW/2, HW,HW], projection=prj) |
| 124 | + axin.set_global() |
| 125 | + lvlset = [np.linspace(0.05, 0.45, 8), lambda x,p:'%.1f'%x] |
| 126 | + sfplt.plotODF(nlm, lm, axin, lvlset=lvlset, cmap=cmap, showcb=False, nchunk=None) |
| 127 | + sfplt.plotcoordaxes(axin, geo, negaxes=False, color=sfplt.c_dred, axislabels='xi') |
| 128 | + return axin |
| 129 | + |
| 130 | +for zi in np.linspace(0.1, 0.9, 4): |
| 131 | + I = np.argmin(np.abs(z-zi)) |
| 132 | + plotCPO(ax, nlm[I], (1.2,zi)) |
| 133 | + |
| 134 | +### Plot observations |
| 135 | + |
| 136 | +df = pd.read_csv('../../../data/icecores/GRIP/orientations.csv') # fetch from github |
| 137 | +zobs = df['zrel'].to_numpy() |
| 138 | +kw = dict(marker='o', facecolor='none', zorder=1) |
| 139 | +ax.scatter(df['lam1'].to_numpy(), zobs, edgecolor=c1, **kw) |
| 140 | +ax.scatter(df['lam2'].to_numpy(), zobs, edgecolor=c2, **kw) |
| 141 | +ax.scatter(df['lam3'].to_numpy(), zobs, edgecolor=c3, **kw) |
| 142 | + |
| 143 | +### Save plot |
| 144 | + |
| 145 | +plt.savefig('GRIP-parcel.png', dpi=175, pad_inches=0.1, bbox_inches='tight') |
0 commit comments