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.

AC1: 'Reply on RC1', David Hindle, 06 Jul 2022
For the final response to reviewer one, we refer to our first comment and reply, but add the following for clarity.
1) As we stated, the paper actually has numerical examples which can determine accuracy by first order principles ("raft" can be deduced from first principles and force balance and broken plate can be compared directly to an analytical solution to an end loaded plate). Both are within small decimals of a percent of the numerical solution. We will however edit the text to make the importance of this point more clear.
2) The comparison to a 2D model was (and remains for us) an issue. It is certainly not clear that a 2D solution to a linear elastic problem is "correct" in any formal sense. It is a completely different approach to the problem, and although 2D, uses a 2nd order formulation. Mathematically, this cannot be the same as a 4th order equation, even if the 4th order equation is derived from a linear elastic starting point.
However, after some thought, we have produced the simplest and most fundamental comparison of the differential part of the equations, and in doing so removed all the additional geometric complexity of 2D (which has nothing to do with the similarity of 1D and 2D results). We have thus produced a simple, elastic beam model without a hydrostatic advective foundation (i.e. the term k u is eliminated leaving only the elastic beam equilibrium). For simplicity, we used FENICS and a modification of the 2D linear elastic tutorial, creating a beam with a "step" in its middle (its thickness reduced by 50% at a single point). The beam is selfloaded by its own weight and clamped at both ends. Hence, we have a step discontinuity in elastic thickness, modelled in 2D, using the linear elastic formulation. This geometry in 2D representation is on a single line, and "instantaneous". In 1D, the "step" must be spread over one single grid spacing. Hence there is always a tiny difference in formulation due to the nature of the approximations. We find a difference of ~34% depending on beam thickness, load and so on. However, large deformations are ruled out a priori, because both methods use the small strain approximation of linear elasticity (the 1D equation is actually less sensitive to this than the 2D finite element equivalent). We also note that no solution exists for such a problem with the whole station 1D discretisation. We would therefore propose to add this model comparison to the appendix, with a technical description, appropriate equations and figure, and reference it in the text, with perhaps a small amount of background.
3) Regarding 2D or 1D representations, the evidence of the simple comparison of the 2D finite element solution above is that there are relatively small differences. We have also noted that empirical comparisons show that even in cases of detached, separate, tilted blocks, the 1D comparison performs very well. It actually remains to test whether a 2D model could do the same. If not, we would strongly suggest the problem lies more with the 2D formulation, and the associated additional geometric and parametric complexity, particularly the implementation of the hydrostatic advection term, than anything else.
The remaining specific comments by the reviewer are addressed in our first response. We will modify the text as stated there to respond to them. We will add the Smoluchowski reference at an appropriate place in the text.
We thank Professor Schmalholz for his time reviewing our manuscript.

AC1: 'Reply on RC1', David Hindle, 06 Jul 2022

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.

AC2: 'Reply on RC2', David Hindle, 06 Jul 2022
For the final response to reviewer two, we refer to our first comment and reply, but add the following for clarity.
The principle points of reviewer two were dealt with in our earlier reply.
We would propose to make the modifications contained therein.
We thank Professor Colli for his time reviewing our manuscript.

AC2: 'Reply on RC2', David Hindle, 06 Jul 2022

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.

AC3: 'Reply on RC3', David Hindle, 06 Jul 2022
For the final response to reviewer three, we refer to our first comment and reply, but add the following for clarity.
The issue of elastic 2d models has been discussed in our final response to reviewer 1. We propose to include a 2D numerical version of a "stepped" beam problem compared to a half station equivalent 1D solution. Differences are small. The two equations (2nd order and 4th order) are actually not strictly mathematically comparable. It is widely assumed in the literature that they are equivalent however. We thank Professor Beaumont for his time reviewing our manuscript.

AC1: 'Reply on RC1', David Hindle, 06 Jul 2022
For the final response to reviewer one, we refer to our first comment and reply, but add the following for clarity.
1) As we stated, the paper actually has numerical examples which can determine accuracy by first order principles ("raft" can be deduced from first principles and force balance and broken plate can be compared directly to an analytical solution to an end loaded plate). Both are within small decimals of a percent of the numerical solution. We will however edit the text to make the importance of this point more clear.
2) The comparison to a 2D model was (and remains for us) an issue. It is certainly not clear that a 2D solution to a linear elastic problem is "correct" in any formal sense. It is a completely different approach to the problem, and although 2D, uses a 2nd order formulation. Mathematically, this cannot be the same as a 4th order equation, even if the 4th order equation is derived from a linear elastic starting point.
However, after some thought, we have produced the simplest and most fundamental comparison of the differential part of the equations, and in doing so removed all the additional geometric complexity of 2D (which has nothing to do with the similarity of 1D and 2D results). We have thus produced a simple, elastic beam model without a hydrostatic advective foundation (i.e. the term k u is eliminated leaving only the elastic beam equilibrium). For simplicity, we used FENICS and a modification of the 2D linear elastic tutorial, creating a beam with a "step" in its middle (its thickness reduced by 50% at a single point). The beam is selfloaded by its own weight and clamped at both ends. Hence, we have a step discontinuity in elastic thickness, modelled in 2D, using the linear elastic formulation. This geometry in 2D representation is on a single line, and "instantaneous". In 1D, the "step" must be spread over one single grid spacing. Hence there is always a tiny difference in formulation due to the nature of the approximations. We find a difference of ~34% depending on beam thickness, load and so on. However, large deformations are ruled out a priori, because both methods use the small strain approximation of linear elasticity (the 1D equation is actually less sensitive to this than the 2D finite element equivalent). We also note that no solution exists for such a problem with the whole station 1D discretisation. We would therefore propose to add this model comparison to the appendix, with a technical description, appropriate equations and figure, and reference it in the text, with perhaps a small amount of background.
3) Regarding 2D or 1D representations, the evidence of the simple comparison of the 2D finite element solution above is that there are relatively small differences. We have also noted that empirical comparisons show that even in cases of detached, separate, tilted blocks, the 1D comparison performs very well. It actually remains to test whether a 2D model could do the same. If not, we would strongly suggest the problem lies more with the 2D formulation, and the associated additional geometric and parametric complexity, particularly the implementation of the hydrostatic advection term, than anything else.
The remaining specific comments by the reviewer are addressed in our first response. We will modify the text as stated there to respond to them. We will add the Smoluchowski reference at an appropriate place in the text.
We thank Professor Schmalholz for his time reviewing our manuscript.

AC3: 'Reply on RC3', David Hindle, 06 Jul 2022
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  

635  261  23  919  11  10 
 HTML: 635
 PDF: 261
 XML: 23
 Total: 919
 BibTeX: 11
 EndNote: 10
Viewed (geographical distribution)
Country  #  Views  % 

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