The model gFlex implements multiple (user-selectable) solution methods to solve for flexure and isostasy due to surface loading in both one (line loads) and two (point loads) dimensions. It works for elastic lithospheric plates of both constant and spatially variable elastic thickness and allows the user to select the solution method.
This page contains some information on gFlex, but its main documentation is in the README.md file, displayed at the gFlex GitHub repository at https://github.com/umn-earth-surface/gFlex.
An analytical approach to solving the flexure equations is computed by the superimposition of analytical solutions for flexure with a constant elastic thickness. Current implementations perform this superposition in the spatial domain. This works both on uniform grids and arbitrary meshes.
- The analytical solution method for an arbitrary mesh is useful for coupling Flexure with finite element models such as CHILD without requiring any regridding.
- Analytical solutions by superposition are an N2 problem, making this method become increasingly problematic for larger grids and numbers of nodes.
- Analytical solutions are computed based on sets of point loads at nodes or centers of cells, so they will fail and show too much isostatic response if the cells become larger than a modest fraction of a flexural wavelength; this is because at too large of a grid size, the approximation of summing immediately adjacent loads breaks down.
- Analytical solutions work only with approximations of constant elastic thickness.
Analytical solutions using spectral techniques are not yet implemented.
For the numerical implementation, gFlex computes a finite difference solution to the flexure equations for a lithospheric plate of nonuniform (or uniform, if one so desires) elastic thickness via a thin plate assumption. It uses the UMFPACK direct solvers to compute the solutions via a lower-upper decomposition of a coefficient matrix. The coefficient matrix in the 1D case is a pentadiagonal sparse matrix that is trivial to generate. In the 2D case, 2D grid is reordered into a 1D vector. UMFPACK solution routines are then able to copmute solutions to flexure in around a second or less, though sometimes up to a minute for very large grids. Iterative solution methods may also be used, but are not tested.
Some advantages of using a numerical (rather than analytical) solution are that:
- This method permits spatial variability in lithospheric elastic thickness. This allows the use of real maps of elastic thickness in models or synthetic maps of elastic thickness variability to test hypotheses.
- This rapid solution once the coefficient matrix is built makes this method a good choice for numerical models that require frequent updating of flexural deformations of the lithosphere.
- The rapid solution technique likewise allows efficient calculations of mixed sediment and/or water loading, or water loading with onlap and offlap, such that a constant fill density cannot be assumed and solutions must be produced iteratively.
- This model will not break down when grid cell sizes are increased greatly, as will happen for superposition of analytical solutions if the grid cells become too large relative to a flexural wavelength for the point load assumption to hold true.
Key physical parameters and equations
The largest component of gFlex is a solution to the PDE for lithospheric flexure in 2 dimensions. Our finite difference solutiosn follow van Wees and Cloetingh (1994):
- Flexure was developed first in MATLAB (Spring / early Summer 2010) and then in Python with Numpy, Scipy, and Matplotlib (translated October 2010).
- As of October 2011, Flexure became IRF- and CMT-compliant, and it was coupled to the landscape evolution model CHILD for the Fall 2011 CSDMS meeting. (Abstract and presentation are here, though I haven't dared to watch myself.)
- Current work is being done to improve boundary condition handling and the speed of finite difference solutions to constant elastic thickness problems (Early 2012).
- The project was revived in fall 2014m renamed gFlex, and brought to completion in late winter 2015.
|Nr. of publications:||2|
|Featured publication(s)||Year||Model described||Type of Reference||Citations|
| Wickert, A.D.; 2012. GFlex, version 0.6.. , , . |
|2012||GFlex||Source code ref.||--|
| Wickert, A. D.; 2016. Open-source modular solutions for flexural isostasy: gFlex v1.0. Geoscientific Model Development, 9, 997–1017. |
|See more publications of GFlex|
Version 0.9 (development) and higher
Instructions for the current version of gFlex are available from the README.md file at: https://github.com/umn-earth-surface/gFlex
(Modified from instructional emails)
This is the old version of Flexure".
You need to have python, numpy, scipy, and matplotlib installed to use flexure.
The version 0.1 structure (prone to change) consists of two main files.
- flexcalc.py is a python module which contains all of the functions needed to execute flexure.
- flexit.py is the frontend that, via optparse, gives you the ability to specify inputs and outputs and run flexural solutions via an interactive command-line interface.
In addition to these files, version 0.1 comes with some basic test loads and elastic thickness maps.
- Trial run
For the basic functionality on the first runthrough, you navigate to the directory with the flexit.py and type something in like:
python flexit.py -vcrp --dx=20000 --Te=Te_sample/bigrange_Te_test.txt --q0=q0_sample/top_bottom_bars.txt
You select which of the sample Te and q0 files you want. Andy_output.png is the flexural response (variable=w) output from:
python flexit.py -vcrp --dx=20000 --Te=Te_sample/bigrange_Te_test.txt --q0=q0_EW_bar.txt
The flags are explained in the help file (python flexit.py -h), as are other options for running the code. Basically you will want to run "-c" the first time, but not again unless you are going to redo the coefficient matrix (i.e., use a different pattern of elastic thickness). This calculation can take a long time for large grids, so you will want to store these files. When running without making a coefficient matrix, unless you're using the default coefficient matrix name (as we do above), you will have to specify its location with "--coeff-file=NAME". For example:
python flexit.py -vrp --Te=Te_sample/bigrange_Te_test.txt --q0=q0_sample/top_bottom_bars.txt --coeff_file=coeffs.txt
If all else fails
Have you typed:
python flexit.py -h
? This gives you all of the in-program help information.
Feel free to email Andy Wickert for anything related to this model (see contact info in box at top).
The coefficient matrix for the 2D finite difference solution with variable elastic thickness requires a map of elastic thicknesses in *.txt / ASCII format. This elastic thickness map must be two cells wider on each side than the map of loads; this is because the finite difference solution must "look" two cells in every direction. It also requires the specification of several parameters, including:
- Young's modulus (defaults to 1011 Pa)
- Poisson's ratio (defaults to 0.25)
- Mantle density (defaults to 3300 kg/m3)
- Density of infilling material (defaults to 0 kg/m3)
This outputs an ASCII sparse matrix file (Matrix Market *.mtx format).
The flexural solution requires the ASCII file for the sparse coefficient matrix generated above and an imposed array of loads (also ASCII), along with the specification of input and output file names.
The coefficient matrix creator writes a *.mtx sparse matrix ASCII file that is used in the direct solution. This matrix is characteristic to a given pattern of elastic thickness, and therefore can be reused if elastic thickness does not change.
The real solver outputs an ASCII grid of deflections due to the load. This is the output that is of scientific interest and/or useful to plug into other modules (e.g., for flexural subsidence).