A corrected finite-difference scheme for the flexure equation with abrupt changes in coefficient

The fourth-order differential equation describing elastic flexure of the lithosphere is one of the cornerstones of geodynamics that is key to understanding topography, gravity, glacial isostatic rebound, foreland basin evolution, and a host of other phenomena. Despite being fully formulated in the 1940s, 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 as they are currently formulated are in general potentially unreliable in an unpredictable manner for cases in which 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 we 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 one-dimensional form and also give the two-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 wherein 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.

1 Introduction

The elastic bending of the lithosphere under crustal loads is a fundamental part of modern geodynamics describing a swathe of processes including glacial isostatic adjustment (Walcott, 1972), foreland basin formation in compression (Beaumont, 1981), and the flexural response of the lithosphere to extension (Egan, 1992). The mathematical theory as applied to vertical loads deflecting the Earth's lithosphere was originally proposed in the pre-plate tectonic era by Gunn (1943a, b, 1947), who was interested in the wider question of compensation of loads by isostatic balance. The original theory of elastic beam bending for engineering from which it was derived is attributed to Leonhard Euler and Daniel Bernoulli in the 1750s. Gunn (1943b) wrote his series of papers at the culmination of a many-decades-long debate between geodesists and geologists concerning how loads on the crust and lithosphere were compensated for by displacement of mantle material (see Barrell, 1914; Gilbert, 1889; Fisher, 1895). Geodesists had long favoured the idea that loads were all locally compensated for (Airy or Pratt isostasy). Gunn (1943b) was the first person to fully realise and formulate the necessary equations describing how a load can be compensated for over a much greater distance than its own width due to the elastic strength of the lithosphere. Before him, it was clear that several authors (Barrell, 1914; Gilbert, 1890) had very similar insights but were never able to quantitatively demonstrate them (see Watts, 2001, for a summary of the history of isostasy).

Figure 1 The force balance across a segment of a plate, Δx, and the derivation of the fourth-order differential equation describing elastic flexure. The segment is supported from below by mantle restoring forces (kΔx) and loaded from above by a distributed load qΔx. Shear stress V is equal to the first derivative of the moment M^′, which in turn is equal to the coefficient D multiplied by the second derivative of deflection of the plate u^′′. Over the small distance, between x−0.5Δx and x+0.5Δx, the flexural bending equation given above arises.

Gunn's work transformed our understanding of lithospheric mechanics, establishing how loads on the crust were balanced by elastic bending of the lithosphere as well as suggesting how this compensation would affect measured gravity anomalies. Subsequently, in the early post-plate tectonic era, Walcott (1970a, b, c) published a series of manuscripts on the question of the elastic thickness of the Earth's lithosphere. Elastic thickness (h) is the key parameter in the flexural coefficient D (the rigidity), which itself is related to the bending moment of the lithosphere (Fig. 1). A high value of D means less bending of the lithosphere under loading. Following Walcott (1970a), estimating the value of elastic thickness h, especially for continental lithosphere, became one of the most strongly debated topics in geodynamics (e.g. McKenzie and Fairhead, 1997; Watts, 2001; Audet and Mareschal, 2004; Kaban et al., 2018; Watts, 1992; Pérez-Gussinyé and Watts, 2005; Burov and Watts, 2006). At present, estimated values of h for the continents vary from 0 to 100 km, presumably due to a combination of factors related to thermal structure and tectonic history of a particular region of the lithosphere.

The fundamental equation describing the elastic bending of the lithosphere is our concern in this paper. Its original form presented to geologists by Gunn (1943b) remains unchanged. It describes the balance of forces (bending moment, vertical shearing forces, any added loads, and restoring forces from the buoyancy of mantle below) in an elastic plate resting on a fluid (inviscid) mantle, a so-called Winkler foundation (Fig. 1). Moment and shearing forces are converted to functions of the deflection of the plate (u(x) in this paper), resulting in a fourth-order differential equation in u given as

$$\begin{array}{}\text{(1)}& D{u}^{\prime \prime \prime \prime }+P{u}^{\prime \prime }-ku=q,\end{array}$$

where u(x) is the vertical deflection of an originally horizontal surface of the plate due to loading, $u^{\prime }={\displaystyle \frac{\mathrm{d}u}{\mathrm{d}x}}$ and so on for higher derivatives; q(x) is a term describing the applied load; P is a plate-wide stress, which is actually a uniform, compressive load in the plane of section of a plate; $k=({\mathit{\rho }}_{\mathrm{M}}-{\mathit{\rho }}_{\mathrm{F}})\cdot g$ is a constant allowing compensation of plate deflection by displacement of underlying “fluid” substrate (mantle); and $D={\displaystyle \frac{E{h}^{\mathrm{3}}}{\mathrm{12}(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})}}$ is the flexural rigidity in which E is the elastic modulus, h the elastic thickness of the lithosphere, and ν Poisson's ratio. For many solutions given in the literature D is assumed to be constant, meaning h does not vary along a plate's length. Should it be the case that h(x) is variable, however, then a modified general form of the equation above is

$$\begin{array}{}\text{(2)}& (D{u}^{\prime \prime }{)}^{\prime \prime }+P{u}^{\prime \prime }-ku=q.\end{array}$$

Equations (1) and (2) can be reformulated in the following ways. We begin by rewriting them explicitly to show the different components of the coefficients, while also setting P=0.

$$\begin{array}{}\text{(3)}& D{u}^{\prime \prime \prime \prime }-({\mathit{\rho }}_{\mathrm{M}}-{\mathit{\rho }}_{\mathrm{F}})gu={\mathit{\rho }}_{\mathrm{L}}gq\left(x\right)\end{array}$$

We note that the load term (ρ_M−ρ_F)gu corresponds to two body forces: ρ_Mgu is the restoring force (per unit length and width with units of Newtons) due to displacement of mantle by deflection of the plate, and ρ_Fgu is a force exerted by material assumed to infill any surface deflections of the plate below an arbitrary reference level, but it is also clear that without modification, any solution of this problem equally assumes that infill forces are removed wherever there is a positive deflection of the plate (Fig. 2). Such infill can range from nothing (empty basins) to water (oceanic cases, for instance) to sedimentary product (foreland basins) or a mixture of any of the above in various combinations. We also note that the load term ρ_Lgq(x) defines a separate load of potentially different density to that assumed for the infill.

Figure 2 Flexure equation and its physical significance. “Restoring” forces (constant multiplied by the plate vertical deflection, u(x)) have differing effects according to whether u(x) is positive when there is, by default, erosion of the plate to level zero; when negative, there may be infill of basins created by flexural subsidence. Mantle forces are always present and damp subsidence due to surface loading but equally damp uplift when flexure bends things above zero reference level.

Figure 3 The two types of flexure model derived from Eq. (2). (a) The “fixed topography” situation in which subsidence matches the prescribed topographic profile (associated with load density). Load thickness is then equal to topography Topo(x) minus subsidence u(x). (b) The fixed load case, equivalent to arbitrary forces. Any “load” may be applied, generating subsidence. Calculating infill of basins generated, should there be any, requires an iterative procedure since the amount of accommodation space must first be calculated explicitly and subsequently filled, hence creating more accommodation space.

Two versions of the equation can be developed from here. The first one (Fig. 3a), which has rarely been explicitly discussed or used in geodynamics, involves separating the infill load term from the mantle restoring force (see also Appendix A).

$$\begin{array}{}\text{(4)}& D{u}^{\prime \prime \prime \prime }-{\mathit{\rho }}_{\mathrm{M}}gu={\mathit{\rho }}_{\mathrm{L}}gq\left(x\right)-{\mathit{\rho }}_{\mathrm{F}}gu\end{array}$$

