gFlex¶

Multiple methods to solve elastic plate flexure, designed for applications to Earth’s lithosphere.

gFlex computes lithospheric flexural isostasy — the bending of Earth’s elastic outer shell under surface loads such as ice sheets, sediment, lava flows, or water. Both one-dimensional (profile) and two-dimensional (map view) solutions are supported, using either a finite-difference method (which handles spatially variable elastic thickness) or superposition of analytical solutions (fast, for constant elastic thickness).

The deflection \(w\) satisfies

\[\nabla^2(D\,\nabla^2 w) - T_e\,\boldsymbol{\sigma} : \nabla\nabla w + \Delta\rho\, g\, w = q,\]

where \(\boldsymbol{\sigma} : \nabla\nabla w\) is the double contraction of the in-plane stress tensor with the Hessian of \(w\), \(q\) [Pa] is the applied surface normal stress, \(\Delta\rho = \rho_m - \rho_\text{fill}\) [kg m⁻³] is the mantle minus infill density, \(g\) [m s⁻²] is gravitational acceleration, and \(D = E T_e^3 / \bigl[12(1 - \nu^2)\bigr]\) is the flexural rigidity (\(E\) = Young’s modulus, \(T_e\) = elastic thickness, \(\nu\) = Poisson’s ratio). See Theory for the full expanded equations and physical interpretation of each term.

Note

When you use gFlex, please cite:

Wickert, A. D. (2016), Open-source modular solutions for flexural isostasy: gFlex v1.0, Geosci. Model Dev., 9(3), 997–1017.

Installation¶

pip install gflex

gFlex requires Python ≥ 3.10 and depends on NumPy, SciPy, and Matplotlib.

Quick start¶

2-D finite-difference deflection under a rectangular load:

import numpy as np
from gflex import F2D

flex = F2D()
flex.Quiet = True
flex.Method = 'FD'
flex.PlateSolutionType = 'vWC1994'
flex.Solver = 'direct'
flex.g = 9.8
flex.E = 65e9
flex.nu = 0.25
flex.rho_m = 3300.
flex.rho_fill = 0.
flex.Te = 35e3 * np.ones((50, 50))   # uniform 35 km elastic thickness
flex.qs = np.zeros((50, 50))
flex.qs[10:40, 10:40] = 1e6          # 150 × 150 km load at 1 MPa
flex.dx = flex.dy = 5000.            # 5 km grid
flex.BC_W = flex.BC_S = flex.BC_N = '0Displacement0Slope'
flex.BC_E = '0Moment0Shear'
flex.initialize()
flex.run()
flex.finalize()

deflection = flex.w   # (50, 50) array; negative values = downward

Configuration files¶

As an alternative to the programmatic API, gFlex can be driven by a configuration file — passed to the gflex CLI or to the F1D / F2D constructor. Two formats are supported:

YAML (.yaml / .yml) — recommended for new workflows. See input/input_f1d.yaml and input/input_f2d.yaml for complete 1-D and 2-D examples.
INI (any extension) — the legacy format; see input/input_f1d and input/input_f2d.

A minimal 2-D YAML configuration:

mode:
  dimension: 2
  method: FD
  PlateSolutionType: vWC1994
parameter:
  YoungsModulus: 6.5e10
  PoissonsRatio: 0.25
  GravAccel: 9.8
  MantleDensity: 3300
  InfillMaterialDensity: 0
input:
  Loads: path/to/loads.txt
  ElasticThickness: path/to/Te.txt
output:
  Plot: both
numerical:
  GridSpacing_x: 4000
  BoundaryCondition_West: 0Moment0Shear
  BoundaryCondition_East: 0Displacement0Slope
  Solver: direct
  ConvergenceTolerance: 1.0e-3
numerical2D:
  GridSpacing_y: 4000
  BoundaryCondition_North: Mirror
  BoundaryCondition_South: 0Slope0Shear

Run either format from the command line:

gflex path/to/config.yaml   # YAML (extension required)
gflex path/to/config        # INI (any extension)

See Configuration Files for a full parameter reference and annotated examples.

The package also provides domain-padding utilities for variable-Te grids (see pad_domain() for 2-D and pad_domain_1d() for 1-D), a flexural wavelength calculator (see flexural_wavelengths()), a Landlab component, a CSDMS Basic Model Interface (BmiGflex), and a command-line entry point (gflex <config_file>).

Contents