Code

import matplotlib.pyplot as plt
import numpy as np
import itertools

crustal_thickness = 100E3 # m
depth = np.linspace(1,crustal_thickness,35) # m

### Ref. stress state ----------------------------------------------------------
density = 2700 # kg/m3
gravitational_accel = 9.8 # ms-2
lith_pressure = depth * gravitational_accel * density # Pa
lith_pressure = lith_pressure*1E-6 # MPa

### Geotherm -------------------------------------------------------------------
k = 2 # J s-1 m-1 K-1
S0 = 3E-6 # Wm-3
h_r = 10E3 # m
q_m = 0.025 # Wm-2
C_2 = 300 # K

term_1 = h_r**2 * S0 * np.exp(- (depth/h_r))
term_2 = depth * h_r * S0 * np.exp(- (crustal_thickness/h_r))
term_3 = q_m * depth

T = (-term_1 - term_2 + term_3 + C_2)/k

### Frictional behaviour -------------------------------------------------------
byerlee_shear_stress = []
for i, value in enumerate(lith_pressure):
    if 3 < lith_pressure[i] < 200:
        shear_stress = 0.85 * lith_pressure[i]
        byerlee_shear_stress.append(shear_stress)
    elif lith_pressure[i] > 200:
        shear_stress = 0.6 * lith_pressure[i] + 50
        byerlee_shear_stress.append(shear_stress)
    else:
        byerlee_shear_stress.append(0)

### Viscous behaviour ----------------------------------------------------------
'''
The parameters (A, n and Q) and the calculation of the stress supported by the
viscous flow are defined further down in the plotting loop.
'''

eps_tectonic_list = [1E-14] # s-1
# eps_tectonic_list = [1E-10,1E-12,1E-14] # s-1
material_list = ['Qtz', 'Qtzite', 'Diorite']
R = 8.314

### Plotting -------------------------------------------------------------------
fig, ax = plt.subplots(1,2, figsize = (8,5))

ax[0].plot(T,depth, 'r', label = 'T (C)')
ax[0].plot(lith_pressure,depth, 'k', label = 'P$_{lith}$ (MPa)')

ax[0].legend()
ax[0].invert_yaxis()
ax[0].set_ylim(crustal_thickness, 0)
ax[0].set_xlim(0,1000)
ax[0].set_ylabel('depth (km)')

### Plotting frictional behaviour
ax[1].plot(byerlee_shear_stress, depth, label = 'Frictional strength')

### Plotting viscous behaviour
linestyle_list = itertools.chain([':','--', '-'])
for eps in eps_tectonic_list:
    eps_line = next(linestyle_list)

    # Flow law parameters from practical table
    A_list = itertools.chain([6.3E-5, 6.3E-8, 3.2E-2]) # MPa-n s-1
    n_list = itertools.chain([5.3, 3.1, 2.4])
    Q_list = itertools.chain([213E3, 135E3, 212E3]) # kJ/mol
    colour_list = itertools.chain(['C2','C3', 'C4'])

    for i in material_list:
        if i == 'Qtz':
            print(i)
            A = next(A_list)
            n = next(n_list)
            Q = next(Q_list)
            colour = next(colour_list)

            inverse_n = 1/n
            Ar = - Q/(R*T)
            B = eps/(A*np.exp(Ar))
            flow_stress = np.power(B,inverse_n)

            ax[1].plot(flow_stress, depth, color = colour,linestyle = eps_line, label = '%s, %s' %(i, eps))
        else:
            continue

ax[1].legend(frameon = False)
ax[1].invert_yaxis()
ax[1].set_ylim(crustal_thickness, 0)
ax[1].set_xlim(0,1000)
ax[1].set_xlabel('strength (MPa)')

plt.tight_layout()

### Saving figure
# plt.savefig('crustal_strength_profile.png', dpi = 350)

plt.show()

"""
Strain model used in Gilgannon et al. (2023)

Script documents an example of the transformation of linear objects in a strain
feild. Equations used come from section 4-5 in Ramsay (1967).

AUTHOR
    James Gilgannon
    email: jamesgilgannon@hotmail.com

REVISION
    1 (2023-02-22)

REFERENCES
    [1] Ramsay, J.G. (1967). Folding and fracturing of rock
"""

import numpy as np
import matplotlib.pyplot as plt
import mplstereonet

# Number of measurements
number_of_points = 18
number_of_plunges = 5