With the load term now entirely on the right-hand side of the equation, we see it consists of two parts with the infill part, ρ_Fgu (right hand side), dependent on the deflection u(x) for which we are solving. However, a fixed load ρ_Lgq(x) is also being applied, which is itself an arbitrary function of x. The physical interpretation of this depends on how this fixed load term is regarded. In most cases, it is probably assumed to be some kind of imposed crustal load, such as a thrust sheet or ice sheet, or, in oceanic cases, a seamount or volcanic island. Under these circumstances, it is equally clear that infill material cannot occupy space wherein the fixed load is applied unless the special circumstance applies that the top of the fixed load at some point lies below the original reference level, in which case a reduced accommodation space is available for infill material defined as the local sum u(x)+q(x). In general, implementing a solution to this form of the problem requires numerical methods, since arbitrary piecewise variations of load density may be required, and a solution will need to be iterative due to the dependence of the infill load term on u(x). It is also important to state that the fixed load term ρ_Fgq(x) is actually an applied force scaled for a particular load thickness. Hence, arbitrary “forces” can equally well be applied to the lithosphere, whatever their origin is assumed to be.

In general, it can be seen by examination that if D=0 (no flexural bending occurs and hence no marginal basin forms beyond the end of the load, so there is no basin infill load) the solution will correspond to “Airy isostasy” (i.e. an “iceberg” model) with

$$\begin{array}{}\text{(5)}& u\left(x\right)=\left\{\begin{array}{ll}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{L}}}{{\mathit{\rho }}_{\mathrm{M}}}}q\left(x\right)& q\left(x\right)\ne \mathrm{0}\\ \mathrm{0}& q\left(x\right)=\mathrm{0}\end{array}\right..\end{array}$$

An alternative development of the equation for solution is by division through by (ρ_M−ρ_F)g, giving

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{D{u}^{\prime \prime \prime \prime }}{({\mathit{\rho }}_{\mathrm{M}}-{\mathit{\rho }}_{\mathrm{F}})g}}-u={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{L}}}{({\mathit{\rho }}_{\mathrm{M}}-{\mathit{\rho }}_{\mathrm{F}})}}q\left(x\right).\end{array}$$

Here we see immediately that in contrast to Eq. (4), a direct, non-iterative solution is possible. This occurs because the physical meaning of the equation in this form is quite different. If we set D=0 once more, then

$$\begin{array}{}\text{(7)}& u\left(x\right)=\left\{\begin{array}{ll}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{L}}}{({\mathit{\rho }}_{\mathrm{M}}-{\mathit{\rho }}_{\mathrm{F}})}}q\left(x\right)& q\left(x\right)\ne \mathrm{0}\\ \mathrm{0}& q\left(x\right)=\mathrm{0}\end{array}\right..\end{array}$$

Equation (7) also represents a case of Airy isostasy. However, in this form, q(x) no longer represents a load thickness, but rather a load surface topography (Fig. 3b), for which the appropriate flexural compensation function u(x) is calculated. Hence, the difference q(x)−u(x) gives the resulting finite load thickness (in Fig. 3b we explicitly show Topo(x) as the “load” term and q(x) as a term derived from Topo(x)−u(x)). As can be seen from Eq. (7), however, a problem with this formulation arises due to the different density terms employed, in particular the different values of ρ_L and ρ_F. Only when ρ_F=ρ_L is a condition of Airy isostatic balance of a load of density ρ_L actually calculated, and hence the appropriate deflection u(x) and resulting load thickness q(x)−u(x). In analytical solutions of this equation, the only way to avoid this problem is by assuming ρ_F=ρ_L everywhere. Otherwise, for piecewise variable density terms, a numerical solution is required.

In summary, the general equation of flexure of the lithosphere can be formulated in two different ways. In the first, a flexure-dependent load term due to infill results, regardless of density variations, and requires an iterative solution due to the fact that regions with an imposed load cannot be simultaneously occupied by infill. Moreover, this form of the equation allows imposition of arbitrary forces to an elastic plate. A second and more common formulation of the problem describes the flexural subsidence required to support a particular surface topography. Analytical solutions to any flexure problem are generally unable to account for variable density of different load components or to differentiate between fill of basins created by flexural subsidence and removal of fill in any positively deflected regions. Hence, we now consider some numerical solutions to flexure problems.

Figure 4 Grid and finite-difference stencil for the solution of the flexure equation with variable elastic thickness, showing the grid elements involved in the discretisation of the problem pertaining to the value of u_i at grid node i. Variation in D (due to variation in h) at any node i is achieved across grid nodes $i-\mathrm{1},i,i+\mathrm{1}$. Hence, for abrupt changes in D and h, the value of h must be adjusted for at least three adjacent nodes in order to take full effect. The discretisation of the derivative of u requires the five nodes $i-\mathrm{2},\mathrm{\dots },i+\mathrm{2}$. Grid spacing dx gives the problem a physical dimension.

2 Numerical flexure solutions

For the most part, finite-difference methods have been applied to solve the flexure equation numerically in both the one-dimensional beam-type situation and for two-dimensional, thin elastic sheets. For the simplest case of constant flexural rigidity, D, the left-hand side of Eq. (3), for instance, is discretised as (see Fig. 4)

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{D}{\mathrm{d}{x}^{\mathrm{4}}}}\left(u_{i+\mathrm{2}}-\mathrm{4}{u}_{i+\mathrm{1}}+\left(\mathrm{6}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{M}}g\mathrm{d}{x}^{\mathrm{4}}}{D}}\right)u_i-\mathrm{4}{u}_{i-\mathrm{1}}+u_{i-\mathrm{2}}\right).\end{array}$$

Over the past 40 years, a number of numerical solutions to flexure problems were proposed. A lot of this effort was aimed at solving problems for the two-dimensional extension of the flexure equation to an elastic sheet with variable elastic thickness and hence flexural rigidity (Van Wees and Cloetingh, 1994). Van Wees and Cloetingh (1994) corrected what was probably the earliest attempt at a numerical solution with variable flexural rigidity (Bodine et al., 1981). Relatively few publications have dealt with the details of the numerical, one-dimensional flexure equation discussed here for both constant and variable elastic thickness cases. It should also be noted that since the Van Wees and Cloetingh (1994) initial publication, most of the succeeding work involving numerical solutions of the flexure equation has been based on their numerical derivation (Stewart and Watts, 1997; Van der Meulen et al., 2000; Govers et al., 2009; Braun et al., 2013; Wickert, 2016).

An early paper using a numerical solution (Stewart and Watts, 1997) illustrates this form of the one-dimensional numerical solution to the flexure equation with variable elastic thickness. Taking Eq. (2), for instance, the product rule of differentiation is applied prior to discretisation, and hence

$$\begin{array}{}\text{(9)}& D{u}^{\prime \prime \prime \prime }+\mathrm{2}{D}^{\prime }{u}^{\prime \prime \prime }+D^{\prime \prime }{u}^{\prime \prime }-ku=q.\end{array}$$

We note that this formulation of the problem, whilst mathematically correct, has a clear, unambiguous physical implication. Any function describing the variation of elastic thickness h(x) (and hence rigidity, D(x)) as a function of position must be continuous and at least twice differentiable (see Appendix C). It will often be possible to get solutions to this problem with other functions of D(x), but they are not consistent with the way the problem is posed in Eq. (9). The use of the product rule derivation of the problem extends to all of the aforementioned publications concerning the two-dimensional sheet-like problem as well (see passing from Eq. 3 to 7 of Van Wees and Cloetingh, 1994). The resulting finite-difference discretisation is given in Appendix B.

