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 two-dimensional (2D) solution for a plate with variable rigidity? The presented thin-plate equation describing the first-order 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 plate-parallel strains. The authors discuss cases for which the flexural rigidity varies laterally significantly. However, it is not demonstrated whether the thin-plate 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 thin-plate equation would provide the same result. I am missing in the manuscript an analysis which clarifies for what type of thickness variations the presented thin-plate 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 thin-plate bending equation for elastic deformation is extremely useful to understand and quantify the first-order 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 thin-plate 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 thin-plate equation for elastic deformation will capture the flexural response of a lithosphere. There are studies which compare qualitatively the results of the thin-plate flexure equation with results of full 2D numerical results of lithosphere flexure and these studies support the first-order applicability of the thin-plate 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 thin-plate 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.1365-246X.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 thin-plate model with variable thickness presented by the authors, and the authors could show the applicability of their model to lithosphere flexure.

Minor comments

Line 1-31 : The elastic flexure equation including plate-parallel compressive stresses has already been applied in pre-plate tectonics times to explain the buckling of the Earth crust (Smoluchowski, 1909).

Line 36-40: 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 plate-parallel 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 “plate-wide stress” is not so clear. P could be referred to as a plate-parallel compressive stress, e.g. when buckling is studied. The term “plate-wide” does not say anything about the orientation of this stress.

Line 54: What do the authors mean by “things”?

Line 54-63: 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 67-78: 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 155-158: 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 194-203: 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.

This manuscript focuses on the numerical solution of the thin-plate equation for the elastic flexure of the lithosphere with variable flexural rigidity. In particular, the authors compare two finite-difference 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 ad-hoc 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 156-158: 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.

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 whole-station and half-station 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 half-station 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 half-station 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 half-station 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 thin-plate 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.