### Generate set lines with uniform distriubtion -------------------------------
trends_1_deg = np.linspace(-180, 180,number_of_points)
plunges_1_deg = np.linspace(0, 90, number_of_plunges)
# Apply transformation to make plunges uniform
plunges_1_deg_unif = np.cos(plunges_1_deg)
plunges_1_deg = np.rad2deg(plunges_1_deg_unif)

# Recast data across a grid
xx, yy = np.meshgrid(trends_1_deg, plunges_1_deg)
trends_1_deg = [item for sublist in xx for item in sublist]
plunges_1_deg = [item for sublist in yy for item in sublist]

# convert to radians
trends_1_rad = np.deg2rad(trends_1_deg)
plunges_1_rad =np.deg2rad(plunges_1_deg)

### Calculate directional cosines ----------------------------------------------
# Convestion for axis association: l = x, m = y, n = z
l = np.sin(trends_1_rad)*np.cos(plunges_1_rad)
m = np.cos(trends_1_rad)*np.cos(plunges_1_rad)
n = -np.sin(plunges_1_rad)

### Defines strain feild -------------------------------------------------------
e_x = 0
e_y = 0
e_z = -0.8

### Calculate principle quadratic elongations ----------------------------------
lambda_1 = (1+e_x)**2
lambda_2 = (1+e_y)**2
lambda_3 = (1+e_z)**2

### Preform deformation and rotation of lines ----------------------------------
sum_directional_cosines_quad_strain = (lambda_1*(l**2)) + (lambda_2*(m**2)) + (lambda_3*(n**2))
l_prime_2 = (lambda_1*(l**2))/sum_directional_cosines_quad_strain
m_prime_2 = (lambda_2*(m**2))/sum_directional_cosines_quad_strain
n_prime_2 = (lambda_3*(n**2))/sum_directional_cosines_quad_strain

# New deformed directional cosines
l_prime = np.sqrt(l_prime_2)
m_prime = np.sqrt(m_prime_2)
n_prime = np.sqrt(n_prime_2)

# Ensure correct directionaility because deformation operation makes all cosines +ve
# and collapses the lines into one quadrant
for i, l_value in enumerate(l):
    if l[i] < 0:
        l_prime[i] = -l_prime[i]

for i, m_value in enumerate(m):
    if m[i] < 0:
        m_prime[i] = -m_prime[i]

for i, n_value in enumerate(n):
    if n[i] < 0:
        n_prime[i] = -n_prime[i]

# Convert back to plunge and trend but enusre that trend azi are respected by adding 180
# where appropriate ---------------------------------------------------------------------------
trends_new = []
for i, val in enumerate(m_prime):
    if m_prime[i] <0:
        trend = np.arctan(l_prime[i]/m_prime[i]) + np.pi
        trends_new.append(trend)
    else:
        trend = np.arctan(l_prime[i]/m_prime[i])
        trends_new.append(trend)

plunges_new = -np.arcsin(n_prime)

# Convert to degrees
plunges_new_deg = np.rad2deg(plunges_new)
trends_new_deg = np.rad2deg(trends_new)
# Convert to circular angles appropriate for azimuth
trends_new_deg = [i % 360 for i in trends_new_deg]

### Plotting -------------------------------------------------------------------

# Pole figure params
sigma=5
vmax =5
vmin =0
levels = 5
cmap = 'Greys'


fig = plt.figure()
ax = fig.add_subplot(221, projection='stereonet')
ax.line(plunges_1_deg,trends_1_deg)
ax.set_azimuth_ticks([0,270])
ax.set_azimuth_ticklabels(['X','Y'])

ax = fig.add_subplot(222, projection='stereonet')
ax.line(plunges_new_deg,trends_new_deg)
ax.set_azimuth_ticks([0,270])
ax.set_azimuth_ticklabels(['',''])


ax = fig.add_subplot(223, projection='stereonet')
u = ax.density_contourf(plunges_1_deg,trends_1_deg,
    measurement='lines', cmap =cmap,
    method='exponential_kamb',
    sigma=sigma,
    vmax =vmax,
    levels = levels
    )

ax.set_azimuth_ticks([0,270])
ax.set_azimuth_ticklabels(['',''])
plt.colorbar(u,ax=ax)

ax = fig.add_subplot(224, projection='stereonet')
u = ax.density_contourf(plunges_new_deg,trends_new_deg,
    measurement='lines', cmap =cmap,
    method='exponential_kamb',
    sigma=sigma,
    vmax =vmax,
    levels = levels
    )
ax.set_azimuth_ticks([0,270])
ax.set_azimuth_ticklabels(['',''])

plt.colorbar(u,ax=ax)

plt.show()