Perhaps surprisingly, there is an alternative and quite different method of discretising Eq. (2) available. This has been termed the “half-station” method (Cyrus and Fulton, 1966, 1968) and avoids a product rule derivation prior to discretisation entirely. Instead, Eq. (2) is directly transformed into a finite-difference approximation, replacing the derivatives both within and outside the brackets with finite-difference approximations to second derivatives (see Appendix B). The main effect of this is that no third or fourth derivative terms in u (and hence also first and second derivative terms in D) explicitly arise. Instead, the fourth-order nature of the differential equation as well as gradients in D are implicitly contained in the numerical scheme. Somewhat remarkably, this means that there is no restriction on the nature of the function D(x). Any piecewise, arbitrary variation of D(x) (and thus h(x)) is consistent with this discretised form of the equation. The half-station method is generalisable to two dimensions, meaning an alternative discretisation arises with Eq. (3) of Van Wees and Cloetingh (1994) as a starting point.

It must be noted here that both the half-station and whole-station methods are proven (Cyrus and Fulton, 1966, 1968) to be second-order accurate, finite-difference solutions to the flexure equation. Both converge to the analytical solution to the problem, effectively by a limited case of the Lax–Milgram theorem (Adams and Fournier, 2003, as demonstrated in Cyrus and Fulton, 1966). However, it is also clear that the whole-station method will never give a correct (convergent) solution for an abrupt variation (piecewise jump) in coefficient (Appendix C). Hence, for any case involving a piecewise linear or abrupt jump in the value of the coefficient, only the half-station method is convergent (see Appendix C).

Figure 5 Comparison of broken and unbroken plate loading simulated with “line load” (for parameters used, see Table A1): (a) the “broken” plate contains three node breaks of h=0.01 m; (b) continuous plate, same load. The maximum subsidence in ratio $\sim \mathrm{4}:\mathrm{1}$ with a small discrepancy due to the fact that the load is not a true “line” load but rather has a finite width equal to the grid spacing. This exactly corresponds to the analytical results of Gunn (1943b). Figure prepared using GMT v6.0.0 (Wessel et al., 2019).

3 Comparison of the different numerical schemes

We now address the question of the different behaviours of the numerical schemes with respect to variations of h(x). As we have already discussed, for the product rule version of the equation (so-called “whole-station” method; Cyrus and Fulton, 1968) the formal derivation of the numerical scheme actually requires a continuous, at least twice differentiable function of h(x). There are a number of interesting and illustrative cases of elastic thickness variation, however, in which abrupt changes are required. We note that both whole-station and half-station discretisations contain, for any grid point or node i, terms involving $D_{i+\mathrm{1}},D_i$, and D_i−1 (see Fig. 4 and Appendix B), meaning that a discontinuity in h(x) must be present across at least three grid nodes to take full effect. This was also commented upon by Wickert (2016) in the case of the whole-station scheme, but without detailed analysis of the results. A recognised difficulty with varying elastic thickness is how to find analytical solutions to which to compare numerical results. An indirect way of doing this, however, comes from Gunn (1943b), who showed analytically (see Eqs. 13 and 16, Gunn, 1943b) that for the same point load, the maximum deflection of a broken plate (i.e. one loaded at its end) is 4 times that of a continuous (infinite) plate equivalently loaded in its centre. Suppose we take a continuous plate and reduce the elastic thickness to zero over three nodes, creating an elastic “break”, and place a “point” load at a single node, directly to the left or right of the elastic break. We note that the load in this case is not at a single point, but instead applied over a finite width equal to the grid spacing used in the numerical scheme, so direct comparison to an analytical solution is difficult. However, we would expect the plate, when loaded just next to the elastic break, to act like the end-loaded or broken plate, whereas the same plate loaded equally but without an elastic break should behave like the continuous one, so the relative maximum subsidence of the two numerical cases should change by a factor of 4. Results of the experiments are shown in Fig. 5. The half-station method gives exactly the result expected, showing that it corresponds to an end-loaded plate when elastic breaks or discontinuities are present within a larger plate. Due to the severe violation of the conditions of continuity and twice differentiability in the function of h(x), the whole-station method gives no meaningful result and is unable to simulate a “broken” plate.

Figure 6 Isostatic “raft” model. With an effectively infinite thickness plate and a symmetrically loaded raft, which is detached at both ends by elastic breaks (h=0.01m for three nodes), the mass of mantle displaced (760 m thick layer) is almost exactly equal to the mass of the applied load, demonstrating isostatic balance without flexure (for parameters used, see Table A1). Figure prepared using GMT v6.0.0 (Wessel et al., 2019).

Figure 7 Tilted block model directly equivalent to the force balance model of McQueen and Beaumont (1989) but also incorporating flexure. Calculations assume that basins are filled with infill material (red shading) and the plate surface is eroded to zero topography, removing material (blue shading). For parameters used, see Table A1. (a) 100 km block length. Blocks are visually close to rigid and tilted (b) 200 km blocks with identical loading, which undergo substantial bending (∼5 times that of the 100 km block). Elastic bending increases as the plate is more strongly held down by the greater length over which mantle resistance forces can act. It is important to note, however, that the plate segments have unconstrained boundaries and are held in place only by their interaction with the mantle. Figure prepared using GMT v6.0.0 (Wessel et al., 2019).

Another interesting case to test the numerical schemes is that of what we can term an “isostatic raft” (Fig. 6). In this case, we simulate an effectively infinitely stiff plate (h≥500 km) with a central region of length ∼200 km, which is bordered at each end by an elastic break. The “raft” is loaded evenly across its centreline by a rectangular-shaped load of width ∼40 km. The relative dimensions are chosen purely to illustrate the point. If the plate segment were truly infinitely stiff, it would undergo no bending at all, and the load mass applied would be compensated for by escape of an equal mass of mantle substrate. The half-station method in this case produces an almost constant subsidence of the plate segment and causes a tiny amount of flexural bending. The resulting mass difference between displaced mantle and load is <0.2 %, showing the expected “raft-like” response. Again, for such a case, the whole-station method produces a spurious result.

A geological application of the raft analogue arises when we consider tilted crustal blocks formed in compression (McQueen and Beaumont, 1989). This refers to the concept of short segments of crust and lithosphere bounded by basement-transecting faults. McQueen and Beaumont (1989) initially created a simple force balance model wherein horizontal stress across a fault-bounded block generates a moment, which turns and tilts the block against the resisting force of the mantle and is potentially augmented by the effects of erosion and sedimentation on the tilted block surface. This idea was principally used to explore the amount of compressive stress required to “break” the lithosphere in plate interiors; however, a corollary of it was to explain subsidence and basin formation as also due to the block-tilting process. The model treats the block as completely rigid by default and assumes that the horizontal compressive force is responsible for the tilting, neglecting the effects of “self-loading” due to one block overriding another and also flexural bending induced in the block.

Whilst our flexure model cannot be directly related to horizontal compressional stresses potentially involved in breaking the lithosphere, it is trivial to produce a succession of adjacent crustal segments by placing elastic breaks across a plate, creating isolated segments of the desired dimensions. By loading each segment at or near its end, thus creating a turning force, and noting that any load can be treated as an arbitrary system of forces which arise for many different reasons, we produce a result similar to that of McQueen and Beaumont (1989) but which also takes into account the flexural bending in the segments. Figure 7 shows two situations of identically loaded blocks with length 100 and 200 km. The longer blocks undergo substantially more bending as a result of the mantle resisting force being spread over a greater length and consequently holding the plate down more firmly, allowing it to bend elastically to a greater degree under loading. We note that a flexural model with a plate containing elastic breaks effectively parameterises lithospheric structure in terms of changes in elastic thickness. In the case of the tilted block model, this parameterisation can be thought of as the net effect of bounding faults, with the applied load characterising any combination of the possible forces acting on a block (e.g. moment on the block due to horizontal stress, friction on the fault resisting tilting, self-loading due to one block overriding another). In our illustrative models here, we have used 4 km thick, 15 km wide, distributed loads, corresponding to a net force of 1.6×10¹² N. This produces relatively small uplifts at block corners, although clearly, the major component of the uplift of block corners is likely to be due to steady transport of basement faults bounding adjacent tilted blocks due to shortening. The subsidence induced in basins, by contrast, can be more directly related to the response to loads on the crust. In the case of the Laramide orogeny, for instance, maximum sedimentary thicknesses in the associated basins are ∼2–4 km (Hagen et al., 1985), which is quite close to the ∼3 km subsidence under the load in our models.

