09 Apr 2021
09 Apr 2021
Numerical solutions of the flexure equation
 ^{1}GeorgAugustUniversität Göttingen, Goldschmidtstr. 3, 37077 Göttingen, Germany
 ^{2}Université de Neuchâtel, Rue EmileArgand 11, CH2000, Neuchâtel
 ^{1}GeorgAugustUniversität Göttingen, Goldschmidtstr. 3, 37077 Göttingen, Germany
 ^{2}Université de Neuchâtel, Rue EmileArgand 11, CH2000, 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)

RC1: 'Comment on se202136', Stefan Markus Schmalholz, 13 May 2021
reply
In the manuscript titled “Numerical solutions of the flexure equation” the authors discuss two finite difference discretization schemes for the elastic flexure equation with spatially varying flexural rigidity. For spatially variable rigidity, they argue that the commonly applied discretization scheme, which they term “whole station method”, is not as accurate and numerically stable than another known, but not frequently applied, scheme, which they term “half station method”. The authors present examples with significantly varying rigidity for which the half station method provides more accurate solutions. From a strict numerical point of view, the comparison of the two discretization schemes is interesting. After reading the manuscript, I have three main questions concerning (1) the numerical accuracy of the presented numerical schemes, (2) the mechanical validity of the studied flexure equation and (3) the applicability of the studied equation to bending of the lithosphere. I have also several minor comments.
(1) Numerical accuracy: How does the accuracy of the whole and half station method increase with increasing numerical resolution? It would be useful if the authors provide at least one more comparison of the numerical results with an analytical solution for a variable thickness, for example for a linearly varying thickness or similar. I guess that the engineering literature has many of such analytical solutions. The authors could quantify how the error of the numerical solution, with respect to the analytical solution, decreases with increasing numerical resolution. Such “resolution or convergence test” could be done for both half station and whole station method and the numerical accuracy of the two methods could be compared.
(2) Mechanical validity: How accurate is the studied flexure equation with variable flexural rigidity compared to the correct full twodimensional (2D) solution for a plate with variable rigidity? The presented thinplate equation describing the firstorder characteristics of bending, or deflection, of an elastic plate (or beam) exposed to plate orthogonal loads is based on many simplifying assumptions. One of these assumptions is that shear strains are negligible so that the lateral variation of the vertical deflection (defining the curvature) can be related to plateparallel strains. The authors discuss cases for which the flexural rigidity varies laterally significantly. However, it is not demonstrated whether the thinplate equation with variable rigidity would still provide correct results if compared to a full 2D elastic solution for a bending elastic plate for which the thickness variation is accurately resolved (e.g. using a 2D finite element numerical solution). Also, the flexural rigidity can vary spatially in the same way, for example, either (1) due to a variation of thickness with constant Youngs modulus or (2) a variation in Youngs modulus with constant thickness. If these two cases are fully resolved with a 2D model, then the flexural response might be different, whereas the thinplate equation would provide the same result. I am missing in the manuscript an analysis which clarifies for what type of thickness variations the presented thinplate equation with variable rigidity provides mechanically accurate, or reasonable, results. If this “validity” issue is not clarified, one might end up with a numerically accurate solution for an equation which is mechanically considerably inaccurate.
(3) Applicability to bending of the lithosphere: The thinplate bending equation for elastic deformation is extremely useful to understand and quantify the firstorder response of the lithosphere to vertical loads, such as the formation of a flexural bulge associated to subduction zones or foreland basins. For the presented thinplate equation the considered thickness variations should be symmetric around the central, or neutral, line of the plate. Such “symmetric” thickness variations are likely not directly applicable to the lithosphere. It is hence not obvious how well a thinplate equation for elastic deformation will capture the flexural response of a lithosphere. There are studies which compare qualitatively the results of the thinplate flexure equation with results of full 2D numerical results of lithosphere flexure and these studies support the firstorder applicability of the thinplate equation with constant thickness (e.g. Funiciello et al. 2003, 10.1029/2001JB000896; Bessat et al. 2020, doi: 10.1093/gji/ggaa092; ). I do not remember studies which compared thinplate results for variable thickness with corresponding full 2D numerical results for a lithosphere exhibiting variable effective flexural rigidity. For example, Hetenyi et al. (2006; https://doi.org/10.1111/j.1365246X.2006.03198.x) show 2D flexural lithosphere models for which they calculate the spatial variation of the effective elastic thickness. The model of Hetenyi et al. (2006) could be an interesting application for the thinplate model with variable thickness presented by the authors, and the authors could show the applicability of their model to lithosphere flexure.
Minor comments
Line 131 : The elastic flexure equation including plateparallel compressive stresses has already been applied in preplate tectonics times to explain the buckling of the Earth crust (Smoluchowski, 1909).
Line 3640: This part is a bit misleading because it might imply that the flexural strength of the lithosphere is largely unknown, since the authors write that present estimates for h are between 0 and 100 km. However, the flexural response and hence the effective elastic thickness, depends on many parameters, such as the thermal structure or the width of the load (e.g. studies by Burov and/or Watts). This dependency may explain why published values of h vary significantly.
Line 41 to equation 1: The authors mention bending, but equation 1 is an equation for buckling, since P represents plateparallel stress which can cause buckling. Bending is usually referred to a situation where only forces, or loads, are considered which act orthogonal to the plate.
Line 48: The term “platewide stress” is not so clear. P could be referred to as a plateparallel compressive stress, e.g. when buckling is studied. The term “platewide” does not say anything about the orientation of this stress.
Line 54: What do the authors mean by “things”?
Line 5463: This part is a unclear with respect to physical units and the use of the term force. A force has units of [N]. A term consisting of the product of density (or density difference) times gravitational acceleration (g) times length (u) has units of [Pa] = [N/m^2]. In equation 1, q is used but in equation 2, q(x) is multiplied by density times g. So q and q(x) have fundamentally different units. In line 57 suddenly a term of density difference times g is termed a force.
Line 6778: This part is confusing. The physical presence of a load, such as a seamount, is discussed which might affect the infill space. So why not discuss how a seamount would affect the total thickness of the plate and hence the flexural rigidity?
Equation 5: From equation 4 to equation 5: when D=0, where did rho_F go?
Line 155158: For the “whole station” method the thickness variation should be simply spread over more numerical grid points. Also, with a higher numerical resolution, the whole station method can model the same thickness variation over the same length but requiring higher resolution. But saying that the whole station method is unable to simulate a "broken" plate is in my opinion incorrect, because the whole station method is able to provide results close to analytical solutions for a broken plate.
Line 194203: The accuracy of the whole station method can be increased by increasing the numerical resolution to correctly resolve the variations in h(x). The authors show that the half station method needs a smaller numerical resolution than the whole station method. But the whole station method can solve the same problems with the same accuracy, it just requires a higher resolution.
Best regards,
Stefan Schmalholz
Reference
Smoluchowski, M.: Über ein gewisses Stabilitätsproblem der Elastizitätslehre und dessen Beziehung zur Entstehung von Faltengbirgen, B. Acad. Sci. Cracov., 2, 3–20, 1909.
David Hindle and Olivier Besson
Model code and software
flexhs.f90 David Hindle, Oliver Besson https://github.com/davidhindle/flexure1dhs
David Hindle and Olivier Besson
Viewed
HTML  XML  Total  BibTeX  EndNote  

225  39  6  270  1  1 
 HTML: 225
 PDF: 39
 XML: 6
 Total: 270
 BibTeX: 1
 EndNote: 1
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1