Thin shell displacements, moho-relief, and strains on Mars

This notebook contains a series of examples that detail how to call several functions that solve the Banerdt (1986) thin shell model under different assumptions. Various improvements have been made to the model including the possibility to account for finite-amplitude correction and filtering (Wieczorek & Phillips, 1998), lateral density variations at any arbitrary depth and within the surface, flexure, or moho reliefs (Wieczorek et al., 2013), and density difference between the surface topography and crust (Broquet & Wieczorek, 2019).

The model links and outputs 8 parameters expressed in spherical harmonics (degree $l$, order $m$): the topography ($H_{lm}$), geoid at the surface ($G_{lm}$), geoid at the moho depth ($Gc_{lm}$), net acting load on the lithosphere ($q_{lm}$), tangential load potential ($\Omega_{lm}$), flexure of the lithosphere ($w_{lm}$), crustal root variations ($\delta c_{lm}$), and internal density variations ($\delta \rho_{lm}$).

These 8 parameters are linked through an initial set of 5 equations for $q_{lm}$, $w_{lm}$, $G_{lm}$, $Gc_{lm}$, and $\Omega_{lm}$ (see Banerdt, 1986). Thus, if 3 conditions are given (e.g., $H_{lm}$, $G_{lm}$, and the assumption that $\delta \rho_{lm}$ = 0), the system will be evenly determined, and the stresses and strains along with the corresponding moho relief can be calculated.

To get started, import the standard matplotlib library for graphics, cmcrameri to use perceptually uniform colormaps, numpy for mathematical extensions to Python, pyshtools (version >= 4.9) for spherical harmonic manipulations and plotting, and the Displacement_strain_planet package.

import os
import matplotlib.pyplot as plt
from cmcrameri import cm
import numpy as np
import pyshtools as pysh
from Displacement_strain_planet import *

For simplicity, we here we make use of the pyshtools package to load the gravity field (spherical harmonic model of Genova et al., 2016) and topography (spherical harmonic model of Wieczorek et al., 2015) of Mars and to perform some spherical harmonics manipulations.

For Mars, the $C_{20}$ coefficient of the topography and gravity are strongly influenced by the planet's rotational flattening. These are generally entirely excluded (or ~95%) for loading calculation as this relief should not be considered as a load.

lmax = 90 # Maximum spherical harmonic degree to perform all calculations
pot_clm   = pysh.datasets.Mars.GMM3(lmax=lmax)  
topo_clm  = pysh.datasets.Mars.MarsTopo2600(lmax=lmax)


R       = topo_clm.coeffs[0,0,0]     # Mean planetary radius
pot_clm = pot_clm.change_ref(r0 = R) # Downward continue to Mean planetary radius

# Compute the geoid as approximated in Banerdt's formulation
geoid_clm = pot_clm * R

# Constants
G      = pysh.constants.G.value        # Gravitational constant
gm     = pot_clm.gm                    # GM given in the gravity model file
mass   = gm / G                        # Mass of the planet
g0     = gm / R**2                     # Mean gravitational attraction of the planet


# Remove 100% of C20
percent_C20 = 0. 
topo_clm.coeffs[0,2,0]  = (percent_C20 /100.) * topo_clm.coeffs[0,2,0]
geoid_clm.coeffs[0,2,0] = (percent_C20 /100.) * geoid_clm.coeffs[0,2,0]

We set some initial parameters for the upcoming inversion and send it to arguments.

# Parameters
c    = 50e3    # Mean Crustal thickness
Te   = 100e3   # Elastic thickness
rhom = 3500.   # Mantle density
rhoc = 2800.   # Crustal density
rhol = 2800.   # Surface density
E    = 100e9   # Young's modulus
v    = 0.25    # Poisson's ratio

args_param_m   = (g0, R, c, Te, rhom, rhoc, rhol, lmax, E, v, mass)
args_expand = dict(lmax = 5*lmax, lmax_calc = lmax)
args_fig = dict(figsize=(12,10), dpi= 100)
args_plot = dict(tick_interval=[45, 30], colorbar = 'bottom', cmap = cm.roma_r)
args_ticksr = dict(ticks = 'wSnE', ylabel = None)

We define some initial arrays for the examples below. zeros is a spherical harmonic array of zeros. drhom_lm_plume is an example array with spherical harmonic coefficients for a mantle plume with a density contrast of -30 kg m$^{-3}$, which tapers at about 1000 km radius, located underneath Tharsis.