Figure 8 Piecewise linear variation of h(x) (for parameters used, see Table A1) with a (a) relatively gentle gradient, for which the whole-station and half-station schemes are in good agreement. (b) Sharper gradient for which there is a substantial difference between the whole-station and half-station methods. Figure prepared using GMT v6.0.0 (Wessel et al., 2019).

All the preceding cases concern situations in which the lithosphere is modelled as segmented or broken. In many cases, however, we consider elastic thickness to vary more steadily (e.g. Van der Meulen et al., 2000; Stewart and Watts, 1997). In such cases, the most straightforward spatial variation of h(x) is described by a piecewise linear function interpolated between a few points of fixed value. We find that in cases in which the gradient of the imposed linear change is not too sharp, the whole-station approximation can return reasonable results, including in cases in which h→0. Below a certain threshold, however, the error (which we take as the difference to the half-station solution) quickly reaches ≥2 %, which for studies fitting flexural curves to gravity signals, for instance (e.g. Stewart and Watts, 1997), will be critical. Figure 8 shows two cases with differing gradients in h(x). While the gentle gradient yields a difference of 6 m maximum subsidence (∼0.5 %) the sharper gradient reaches 340 m (∼30 %). As the gradient in h(x) increases still further, the whole-station scheme will ultimately reach a point at which it returns no solution at all, whilst the half-station method is stable for any combination of loads and variations of h(x).

4 Other numerical issues

A more general mathematical issue concerns the positivity of any numerical solution to a fourth-order differential equation of this kind. As can be verified (see Appendix B), all discretisations of the flexure equation using second-order finite-difference approximations will yield an identical set of linear equations when the value of h is a constant. The discretised form then becomes that shown in Eq. (8). It can also be seen that for constant values of D at least, for each line, the sum across the columns ($\sum _j{a}_{i,j}$) is ρ_Mg, the mantle restoring force, since all other terms involving D sum to zero. As a result, the residual term due to the mantle restoring force is, according to the maximum principle (Axelsson, 1994), necessary for maintaining the positivity of the system of equations represented by Eq. (8). Hence, the composition of the main diagonal, which itself consists of a sum of two terms, ${\displaystyle \frac{\mathrm{6}D}{\mathrm{d}{x}^{\mathrm{4}}}}$ and ρ_Mg, becomes of critical importance. This is due to the issue of round-off, whereby the capacity of 64 bit representations of numbers to sum terms with large contrasts in magnitude leads to the smaller term being partly or entirely lost as the maximum number of significant figures available in arithmetic operations (approximately 15) is exceeded. In the large term, D will vary as a function of h³, where elastic thickness h may reach values of 100 km on Earth (Kaban et al., 2018) and possibly even 300 km on Mars (Thor, 2016). Grid spacing dx requires values of ∼100 m or less to ensure convergence and also to allow reasonable resolution in the representation of loads. It should be noted that currently, the highest-resolution, public, and globally available topographic databases, SRTM (Earth Resources Observation And Science (EROS) Center, 2017) and ASTER GDEM (NASA/METI/AIST/Japan Spacesystems And U.S./Japan ASTER Science Team, 2009), are both on 1 arcsec (∼30 m) grids, whilst TANDEM-X data (DLR – German Aerospace Center, 2018) are relatively freely accessible to scientists on a 0.4 arcsec (∼10–12 m) grid. Hence, ${\displaystyle \frac{\mathrm{6}D}{\mathrm{d}{x}^{\mathrm{4}}}}$ will quite conceivably be of magnitude 10¹⁹ or more in “real” problems in geodynamics. The small term ρ_Mg will always be of order 10⁴. Under such conditions, round-off will lead to an “effectively singular” matrix and the numerical problem will fail. One work-around is relatively easily available. The use of quadruple-precision (128 bit) representation of floating-point numbers allows >30 significant figures to be taken account of in arithmetic operations. Although this costs additional memory and some speed, it ensures that any conceivable problem of flexure with “real-world” dimensions and parameters will be correctly dealt with by the numerical algorithm.

Besides the general issue described above, the additional term for the plate-wide stress P (see Eq. 2 and Appendix B) also has the potential to cause problems for the numerical solution. In particular, P, which has always been assumed to be constant throughout a plate, will interact with regions of variable elastic thickness in a potentially problematic way. For cases with elastic breaks, where D→0, for instance, a constant value of P will often lead to failure of the numerical solution. In such cases, it is probably reasonable to set P→0 on the three nodes of the discontinuity, thereby treating these as if they were an infinitely thin fault. Because the nodes either side of the break have the normal value of P applied, the continuity of P is respected to some degree.

5 Implications

The half-station method of discretisation we have presented here is clearly able to deal with a complete set of possible variations of elastic thickness in the flexure equation. The whole-station method (applying the product rule first), by contrast, is unable to be relied upon to do so. Perhaps the most worrying aspect of the whole-station method is that in some circumstances, it will give results that appear plausible but are in fact in error by significant amounts (Fig. 8b) and actually represent solutions among the set of transitional variations of h just before the method fails completely. The principal reason for the behaviour of the whole-station discretisation is the fact that, as posed, a condition of the equation is a twice differentiable, continuous function of D. Were we to apply such functions, the whole-station discretisation would perform safely. However, it is also known that the half-station method works with smoothly varying functions of D, and in general, the errors associated with the half-station method are always smaller than those of the whole-station method (Cyrus and Fulton, 1968). It is also unjustifiable to impose any such constraint on the nature of variations of elastic thickness of the lithosphere. Consequently, it appears clear that the finite-difference discretisation of the flexure equation should be carried out using the half-station method.

The wider application of the half-station method to many other differential equations with variable coefficients was originally noted by Cyrus and Fulton (1966, 1968). The specific application of it to the one-dimensional wave equation and the general rarity of its use have also been discussed by Langtangen (2016, p. 44). It seems clear that this form of finite-difference discretisation, which allows arbitrary piecewise variations of coefficients, is a potentially significant and generally overlooked method for the wider spectrum of the physical sciences.

The use of numerical solutions for the flexure equation covers many aspects of geodynamics. On the one hand, the determination of the elastic thickness of the lithosphere can be done using a forward modelling approach with flexure models used to match gravity data and topography (e.g. Walcott, 1970a; Karner and Watts, 1983; Watts, 1992; Stewart and Watts, 1997). In such cases it may well be necessary to look for solutions incorporating variable elastic thickness, especially around mountain fronts in foreland basins. Flexure models may also be used to study the dynamics of past flexural events (Burkhard and Sommaruga, 1998; DeCelles and Giles, 1996; Horton and DeCelles, 1997; Beaumont, 1981; Hagen et al., 1985; Hindle and Kley, 2021) for which they are often used to model subsidence patterns and explain basin formation. In many of these cases too, variable elastic thickness is likely to need taking account of. Increasingly, topics relating to global sea level rise and the melting of the polar ice caps will demand high-resolution models of flexural responses, which may require taking account of changes in elastic thickness of the lithosphere. More generally, the issue of elastic breaks within continental lithosphere has yet to be substantially explored and could have significant consequences for topics such as intraplate seismicity and seismic hazard. In short, it seems very important to make such numerical approximations in as accurate a way as possible. Current flexure models (Wickert, 2016) are based on the whole-station (product rule) derivation of the numerical scheme and should be revised.

