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: final response (author comments only)

RC1: 'Comment on se202136', Stefan Markus Schmalholz, 13 May 2021
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.

RC2: 'Comment on se202136', Lorenzo Colli, 27 May 2021
This manuscript focuses on the numerical solution of the thinplate equation for the elastic flexure of the lithosphere with variable flexural rigidity. In particular, the authors compare two finitedifference schemes: the commonly used “whole station method” and an already known but overlooked “half station method”. They conclude that the half station method is more stable and should be preferred to the whole station method.
Overall, I think that this manuscript is well written and its main conclusion regarding the superiority of the half station method is broadly correct and of practical use for studies of lithospheric flexure. However, the title of the manuscript is overly general with respect to its content. Moreover, I think the manuscript should characterize better the difference between the two methods in terms of convergence to the analytical solution and in terms of computational efficiency.
The whole station method is able to handle smooth variations of flexural rigidity. But, with the exception of idealized infinitely sharp discontinuities, any variation of flexural rigidity can be represented smoothly if the computational grid is sufficiently fine. So the whole station method should be able to handle even large spatial gradients of flexural rigidity, if the gradient per gridpoint is small enough. The difference between the two methods becomes one of degrees of accuracy and of computational efficiency. In many instances accurate solutions are possible even with the whole station method. On the other hand, there are geophysically interesting problems (some of which are presented by the authors) that feature exceedingly sharp variations of flexural rigidity, such that an accurate solution with the whole station method is impractical. This is an important distinction. The authors make some overly general statements that may lead to a fundamental misunderstanding of the validity of published results based on the whole station method.
Some more tests and examples featuring rapid variations of flexural rigidity should be carried out where analytical solutions are compared against numerical results using the two methods with coarser and finer discretizations. This would strengthen the manuscript considerably. On this note, while analytical solutions for arbitrary loads are hard to come by, it is relatively easy to “backward engineer” an analytical solution by plugging into equation 2 arbitrary adhoc functions of deflection and flexural rigidity, obtaining the corresponding causative load.
Regarding the title, if the authors want to keep the current one they should expand the scope of the manuscript considerably, including a discussion of other numerical methods (e.g., finite elements) and of 2D and 3D systems. Otherwise, they should modify the title to reflect better the contents of the manuscript.
Regards,
Lorenzo Colli
Minor comments and typos:
Line 56: The q in eq. 3 is not the same q of equations 1 and 2. Please use a different letter/symbol.
Line 131: “no third OR fourth derivative terms”
Line 134: typo in “arbitrary”
Lines 156158: please plot the resulting solution for the whole station method in a supplementary figure.
Figure 1: the restoring force from the mantle is k * u * Delta x. The missing u should appear in the equations at the bottom of the figure, too.

RC3: 'Comment on se202136', Christopher Beaumont, 09 Jun 2021
I have read the assessments posted by Stefan Schmalholz and Lorenzo Colli and agree that the points they raise need to be addressed. There is no point me repeating them in detail.
However, the points ‘2) Mechanical Validity’ made by Stefan Schmalholz are exactly my main concern. The authors assume the thin elastic plate flexure equation correctly characterizes flexure of thin uniform plates but do not discuss whether the underlying theory in deriving the equation remains valid when extending the use of the equation to plates with lateral variations in flexural rigidity (D). As pointed out by Stefan there are several simplifying assumptions made in deriving the equation for uniform thickness plates.
The authors of this paper need to add an assessment of when these conditions and assumptions remain valid when D varies laterally. It appears both the wholestation and halfstation methods give reasonably accurate answers when D varies ‘slowly’ along the plate in relation to the plate thickness (h). If so, is there any advantage to the halfstation method in such cases?
More interesting is the question of whether the derivation of the thin plate equation is valid when D varies ‘rapidly’, for example in the form of a step function (which is also relevant to the end loaded truncated plate problem). The halfstation method may be able to solve the equation accurately but is the equation itself correct? It seems that these questions need to be answered near the beginning of the paper to reassure readers.
The authors also discuss the McQeen and Beaumont 1989 paper on Tilted Block Basins and comment that the crustal/lithospheric blocks are rigid and just tilt with no flexure. I have already contacted David Hindle and given him a copy of a related 1986 MSc study by Paul Lloyd which addresses essentially the same tilted block basin problem with elastic flexure and in several different configurations. This thesis should be referenced and could be used for comparisons with models solved with the halfstation method if the lateral forcing can be converted to equivalent vertical forces or torques.
Overall, the Hindle and Besson paper is interesting from an academic point of view and will be strengthened if the ‘mechanical validity’ question is answered in some detail. When reviewing the model used by Paul Lloyd it occurred to me that for awkward geometries why would anyone try to use a thinplate model? If problems like the tilted bock basin one could be solved using a 2D or 3D elastic finite element approach 35 years ago, doing the same thing today must be relatively easy?
My recommendation is that the paper be published after modifications that adequately address the points raised by the reviewers.
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  

312  70  12  394  4  4 
 HTML: 312
 PDF: 70
 XML: 12
 Total: 394
 BibTeX: 4
 EndNote: 4
Viewed (geographical distribution)
Country  #  Views  % 

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