path           = '%s/data' %(os.getcwd())
zeros          = pysh.SHCoeffs.from_zeros(lmax = lmax).coeffs
drhom_lm_plume = pysh.SHCoeffs.from_file('%s/Example_plume_coeffs_clm.txt' %(path))
drhom_lm_plume = drhom_lm_plume.pad(lmax = lmax).coeffs

Crustal thickness and displacements using global present-day gravity and topography

Compensation by crustal root variations and displacements

In this example, we solve for the displacement of the surface of Mars by calling the function Thin_shell_matrix_nmax, assuming that the gravity and topography of the planet are compensated by a combination of crustal root variations and flexure. 3 assumptions are required to solve the system, and we here assume that the observed topography and geoid are known (H_lm = topo_clm.coeffs, G_lm = geoid_clm.coeffs), and that there are no density variations in the interior (drhom_lm = zeros.copy()).

By default, the code accounts for finite-amplitude corrections and solve the model in an iterative manner. More informations, including crustal thickness values at each iterations, can be printed by setting the optional parameter quiet to False. We will also make use of a downward continuation minimum-amplitude filter to damp unrealistic oscillations of the moho-relief (Wieczorek & Phillips, 1998). For this, we call the optional argument filter, set it to 'Ma', and set the degree at which the filter equals 0.5 to 50 with a call to filter_half. The user can also input any given filter using the argument filter_in.

The function outputs the following spherical harmonic coefficients: \ w_lm flexure (defined positive upward), \ A_lm poloidal term of the tangential displacement, \ moho_relief_lm moho relief, \ dc_lm crustal root variations (defined negative upward), \ drhom_lm internal density variations, \ omega_lm tangential load potential, \ q_lm net load on the lithosphere, \ Gc_lm geoid at the moho depth, \ G_lm geoid at the surface, and \ H_lm topography.

And the linear solution sols are expressed as lambda functions (size lmax+1) of all components. Lambda functions can be used to re-calculate the same problem with different inputs very fast.

(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm,
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m, G_lm = geoid_clm.coeffs,
        H_lm = topo_clm.coeffs, drhom_lm = zeros.copy(), filter = 'Ma', filter_half = 50)

We next plot $w_{lm}$, $\delta c_{lm}$, $\delta \rho_{lm}$, and the moho-relief ($w_{lm} - \delta c_{lm}$). To obtain crustal thickness from the moho-relief, we add the topography coefficients (H_lm).

For the plots, we make use of the pyshtools software to convert the spherical harmonic arrays to a SHCoeffs class, expand it to a grid, and plot the results. For the upper left figure plot, we added a zero contour to delimit where the displacement is negative and positive.

fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2, **args_fig)

grid_W = pysh.SHCoeffs.from_array(w_lm/1e3).expand(**args_expand) - R/1e3
grid_W.plot(ax=ax1, cb_label = 'Upward displacement (km)', **args_plot)
ax1.contour(grid_W.data > 0, levels = [0.99], extent = (0,360,-90,90), 
        colors = 'k', origin = 'upper') # Add zero displacement contour
pysh.SHCoeffs.from_array(dc_lm/1e3).expand(**args_expand).plot(ax=ax2, 
        cb_label = 'Crustal root variations (km)', **args_ticksr, 
        **args_plot, cmap_limits = [-25,25])
(pysh.SHCoeffs.from_array((H_lm-moho_relief_lm)/1e3).expand(**args_expand)).\
        plot(ax=ax3, cb_label = 'Crustal thickness (km)', 
        show = True, **args_plot)
ax4.set_visible(False)

Comparison with the mass-sheet case

By default, the Thin_shell_matrix_nmax function computes the non-linear finite-amplitude corrections in an iterative manner. Let's compare what the difference would be assuming a mass-sheet case, which can be obtained by setting the optional nmax parameter to 1. Setting nmax to 1 significantly fasten all computations.

(w_lm_n1, A_lm_n1, moho_relief_lm_n1, dc_lm_n1, drhom_lm_n1, omega_lm_n1, 
q_lm_n1, Gc_lm_n1, G_lm_n1, H_lm_n1, sols) = Thin_shell_matrix_nmax(*args_param_m, 
        G_lm = geoid_clm.coeffs, H_lm = topo_clm.coeffs, drhom_lm = zeros.copy(), nmax = 1, 
        filter = 'Ma', filter_half = 50)