6 Conclusions

Despite a long history of use in the literature and an apparent sense of being work completed, in fact a host of problems arising from a simple numerical analysis of the discretised flexure equation have remained untouched. When we examine these, we find there are significant issues with the method of discretisation used. It is not advisable under any circumstances to use a product rule derivation of an equation of this type when the coefficient varies as a function of coordinate. Realistic models of natural variations in elastic thickness (and many other coefficients in many other equations arising in natural sciences in general) will require sharp changes in those coefficients to be taken account of. A product rule scheme cannot do this successfully, especially for fourth-order differential equations.

Fundamental problems relating to the nature of the system of linear equations arising from discretisations of the flexure equation have also gone unnoticed so far. For small grid spacings, something which will inevitably become increasingly common as computer power increases, the numerical solution will rapidly become unstable and fail unless a 128 bit floating-point representation is employed. If this is used, the problem will probably be avoided at grid spacings ≥1 m, but below this threshold, instability could once more arise quite easily. Although we have presented only one-dimensional problems in this paper, it is nevertheless clear that everything shown here extends to two-dimensional, thin elastic plate formulations as well. To this end, we give the two-dimensional half-station formulation and discretisation of the problem (Appendix B). We will discuss two-dimensional solutions using this scheme in forthcoming papers.

Appendix A: General aspects of the equation when solved numerically

We begin with a general formulation of the flexure problem with variable coefficient and specified load (not topography) as follows. This form of the problem requires an iterative solution. We give the form of the equation for the case in which the gravitational constant is negative, i.e. $g=-\mathrm{9.81}$, meaning load thickness q(x)>0 acts as a downwards force on the lithosphere.

Table A1 Parameters used for Figs. 5–8.

$$\begin{array}{}\text{(A1)}& \left(D\right(x)u^{\prime \prime }{)}^{\prime \prime }+P{u}^{\prime \prime }-{\mathit{\rho }}_{\mathrm{M}}gu={\mathit{\rho }}_{\mathrm{L}}gq(x)-{\mathit{\rho }}_{\mathrm{F}}gu\end{array}$$

u is the deflection of the plate at position x along its length. D(x) (the flexural rigidity) varies in space and is explicitly written as a function of x. The value of D(x) is given by $Eh(x{)}^{\mathrm{3}}/(\mathrm{12}\cdot (\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}))$, where E is the elastic modulus of the lithosphere, h(x) is the effective elastic thickness of the lithosphere and is the parameter in D(x) that varies in space, and ν is Poisson's ratio. P is a constant representing a plate-wide horizontal stress. The term ρ_M⋅g (left hand side) represents a restoring force due to displaced mantle. On the right-hand side of the equation, ρ_Lgq(x) is the imposed load term which is chosen arbitrarily and has a density ρ_L that can be set for whatever load is being modelled. However, this load term can also be thought of as representing any type of force loading the plate (for instance, forces across a fault resolved in the vertical direction or torques from horizontal loads on rigid blocks). ρ_Fgu is the load force due to “infill” of basins. However, due to its dependence on u, the term acts as a load when u is negative and a positive force (pushing or pulling the plate upwards) when u is positive. This upward pull can be thought of as a force due to erosional removal of material, and by default, the amount of erosion is equal to the value of u as if the uplifted segment of plate were eroded to 0 m above reference level (Fig. 2, main text). We may wish to make the density of eroded material different to that of infill, for instance a “crustal” density ρ_C. In this case, the full equation is dependent on the sign of u(x) and can be written as

$$\begin{array}{}\text{(A2)}& \left(D\right(x)u^{\prime \prime }{)}^{\prime \prime }+P{u}^{\prime \prime }-{\mathit{\rho }}_{\mathrm{M}}gu=\left\{\begin{array}{ll}{\mathit{\rho }}_{\mathrm{L}}gq\left(x\right)-{\mathit{\rho }}_{\mathrm{F}}gu& u\left(x\right)<\mathrm{0}\\ {\mathit{\rho }}_{\mathrm{L}}gq\left(x\right)-{\mathit{\rho }}_{C}gu& u\left(x\right)\ge \mathrm{0}.\end{array}\right.\end{array}$$

Equally, basin fill is assumed to fill basins completely to the same reference level. We note that different values for erosion and fill levels (even spatially variable and piecewise) can be implemented relatively easily with a numerical code. For analytical solutions to the problem, it is implicit that there is infill and erosion and that the density of all materials is the same. We also state again that the first aim of the iterative scheme is to separate regions filled with fixed load, where there is subsidence given by u(x) but clearly no accommodation space exists for infill, from the basins created outside the regions occupied by the load. With a numerical method, it is relatively easy to ensure that this is the case.

Appendix B: Discretisation schemes

B1 Half-station discretisation

We apply second-order finite-difference operators simultaneously for both second derivatives inside and outside the brackets in Eq. (A1), which is something referred to as the half-station method (Cyrus and Fulton, 1968).

Hence, if

$$\begin{array}{}\text{(B1)}& f^{\prime \prime }\approx {\mathit{\delta }}^{\mathrm{2}}f=(f_{i+\mathrm{1}}-\mathrm{2}{f}_i+f_{i-\mathrm{1}})/\mathrm{d}{x}^{\mathrm{2}},\end{array}$$

where dx is the grid spacing, and i is the node number, then

$$\begin{array}{}\text{(B2)}& (D{u}^{\prime \prime }{)}^{\prime \prime }\approx {\mathit{\delta }}^{\mathrm{2}}(D{\mathit{\delta }}^{\mathrm{2}}u).\end{array}$$

We discretise the whole term in brackets first on a grid $i=\mathrm{1},\mathrm{\dots },N$:

$$\begin{array}{}\text{(B3)}& {\mathit{\delta }}^{\mathrm{2}}\left(D{\mathit{\delta }}^{\mathrm{2}}u\right)=\left(\right(D{\mathit{\delta }}^{\mathrm{2}}u{)}_{i+\mathrm{1}}-\mathrm{2}(D{\mathit{\delta }}^{\mathrm{2}}u{)}_i+(D{\mathit{\delta }}^{\mathrm{2}}u{)}_{i-\mathrm{1}})/\mathrm{d}{x}^{\mathrm{2}}.\end{array}$$

Then, substituting the terms in brackets and advancing the indices gives us

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{2}}\left(D{\mathit{\delta }}^{\mathrm{2}}u\right)& =\left(\right(D_{i+\mathrm{1}}(u_{i+\mathrm{2}}-\mathrm{2}{u}_{i+\mathrm{1}}+u_i))\\ & -\mathrm{2}\left(D_i\right(u_{i+\mathrm{1}}-\mathrm{2}{u}_i+u_{i-\mathrm{1}}\left)\right)\\ & +\left(D_{i-\mathrm{1}}\right(u_i-\mathrm{2}{u}_{i-\mathrm{1}}+u_{i-\mathrm{2}}\left)\right)/\mathrm{d}{x}^{\mathrm{4}}.\end{array}\end{array}$$

Collecting terms, we obtain

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{2}}\left(D{\mathit{\delta }}^{\mathrm{2}}u\right)& =(D_{i+\mathrm{1}}{u}_{i+\mathrm{2}}-\mathrm{2}(D_{i+\mathrm{1}}+D_i)u_{i+\mathrm{1}}\\ & +(D_{i+\mathrm{1}}+\mathrm{4}{D}_i+D_{i-\mathrm{1}})u_i\\ & -\mathrm{2}(D_{i-\mathrm{1}}+D_i)u_{i-\mathrm{1}}+D_{i-\mathrm{1}}{u}_{i-\mathrm{2}})/\mathrm{d}{x}^{\mathrm{4}}.\end{array}\end{array}$$

