on se-2021-36

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 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 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.

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 plateparallel 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?
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.