### Plotting
fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2, **args_fig)
W_diff = pysh.SHCoeffs.from_array((w_lm-w_lm_n1)/1e3).expand(**args_expand)
W_diff.plot(ax=ax1, cb_label = '$\\Delta$ Upward displacement (km)', **args_plot)
pysh.SHCoeffs.from_array((dc_lm-dc_lm_n1)/1e3).expand(**args_expand).plot(ax=ax2, 
        cb_label = '$\\Delta$ crustal root variations (km)', **args_ticksr, **args_plot)
(pysh.SHCoeffs.from_array(((moho_relief_lm-moho_relief_lm_n1))/1e3).expand(**args_expand)).\
        plot(ax=ax3, cb_label = '$\\Delta$ Moho relief (km)', show = True, **args_plot)
ax4.set_visible(False)

Pratt compensation (crustal density variations)

In this example, we will assume that the observed topography and geoid of Mars are compensated by density variations in the crust only (i.e., Pratt isostasy). In that case, we need to input 2 parameters: $H_{lm}$, $G_{lm}$, and set both $\delta c_{lm}$ and $w_{lm}$ to zero. The density contrast is located between the moho depth and surface, and we thus set the 2 optional parameters top_drho and base_drho to 0 (surface depth) and c (moho depth).

With these 4 inputs, the model becomes overdetermined (4 unknowns, 5 equations). One option is to remove an equation from the system. For consistency, we will here remove one of the five equations that is not that of an input parameter. The simplest way to do that is to use the argument remove_equation of the function and set it to, for example, 'Gc_lm'.

Notes

• It is possible that the system gives different outputs as a function of the equation that is removed (e.g., if we remove G_lm).
• Pratt isostasy could also be viewed as only requiring that there is no crustal root variations (i.e., only setting $\boldsymbol{\delta c_{lm}}$ to zero).
• The script below assumes isostasy, but uses observed gravity. Pratt isostasy can be derived without knowledge of the gravity field as a simple pressure balance between columns of crust. See below for the Airy case.

A simple sanity check can be performed for the case of Pratt compensation. This requires setting the elastic thickness (Te) to a large value in order to have almost zero displacement, in which case compensation would only occur by density variations. For large elastic thicknesses, it is thus not necessary to set $w_{lm}$ to zero and to remove an equation.

As demonstrated below, both approaches give roughly the same results. Slight differences are due to a long-wavelength residual displacement (due to the elastic thickness not being infinite), finite-amplitude corrections, and variations in density within the surface relief making the equations non-linear (short-wavelength differences).

top_drho = 0.
base_drho = c

(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm, 
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m, top_drho = top_drho, 
        base_drho = base_drho, G_lm = geoid_clm.coeffs, H_lm = topo_clm.coeffs, 
        dc_lm = zeros.copy(), w_lm = zeros.copy(), remove_equation = 'Gc_lm',
        filter = 'Ma', filter_half = 50)

# Alternative Te = 2000 km comparing drhom_lm
drhom_lm_Te = Thin_shell_matrix_nmax(g0, R, c, 2000e3, rhom, rhoc, rhol, 
        lmax, E, v, mass, top_drho = top_drho, base_drho = base_drho, 
        G_lm = geoid_clm.coeffs, H_lm = topo_clm.coeffs, dc_lm = zeros.copy(),                      
        filter = 'Ma', filter_half = 50)[4]

### Plotting
fig, (ax1,ax2) = plt.subplots(1,2, **args_fig)
# Adding in the mean crustal density
(pysh.SHCoeffs.from_array(drhom_lm).expand(**args_expand)+rhoc).plot(ax=ax1, 
        cb_label = 'Lateral density variations in the crust (kg m$^{-3}$)', **args_plot)
pysh.SHCoeffs.from_array(drhom_lm-drhom_lm_Te).expand(**args_expand).plot(ax=ax2, 
        cb_label = '$\\Delta$ Lateral density variations (kg m$^{-3}$)', show = True,
        **args_plot, **args_ticksr)

Airy compensation (crustal thickness variations) when knowing the planetary gravity field

In this example, we will assume that the observed topography and geoid of Mars are compensated by crustal thickness variations (i.e., Airy isostasy). Airy isostasy can be achieved when there is no internal density variations and when the elastic thickness is zero.

Notes

• In order to avoid numerical problems, when the elastic thickness is zero, the function internally sets it to 1 m.
• Airy isostasy could also be viewed as only requiring that there is no internal density variations (i.e., only setting $\boldsymbol{\delta \rho_{lm}}$ to zero and not assuming an elastic thickness).
• The script below assumes isostasy, but uses observed gravity. Airy isostasy can be derived without knowledge of the gravity field as a simple pressure balance between columns of crust, see below.