Discretising the remaining parts of the equation then gives

$$\begin{array}{}\text{(B6)}& \begin{array}{l}{\displaystyle \frac{D_{i+\mathrm{1}}}{\mathrm{d}{x}^{\mathrm{4}}}}{u}_{i+\mathrm{2}}\\ \\ +\left(-\mathrm{2}{\displaystyle \frac{(D_{i+\mathrm{1}}+D_i)}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{P_i}{\mathrm{d}{x}^{\mathrm{2}}}}\right)u_{i+\mathrm{1}}\\ \\ +\left({\displaystyle \frac{(D_{i+\mathrm{1}}+\mathrm{4}{D}_i+D_{i-\mathrm{1}})}{\mathrm{d}{x}^{\mathrm{4}}}}-\mathrm{2}{\displaystyle \frac{P_i}{\mathrm{d}{x}^{\mathrm{2}}}}-{\mathit{\rho }}_{\mathrm{M}}\cdot g\right)u_i\\ \\ +\left(-\mathrm{2}{\displaystyle \frac{(D_{i-\mathrm{1}}+D_i)}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{P_i}{\mathrm{d}{x}^{\mathrm{2}}}}\right)u_{i-\mathrm{1}}\\ \\ {\displaystyle \frac{D_{i-\mathrm{1}}}{\mathrm{d}{x}^{\mathrm{4}}}}{u}_{i-\mathrm{2}}\\ \\ =(q_o{)}_i+q(u_i),\end{array}\end{array}$$

where the two load terms, (q_o)_i and q(u_i), represent the static, fixed load and the iteratively calculated infill load, respectively.

If we gather all coefficients into a matrix A and form a matrix equation, the resulting system is of the form

$$\begin{array}{}\text{(B7)}& \mathbf{A}u=q\left(u\right),\end{array}$$

which is a non-linear series of equations in u. We reformulate this as a recursive matrix fixed-point problem, which we solve using a pentadiagonal matrix algorithm (Sebben and Baliga, 1995).

A similar procedure is used to discretise the specified topography formulation (Eq. 3).

B2 Whole-station discretisation

The whole-station discretisation begins from the result of applying the product rule to Eq. (B1), giving us

$$\begin{array}{}\text{(B8)}& D{u}^{\prime \prime \prime \prime }+\mathrm{2}{D}^{\prime }{u}^{\prime \prime \prime }+D^{\prime \prime }{u}^{\prime \prime }+P{u}^{\prime \prime }-{\mathit{\rho }}_{\mathrm{M}}gu={\mathit{\rho }}_{\mathrm{L}}gq\left(x\right)-{\mathit{\rho }}_{\mathrm{F}}gu.\end{array}$$

Discretisation involves applying second-order finite-difference schemes directly to all derivatives.

Hence,

$$\begin{array}{}\text{(B9)}& \begin{array}{l}D{u}^{\prime \prime \prime \prime }\approx {\displaystyle \frac{D_i}{\mathrm{d}{x}^{\mathrm{4}}}}(u_{i+\mathrm{2}}-\mathrm{4}{u}_{i+\mathrm{1}}+\mathrm{6}{u}_i-\mathrm{4}{u}_{i-\mathrm{1}}+u_{i-\mathrm{2}})\\ \\ \mathrm{2}{D}^{\prime }{u}^{\prime \prime \prime }\approx {\displaystyle \frac{(D_{i+\mathrm{1}}-D_{i-\mathrm{1}})}{\mathrm{2}\cdot \mathrm{d}{x}^{\mathrm{4}}}}(u_{i+\mathrm{2}}-\mathrm{2}{u}_{i+\mathrm{1}}+\mathrm{2}{u}_{i-\mathrm{1}}-u_{i-\mathrm{2}})\\ \\ D^{\prime \prime }{u}^{\prime \prime }\approx {\displaystyle \frac{(D_{i+\mathrm{1}}-\mathrm{2}{D}_i+D_{i-\mathrm{1}})}{\mathrm{d}{x}^{\mathrm{4}}}}(u_{i+\mathrm{1}}-\mathrm{2}{u}_i+u_{i-\mathrm{1}}).\\ \end{array}\end{array}$$

Adding the remaining terms from Eq. (B8) and explicitly writing to show the relationship to Eq. (B9), we have the following.

$$\begin{array}{}\text{(B10)}& \begin{array}{l}\left({\displaystyle \frac{D_i}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{(D_{i+\mathrm{1}}-D_{i-\mathrm{1}})}{\mathrm{2}\cdot \mathrm{d}{x}^{\mathrm{4}}}}\right)u_{i+\mathrm{2}}\\ \\ +\left(-\mathrm{4}{\displaystyle \frac{D_i}{\mathrm{d}{x}^{\mathrm{4}}}}-{\displaystyle \frac{(D_{i+\mathrm{1}}-D_{i-\mathrm{1}})}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{(D_{i+\mathrm{1}}-\mathrm{2}{D}_i+D_{i-\mathrm{1}})}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{P_i}{\mathrm{d}{x}^{\mathrm{2}}}}\right)u_{i+\mathrm{1}}\\ \\ +\left(\mathrm{6}{\displaystyle \frac{D_i}{\mathrm{d}{x}^{\mathrm{4}}}}-\mathrm{2}{\displaystyle \frac{(D_{i+\mathrm{1}}-\mathrm{2}{D}_i+D_{i-\mathrm{1}})}{\mathrm{d}{x}^{\mathrm{4}}}}-\mathrm{2}{\displaystyle \frac{P_i}{\mathrm{d}{x}^{\mathrm{2}}}}-{\mathit{\rho }}_{\mathrm{M}}\cdot g\right)u_i\\ \\ +\left(-\mathrm{4}{\displaystyle \frac{D_i}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{(D_{i+\mathrm{1}}-D_{i-\mathrm{1}})}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{(D_{i+\mathrm{1}}-\mathrm{2}{D}_i+D_{i-\mathrm{1}})}{\mathrm{d}{x}^{\mathrm{4}}}}+{\displaystyle \frac{P_i}{\mathrm{d}{x}^{\mathrm{2}}}}\right)u_{i-\mathrm{1}}\\ \\ +\left({\displaystyle \frac{D_i}{\mathrm{d}{x}^{\mathrm{4}}}}-{\displaystyle \frac{(D_{i+\mathrm{1}}-D_{i-\mathrm{1}})}{\mathrm{2}\cdot \mathrm{d}{x}^{\mathrm{4}}}}\right)u_{i-\mathrm{2}}\\ \\ =(q_o{)}_i+q(u_i)\end{array}\end{array}$$

The same procedure as for the half-station discretisation is employed to solve these equations. The iteration could be made more efficient.

B3 Half-station, two-dimensional discretisation

As has been discussed, all two-dimensional, thin elastic sheet type of solutions used up to the present day have been based on the same product rule derivation of the numerical scheme (whole-station). It is equally possible to apply a half-station derivation, however. We start from Eq. (3) of Van Wees and Cloetingh (1994), and assuming that ν is constant, we can write

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}(D{u}_{xx}{)}_{xx}& +(D{u}_{yy}{)}_{yy}+\mathit{\nu }((D{u}_{xx}{)}_{yy}+(D{u}_{yy}{)}_{xx})\\ & +\mathrm{2}(\mathrm{1}-\mathit{\nu })(D{u}_{xy}{)}_{xy}+P\mathrm{\Delta }u-{\mathit{\rho }}_{\mathrm{M}}gu\\ & ={\mathit{\rho }}_{\mathrm{L}}gq\left(x\right)-{\mathit{\rho }}_{\mathrm{F}}gu,\end{array}\end{array}$$

