09 Apr 2021

09 Apr 2021

Review status: this preprint is currently under review for the journal SE.

Numerical solutions of the flexure equation

David Hindle1 and Olivier Besson2 David Hindle and Olivier Besson
  • 1Georg-August-Universität Göttingen, Goldschmidtstr. 3, 37077 Göttingen, Germany
  • 2Université de Neuchâtel, Rue Emile-Argand 11, CH-2000, Neuchâtel

Abstract. The 4th order differential equation describing elastic flexure of the lithosphere is one of the cornerstones of geodynamics, key to understanding topography, gravity, glacial isostatic rebound, foreland basin evolution and a host of other phenomena. Despite being fully formulated in the 1940’s, a number of significant issues concerning the basic equation have remained overlooked to this day. We first explain the different fundamental forms the equation can take and their difference in meaning and solution procedures. We then show how numerical solutions to flexure problems in general as they are currently formulated, are potentially unreliable in an unpredictable manner for cases where the coefficient of rigidity varies in space due to variations of the elastic thickness parameter. This is due to fundamental issues related to the numerical discretisation scheme employed. We demonstrate an alternative discretisation that is stable and accurate across the broadest conceivable range of conditions and variations of elastic thickness, and show how such a scheme can simulate conditions up to and including a completely broken lithosphere more usually modelled as an end loaded, single, continuous plate. Importantly, our scheme will allow breaks in plate interiors, allowing for instance, the creation of separate blocks of lithosphere which can also share the support of loads. The scheme we use has been known for many years, but remains rarely applied or discussed. We show that it is generally the most suitable finite difference discretisation of fourth order, elliptic equations of the kind describing many phenomena in elasticity, including the problem of bending of elastic beams. We compare the earlier discretisation scheme to the new one in 1 dimensional form, and also give the 2 dimensional discretisation based on the new scheme.We also describe a general issue concerning the numerical stability of any second order finite difference discretisation of a fourth order differential equation like that describing flexure where contrasting magnitudes of coefficients of different summed terms lead to round off problems which in turn destroy matrix positivity. We explain the use of 128 bit, floating point storage for variables to mitigate this issue.

David Hindle and Olivier Besson

Status: open (until 04 Jun 2021)

Comment types: AC – author | RC – referee | CC – community | EC – editor | CEC – chief editor | : Report abuse
  • RC1: 'Comment on se-2021-36', Stefan Markus Schmalholz, 13 May 2021 reply

David Hindle and Olivier Besson

Model code and software

flex-hs.f90 David Hindle, Oliver Besson

David Hindle and Olivier Besson


Total article views: 270 (including HTML, PDF, and XML)
HTML PDF XML Total BibTeX EndNote
225 39 6 270 1 1
  • HTML: 225
  • PDF: 39
  • XML: 6
  • Total: 270
  • BibTeX: 1
  • EndNote: 1
Views and downloads (calculated since 09 Apr 2021)
Cumulative views and downloads (calculated since 09 Apr 2021)

Viewed (geographical distribution)

Total article views: 261 (including HTML, PDF, and XML) Thereof 261 with geography defined and 0 with unknown origin.
Country # Views %
  • 1
Latest update: 17 May 2021
Short summary
By making a change to the way we solve the flexure equation, that describes how the earth's outer layer bends when it is subjected to loading by ice sheets or mountains , we develop new ways of using an old method from geodynamics. This lets us study the earth's outer layer by measuring a parameter called the elastic thickness, effectively how stiff and springy the outer layer is when it gets loaded, and and also how the earth's outer layer gets broken around its edges and in its interior.