args_param_m_Te0   = (g0, R, c, 0, rhom, rhoc, rhol, lmax, E, v, mass)

(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm, 
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m_Te0, 
        G_lm = geoid_clm.coeffs, H_lm = topo_clm.coeffs, 
        drhom_lm = zeros.copy(), filter = 'Ma', filter_half = 50)

# Difference without filter
moho_relief_lm_no_filter = Thin_shell_matrix_nmax(*args_param_m_Te0, 
        G_lm = geoid_clm.coeffs, H_lm = topo_clm.coeffs, drhom_lm = zeros.copy())[2]

### Plotting
fig, (ax1,ax2) = plt.subplots(1,2, **args_fig)
pysh.SHCoeffs.from_array((H_lm-moho_relief_lm)/1e3).expand(**args_expand).plot(ax=ax1, 
        cb_label = 'Crustal thickness variations (km)', **args_plot)
(pysh.SHCoeffs.from_array((moho_relief_lm-moho_relief_lm_no_filter)/1e3).expand(**args_expand)).\
        plot(ax=ax2, cb_label = ('$\\Delta$ Crustal thickness variations with' +
        ' vs without filter (km)'), show = True,**args_plot, **args_ticksr)

! Warning:Thin_shell_matrix, elastic thickness set to 1 m to avoid numerical problems !
! Warning:Thin_shell_matrix, elastic thickness set to 1 m to avoid numerical problems !

Airy compensation without making use of the planetary gravity field

Airy isostasy is defined as having equal pressures in equal columns of crust. This scenario can be achieved by using the optional add_equation argument. add_equation takes a symbolic string equation as input and can be used to link variables in the system of equations that is solved.

In our case, Airy isostasy translates to 'H_lm * g0 * rhoc = dc_lm * gmoho * (rhom-rhoc)' or add_equation = 'H_lm * g0 * rhoc - dc_lm * gmoho * (rhom-rhoc)' (the system automatically assumes it is equal to zero).

As this equation provides a new constraint to the system, we only need to specify 2 constraints in addition to add_equation.

#Define some constants to compute gmoho
rhobar = mass * 3.0 / 4.0 / np.pi / R ** 3
gmoho = (g0 * (1.0 + (((R - c) / R) ** 3 - 1.0) * rhoc / rhobar) / ((R - c) / R) ** 2)

add_equation = 'H_lm * %s - dc_lm * %s' %(g0*rhoc, gmoho*(rhom-rhoc))

(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm,
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m,
        H_lm = topo_clm.coeffs, drhom_lm = zeros.copy(), add_equation = add_equation)

### Plotting
fig, (ax1,ax2) = plt.subplots(1,2, **args_fig)

pysh.SHCoeffs.from_array(dc_lm/1e3).expand(**args_expand).plot(ax=ax1, 
        cb_label = 'Crustal root variations variations (km)', **args_ticksr, 
        **args_plot)
(pysh.SHCoeffs.from_array((H_lm-moho_relief_lm)/1e3).expand(**args_expand)).\
        plot(ax=ax2, cb_label = 'Crustal thickness (km)', 
        show = True, **args_plot)

Displacements only due to a plume underneath Tharsis

In this example, we will solve for the displacement assuming that a plume is located underneath an elastic shell of thickness Te. All parameters but the density variations associated with the plume (drhom_lm) are unknown. In that case, the system is largely underdetermined and we need a few assumptions.

1) The first assumption is that there will be no crustal root variations, and we set $\delta c_{lm}$ to zero.

2) A second assumption is that the plume will give rise to dynamic topography only. In that case, we want $H_{lm} = w_{lm}$ or $H_{lm} - w_{lm} = 0$. The function allows the user to add an equation to the system, and the required assumption is obtained by setting the optional argument add_equation to be equal to 'H_lm - w_lm' (the system automatically assumes it is equal to zero). We note that any linear equation that links one (or more) of the 8 parameters can be input.

plume_thickness = 150e3 # Thickness of the mantle plume
top_drho = Te
base_drho = Te + plume_thickness

(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm, 
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m, top_drho = 
        top_drho, base_drho = base_drho, drhom_lm = drhom_lm_plume, 
        dc_lm = zeros.copy(), add_equation = 'H_lm - w_lm')

### Plotting
fig, (ax1,ax2) = plt.subplots(1,2, **args_fig)
(pysh.SHCoeffs.from_array(w_lm/1e3).expand(**args_expand) - R/1e3).plot(ax=ax1, 
        cb_label = 'Upward displacement (km)', **args_plot)
