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.

Highlights

  • We 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

  • We 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

  • The elastic bending of the lithosphere under crustal loads is a fundamental part of modern geodynamics, describing a swathe of processes including glacial isostatic adjustment (Walcott, 1972), foreland basin formation in compression (Beaumont, 1981), and the flexural response of the lithosphere to extension (Egan, 1992)

Read more

Summary

Introduction

The elastic bending of the lithosphere under crustal loads is a fundamental part of modern geodynamics, describing a swathe of processes including glacial isostatic adjustment (Walcott, 1972), foreland basin formation in compression (Beaumont, 1981), and the flexural response of the lithosphere to extension (Egan, 1992). A flexure dependent load term due to infill results, regardless of density variations, and requires an iterative solution due to the fact that regions with an imposed load cannot be simultaneously occupied by infill This form of the equation allows imposition of arbitrary forces to an elastic plate. There are a number of interesting and illustrative cases of elastic thickness variation where abrupt changes are required We note that both whole and half station discretisations contain for any grid point or node i, terms involving Di+1, Di, and Di−1 (see figure 4 and appendix B) meaning that a discontinuity in h(x) must be present across at least 3 grid nodes to take full effect. Current flexure models (Wickert, 2016) are based on the whole station (product rule) derivation of the numerical scheme, and should be revised

Conclusions
Half station discretisation
Whole station discretisation
355 Acknowledgements

Talk to us

Join us for a 30 min session where you can share your feedback and ask us any queries you have

Schedule a call

Disclaimer: All third-party content on this website/platform is and will remain the property of their respective owners and is provided on "as is" basis without any warranties, express or implied. Use of third-party content does not indicate any affiliation, sponsorship with or endorsement by them. Any references to third-party content is to identify the corresponding services and shall be considered fair use under The CopyrightLaw.