where $u_{xx}={\displaystyle \frac{{\partial }^{\mathrm{2}}u(x,y)}{\partial x^{\mathrm{2}}}}$ and so on. Defining partial finite-difference operators on a two-dimensional grid (i,j),

$$\begin{array}{}\text{(B12)}& \begin{array}{l}{\displaystyle \frac{\partial f}{\partial x}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{d}x}}(f_{i+\mathrm{1}/\mathrm{2}}-f_{i-\mathrm{1}/\mathrm{2}})\\ \\ {\displaystyle \frac{\partial f}{\partial y}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{d}y}}(f_{j+\mathrm{1}/\mathrm{2}}-f_{j-\mathrm{1}/\mathrm{2}}),\\ \end{array}\end{array}$$

and assuming $\mathit{\lambda }:=\mathrm{d}x=\mathrm{d}y$, Eq. (B11) can be discretised term by term in a way similar to the half-station method applied to the one-dimensional form. We have used maxima to derive the solution.

$$\begin{array}{}\text{(B13)}& \begin{array}{l}{u}_{i+\mathrm{2},j}\mathrm{2}{D}_{i+\mathrm{1},j}\\ +u_{i+\mathrm{1},j+\mathrm{1}}\left[-(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j+\mathrm{1}}\right.\\ \left.+(\mathit{\nu }+\mathrm{1})D_{i+\mathrm{1},j}+(\mathit{\nu }+\mathrm{1})D_{i,j+\mathrm{1}}-(\mathit{\nu }-\mathrm{1})D_{i,j}\right]+\\ u_{i+\mathrm{1},j}\left[(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j+\mathrm{1}}-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i+\mathrm{1},j}\right.\\ +(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j-\mathrm{1}}+(\mathit{\nu }-\mathrm{1})D_{i,j+\mathrm{1}}-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i,j}\\ \left.+(\mathit{\nu }-\mathrm{1})D_{i,j-\mathrm{1}}+\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}}{P}_{i,j}\right]+u_{i,j+\mathrm{1}}\left[(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j+\mathrm{1}}\right.\\ +(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j}-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i,j+\mathrm{1}}-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i,j}\\ \left.+(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j+\mathrm{1}}+(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j}+\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}}{P}_{i,j}\right]\\ +u_{i,j}\left[-(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j+\mathrm{1}}-\mathrm{2}(\mathit{\nu }-\mathrm{2})D_{i+\mathrm{1},j}-(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j-\mathrm{1}}\right.\\ -\mathrm{2}(\mathit{\nu }-\mathrm{2})D_{i,j+\mathrm{1}}+\mathrm{4}(\mathrm{3}\mathit{\nu }+\mathrm{5})D_{i,j}-\mathrm{2}(\mathit{\nu }-\mathrm{2})D_{i,j-\mathrm{1}}\\ -(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j+\mathrm{1}}-\mathrm{2}(\mathit{\nu }-\mathrm{2})D_{i-\mathrm{1},j}-(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j-\mathrm{1}}\\ \left.-\mathrm{8}{\mathit{\lambda }}^{\mathrm{2}}{P}_{i,j}-\mathrm{2}{\mathit{\lambda }}^{\mathrm{4}}{\mathit{\rho }}_{\mathrm{M}}g\right]+\\ u_{i+\mathrm{1},j-\mathrm{1}}\left[(\mathit{\nu }+\mathrm{1})D_{i+\mathrm{1},j}-(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j-\mathrm{1}}\right.\\ \left.-(\mathit{\nu }-\mathrm{1})D_{i,j}+(\mathit{\nu }+\mathrm{1})D_{i,j-\mathrm{1}}\right]+\\ u_{i,j-\mathrm{1}}\left[(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j}+(\mathit{\nu }-\mathrm{1})D_{i+\mathrm{1},j-\mathrm{1}}-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i,j}\right.\\ \left.-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i,j-\mathrm{1}}+(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j}+(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j-\mathrm{1}}+\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}}{P}_{i,j}\right]\\ +u_{i,j+\mathrm{2}}\mathrm{2}{D}_{i,j+\mathrm{1}}+\\ u_{i-\mathrm{1},j+\mathrm{1}}\left[(\mathit{\nu }+\mathrm{1})D_{i,j+\mathrm{1}}-(\mathit{\nu }-\mathrm{1})D_{i,j}-(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j+\mathrm{1}}\right.\\ \left.+(\mathit{\nu }+\mathrm{1})D_{i-\mathrm{1},j}\right]+\\ u_{i-\mathrm{1},j}\left[(\mathit{\nu }-\mathrm{1})D_{i,j+\mathrm{1}}-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i,j}\right.\\ +(\mathit{\nu }-\mathrm{1})D_{i,j-\mathrm{1}}+(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j+\mathrm{1}}-\mathrm{2}(\mathit{\nu }+\mathrm{3})D_{i-\mathrm{1},j}\\ \left.+(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j-\mathrm{1}}+\mathrm{2}{\mathit{\lambda }}^{\mathrm{2}}{P}_{i,j}\right]+\\ u_{i-\mathrm{1},j-\mathrm{1}}\left[-(\mathit{\nu }-\mathrm{1})D_{i,j}+(\mathit{\nu }+\mathrm{1})D_{i,j-\mathrm{1}}+(\mathit{\nu }+\mathrm{1})D_{i-\mathrm{1},j}\right.\\ \left.-(\mathit{\nu }-\mathrm{1})D_{i-\mathrm{1},j-\mathrm{1}}\right]+\\ u_{i,j-\mathrm{2}}\mathrm{2}{D}_{i,j-\mathrm{1}}+\\ u_{i-\mathrm{2},j}\mathrm{2}{D}_{i-\mathrm{1},j}\\ =\mathrm{2}{\mathit{\lambda }}^{\mathrm{4}}\left[\right(q_o{)}_{i,j}+q\left(u_{i,j}\right)]\end{array}\end{array}$$

Appendix C: Convergence, numerical comparisons

Here we give a general proof of the limitations of the whole-station discretisation as well as a comparison of a single case of the application of a 1D half-station problem to an (approximately) equivalent, 2D, finite-element, numerical solution.

C1 Proof of convergence

Consider the simplified beam equation:

$$\begin{array}{}\text{(C1)}& \left(D\right(x\left)u^{\prime \prime }\right(x){)}^{\prime \prime }+\mathit{\alpha }(x\left)u\right(x)=q(x)\end{array}$$

for $x\in \mathrm{\Omega }=(\mathrm{0},L)$, and with the boundary conditions

$$\begin{array}{}\text{(C2)}& u\left(\mathrm{0}\right)=u^{\prime }\left(\mathrm{0}\right)=u\left(L\right)=u^{\prime }\left(L\right)=\mathrm{0}.\end{array}$$

The following result is a simple application of the Lax–Milgram theorem (e.g. Adams and Fournier, 2003).

Theorem 1. Assume that

• D∈L^∞(Ω), and there is some D₀>0 such that D(x)≥D₀ a.e.,

• α∈L²(Ω), with α(x)≥0 a.e.

Then the problem in Eqs. (C1)–(C2) admits a unique solution in the Sobolev space $H_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{\Omega }\right)$.

This theorem shows that the function D can be very rough.

Now consider the following finite-difference operators.