(pysh.SHCoeffs.from_array(drhom_lm).expand(**args_expand)+rhom).plot(ax=ax2, 
        cb_label = 'Input density variations in the mantle (kg m$^{-3}$)', show = True, 
        **args_plot, **args_ticksr)

More complex add_equation option

Degree-dependent compensation

Let's make a funny example where we force the model to compensate the degree-1 loads by crustal density variations and the rest by crustal thickness variations to show an additional possibility in the function. In that case, the additional equation becomes '(w_lm - dc_lm) * add_array1 + drhom_lm * add_array2', where 'add_array1' and 'add_array2' contain the necessary spherical harmonic coefficients (see below), which are passed through the argument add_arrays. More arrays can be input as 'add_array3', 'add_array4' and so on. All must be with shape equal to the input constraints (2,lmax+1,lmax+1).

Notes

• For multiplications in add_equation, and in order to avoid non-linearities we only use the (0,lmax+1,0) terms of each add_arrays. All terms (2,lmax+1,lmax+1) are used with additions or substractions.

• Multiplications in add_equation can generate problems when used with finite-amplitude correction (nmax != 1). For example, this will introduce correction terms beyond degree-1 for drhom_lm, which is not consistent with the input equation.

top_drho = 0.
base_drho = c
add_equation = '(w_lm - dc_lm) * add_array1 + drhom_lm * add_array2'

add_arrays = np.zeros((2, 2,lmax+1,lmax+1))
add_arrays[0,0,1,0] = 1 # add_array1 is 1 at degree 1 implying w_lm - dc_lm = 0
add_arrays[1,:,:,:] = 1 # add_array is 0 everywhere (implying drhom_lm = 0) but at degree 1.
add_arrays[1,0,1,0] = 0

(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm,
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m, G_lm = geoid_clm.coeffs, 
        H_lm = topo_clm.coeffs, base_drho = base_drho, top_drho = top_drho, 
        add_equation = add_equation, 
        add_arrays = add_arrays, filter = 'Ma', filter_half = 50)

### Plotting
fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2, **args_fig)
(pysh.SHCoeffs.from_array(w_lm/1e3).expand(**args_expand) - R/1e3).plot(ax=ax1,
        cb_label = 'Upward displacement (km)', **args_plot)
pysh.SHCoeffs.from_array(dc_lm/1e3).expand(**args_expand).plot(ax=ax2,
        cb_label = 'Crustal root variations (km)',
        **args_ticksr, **args_plot)
(pysh.SHCoeffs.from_array(drhom_lm).expand(**args_expand)+rhoc).plot(ax=ax3,
        cb_label = 'Lateral density variations in the crust (kg m$^{-3}$)', **args_plot)
(pysh.SHCoeffs.from_array((H_lm-moho_relief_lm)/1e3).expand(**args_expand)).\
        plot(ax=ax4, cb_label = 'Crustal thickness (km)', show = True, 
        **args_plot, **args_ticksr)

Internal load in phase with surface topography

In this last example, we will assume that part of the bottom loading is due to an internal density anomaly located between the surface and 100 km depth. As the surface load depends on the unknown flexure, the easiest is to use the argument add_equation and set it to 'drhom_lm * M - L * rhol * (H_lm - G_lm - w_lm)' where L is the load ratio set to 0.10, and M is the internal load thickness set to 100 km.

load_thickness = 100e3
L_ratio = 0.1
add_equation = 'drhom_lm - (%s * %s * (H_lm - G_lm - w_lm)) / %s' %(L_ratio, rhol, load_thickness)

(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm, 
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m, base_drho = Te+load_thickness, 
        top_drho = Te, G_lm = geoid_clm.coeffs, H_lm = topo_clm.coeffs, 
        add_equation = add_equation, filter = 'Ma', filter_half = 50)

### Plotting
fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2, **args_fig)
(pysh.SHCoeffs.from_array(w_lm/1e3) - R/1e3).expand(**args_expand).plot(ax=ax1, 
        cb_label = 'Upward displacement (km)', **args_plot)
(pysh.SHCoeffs.from_array(drhom_lm)+rhoc).expand(**args_expand).plot(ax=ax2, 
        cb_label = 'Lateral density variations in the crust (kg m$^{-3}$)',  
        **args_ticksr , **args_plot)
(pysh.SHCoeffs.from_array((H_lm-moho_relief_lm)/1e3).expand(**args_expand)).\
        plot(ax=ax3, cb_label = 'Crustal thickness (km)', show = True, **args_plot)
ax4.set_visible(False)

Stresses and strains

Here, we make use of the predicted displacement w_lm and the poloidal term of the tangential displacement A_lm, from the first example, to predict stresses and strains on the lithosphere using the function call Displacement_strains.

The Displacement_strains_shtools function outputs stresses as a function of colatitude ($\theta$) and longitude ($\phi$) as: $\sigma_{\theta}$ = stress_theta, $\sigma_{\phi}$ = stress_phi, and $\sigma_{\theta\phi}$ = stress_theta_phi. Elongation and sheering deformations of the middle surface of the shell are given by: $\varepsilon_{\theta}$ = eps_theta, $\varepsilon_{\phi}$ = eps_phi, and $\varepsilon_{\theta\phi}$ = omega. Bending and twisting deformation of the shell are given by: $\kappa_{\theta}$ = kappa_theta, $\kappa_{\phi}$ = kappa_phi, $\kappa_{\theta\phi}$ = tau. The total deformation is given by tot_theta, tot_phi, and tot_thetaphi.

In these scripts, compression is defined to be negative.

The strains are then rotated into their principal components using the function call Principal_strainstress_angle to output the minimum, maximum, and sum of the principal horizontal strains (min_strain, max_strain, sum_strain). The function also outputs the principal strain angle principal_angle.

# Recompute parameters assuming compensation due to crustal root
# variations and displacement (first example)
(w_lm, A_lm, moho_relief_lm, dc_lm, drhom_lm, omega_lm, q_lm, Gc_lm, 
G_lm, H_lm, sols) = Thin_shell_matrix_nmax(*args_param_m, G_lm = geoid_clm.coeffs, 
        H_lm = topo_clm.coeffs, drhom_lm = zeros.copy(), filter = 'Ma', filter_half = 50)

# Strain
args_param_s = (E, v, R, Te, lmax)
(stress_theta, stress_phi, stress_theta_phi, eps_theta, eps_phi, 
omega, kappa_theta, kappa_phi, tau, tot_theta, tot_phi, tot_thetaphi) = Displacement_strains_shtools(A_lm, w_lm, 
        *args_param_s)

# Principal strains
min_strain, max_strain, sum_strain, principal_angle = \
        Principal_strainstress_angle(-tot_theta, -tot_phi, -tot_thetaphi)

Plotting

In ax3 and ax4, we plot the extensional tectonic features mapped by Knampeyer et al. (2006). They are seen to nicely cluster where the model predicts extension (plotted negative), and follow the rotational pattern given by the principal strain angle.

args_plot = dict(tick_interval=[45, 30], colorbar = 'bottom', cmap = cm.vik)
fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2, **args_fig)

pysh.SHGrid.from_array(min_strain * 1e3).plot(ax=ax1, 
        cb_label = 'Minimum principal horizontal strain ($\\times 10^{-3}$)', 
        cmap_limits = [-4, 4], **args_plot)
pysh.SHGrid.from_array(max_strain * 1e3).plot(ax=ax2,
        cb_label = 'Maximum principal horizontal strain ($\\times 10^{-3}$)', 
        cmap_limits = [-4, 4], **args_plot, **args_ticksr)
pysh.SHGrid.from_array(sum_strain * 1e3).plot(ax=ax3, 
        cb_label = 'Sum of principal horizontal strains ($\\times 10^{-3}$)', 
        cmap_limits = [-3, 3], **args_plot)
pysh.SHGrid.from_array(principal_angle).plot(ax=ax4,
        cb_label = 'Principal angle (°)', 
        cmap_limits = [-90, 90], **args_ticksr, tick_interval=[45, 30], 
        colorbar = 'bottom', cmap = cm.vikO)

# Plot strain direction
skip_i = int(lmax / 5)
skip =(slice(None,None,skip_i),slice(None,None,skip_i))
grid_long, grid_lat = np.meshgrid(np.linspace(0,360,np.shape(principal_angle)[1]),
        np.linspace(90,-90,np.shape(principal_angle)[0]))
ones = np.ones(np.shape(principal_angle))
ax4.quiver(grid_long[skip], grid_lat[skip], ones[skip], ones[skip], scale = 5e1, 
           angles = principal_angle[skip], color = 'g')

# Add the Knapmeyer et al. (2006) mapping of extensional tectonic features
Plt_tecto_Mars(path, ax = [ax3,ax4], extension = True, compression = False)
plt.show()