$$\begin{array}{c}{\displaystyle }{\mathit{\delta }}^{\mathrm{1}}f\left(x\right)={\displaystyle \frac{f(x+\mathrm{d}x)-f(x-\mathrm{d}x)}{\mathrm{2}\mathrm{d}x}}\\ {\displaystyle }{\mathit{\delta }}^{\mathrm{2}}f\left(x\right)={\displaystyle \frac{f(x+\mathrm{d}x)-\mathrm{2}f\left(x\right)+f(x-\mathrm{d}x)}{\mathrm{d}{x}^{\mathrm{2}}}}\\ {\displaystyle }{\mathit{\delta }}^{\mathrm{3}}f\left(x\right)={\displaystyle \frac{f(x+\mathrm{2}\mathrm{d}x)-\mathrm{2}f(x+\mathrm{d}x)+\mathrm{2}f(x-\mathrm{d}x)-f(x-\mathrm{2}\mathrm{d}x)}{\mathrm{2}\mathrm{d}{x}^{\mathrm{3}}}}\\ {\displaystyle }{\mathit{\delta }}^{\mathrm{4}}f\left(x\right)={\displaystyle \frac{f(x+\mathrm{2}\mathrm{d}x)-\mathrm{4}f(x+\mathrm{d}x)+\mathrm{6}f\left(x\right)-\mathrm{4}f(x-\mathrm{d}x)+f(x-\mathrm{2}\mathrm{d}x)}{\mathrm{d}{x}^{\mathrm{4}}}}\end{array}$$

When numerically solving Eq. (C1), two different methods can be used.

• One is the half-station method in which Eq. (C1) is replaced by the difference equation.

  $$\begin{array}{}\text{(C3)}& {\mathit{\delta }}^{\mathrm{2}}\left(D\right(x\left){\mathit{\delta }}^{\mathrm{2}}u\right(x\left)\right)+\mathit{\alpha }\left(x\right)u\left(x\right)=q\left(x\right)\end{array}$$

• Another is the whole-station method in which Eq. (C1) is first developed as

  $$\begin{array}{}\text{(C4)}& \begin{array}{rl}D\left(x\right)u^{\prime \prime \prime \prime }\left(x\right)& +\mathrm{2}{D}^{\prime }\left(x\right)u^{\prime \prime \prime }\left(x\right)+D^{\prime \prime }\left(x\right)u^{\prime \prime }\left(x\right)\\ & +\mathit{\alpha }\left(x\right)u\left(x\right)=q\left(x\right)\end{array}\end{array}$$

  and then replaced by the difference equation

  $$\begin{array}{}\text{(C5)}& \begin{array}{rl}D\left(x\right){\mathit{\delta }}^{\mathrm{4}}u\left(x\right)& +\mathrm{2}{\mathit{\delta }}^{\mathrm{1}}D\left(x\right){\mathit{\delta }}^{\mathrm{3}}u\left(x\right)+{\mathit{\delta }}^{\mathrm{2}}D\left(x\right){\mathit{\delta }}^{\mathrm{2}}u\left(x\right)\\ & +\mathit{\alpha }\left(x\right)u\left(x\right)=q\left(x\right).\end{array}\end{array}$$

The results given in the paper of Cyrus and Fulton (1966) show that both methods are of the same order dx² when the function D is twice differentiable.

C2 Whole-station approach and regularity of D

Let us show that this formula makes sense only if D is at least twice differentiable. The whole-station discretisation can be written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{d}{x}^{\mathrm{4}}}}{D}_i(u_{i+\mathrm{2}}-\mathrm{4}{u}_{i+\mathrm{1}}+\mathrm{6}{u}_i-\mathrm{4}{u}_{i-\mathrm{1}}+u_{i-\mathrm{2}})\\ {\displaystyle }& {\displaystyle }+\mathrm{2}{\displaystyle \frac{D_{i+\mathrm{1}}-D_{i-\mathrm{1}}}{\mathrm{2}\mathrm{d}x}}{\displaystyle \frac{-u_{i-\mathrm{2}}+\mathrm{2}{u}_{i-\mathrm{1}}-\mathrm{2}{u}_{i+\mathrm{1}}+u_{i+\mathrm{2}}}{\mathrm{2}\mathrm{d}{x}^{\mathrm{3}}}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{D_{i-\mathrm{1}}-\mathrm{2}{D}_i+D_{i+\mathrm{1}}}{\mathrm{d}{x}^{\mathrm{2}}}}{\displaystyle \frac{u_{i-\mathrm{1}}-\mathrm{2}{u}_i+u_{i+\mathrm{1}}}{\mathrm{d}{x}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+{\mathit{\alpha }}_i{u}_i\\ {\displaystyle }& {\displaystyle }=q_i.\end{array}$$

Now consider the academic case L=1 with two different simple functions D:

1. $D\left(x\right)=\mathrm{1}/\mathrm{2}$ for $\mathrm{0}<x\le \mathrm{1}/\mathrm{2}$ and D(x)=1 for $x>\mathrm{1}/\mathrm{2}$, so D is discontinuous at $x=\mathrm{1}/\mathrm{2}$;

2. $D\left(x\right)=\mathrm{1}/\mathrm{2}+x$ for $\mathrm{0}<x\le \mathrm{1}/\mathrm{2}$ and D(x)=1 for $x>L/\mathrm{2}$, so D is continuous and piecewise affine.

For case 1, when dx→0 the term ${\displaystyle \frac{D_{i+\mathrm{1}}-D_{i-\mathrm{1}}}{\mathrm{2}\mathrm{d}x}}\to \mathrm{\infty }$ in the neighbourhood of $x=\mathrm{1}/\mathrm{2}$.

For case 2, when dx→0 the term ${\displaystyle \frac{D_{i-\mathrm{1}}-\mathrm{2}{D}_i+D_{i+\mathrm{1}}}{\mathrm{d}{x}^{\mathrm{2}}}}\to \mathrm{\infty }$ in the neighbourhood of $x=\mathrm{1}/\mathrm{2}$.

Figure C1 Comparison of a 2D linear elastic beam model with fourth-order, 1D, half-station approximation for a case with abrupt jumps in elastic thickness. The maximum difference between the two models is ∼3.5 % for a case with much larger strain than lithospheric examples.

Conclusion. In both cases it is impossible to verify the numerical accuracy of the whole-station discretisation. Moreover, this scheme leads to incorrect numerical solutions when D is not sufficiently regular.

Remark. In the half-station case, D^′(x) and D^′′(x) are not explicitly computed. So the formula (C3) works even when the function D is discontinuous in some places.

Remark. When the function D is at least twice continuously differentiable, the whole-station and the half-station methods give the same results.

C3 2D numerical comparison

In Fig. C1 we show a case in which only the half-station method is able to give a result. We demonstrate a simple, 25 m long elastic beam pinned at both ends with no substrate. Hence, the equation only concerns the elastic part of the flexure problem and therefore the elliptic part of the differential equation. The beam has a symmetrical shape, with a central section, L=15 m, half the thickness of its ends (h=1). Hence, there are two abrupt jumps in elastic thickness. We model this in 1D with the half-station method using 10 001 nodes over the model length with changes in thickness instantaneous over 1 node. We simulate the same problem with a finite-element solution of the equations of elastic equilibrium in a 2D elastic beam, assuming an isotropic, elastic material, using the FENICS finite-element code, and adapting the 2D linear elasticity tutorial example for our case (Bleyer, 2018). The change in beam elastic thickness is simulated in this case by a geometric reduction in finite beam thickness across its central section. Both beams are self-loaded by their own weight and associated body forces.

The results show a ∼3.5 % maximum difference between results. It should be noted that the discretisation of the finite-element solution is not exactly the same as the finite-difference one due to the complexities of gridding. The strain in this example is also several orders of magnitude larger than lithospheric strains (see the bending value u^′′).