the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A corrected finitedifference scheme for the flexure equation with abrupt changes in coefficient
Olivier Besson
The fourthorder 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 endloaded, 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 finitedifference discretisation of fourthorder, 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 onedimensional form and also give the twodimensional discretisation based on the new scheme. We also describe a general issue concerning the numerical stability of any secondorder finitedifference discretisation of a fourthorder differential equation like that describing flexure wherein contrasting magnitudes of coefficients of different summed terms lead to roundoff problems, which in turn destroy matrix positivity. We explain the use of 128 bit floatingpoint storage for variables to mitigate this issue.
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 preplate 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 manydecadeslong 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).
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 postplate 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érezGussinyé 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 socalled 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 fourthorder differential equation in u given as
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 platewide 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
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.
We note that the load term (ρ_{M}−ρ_{F})gu corresponds to two body forces: ρ_{M}gu is the restoring force (per unit length and width with units of Newtons) due to displacement of mantle by deflection of the plate, and ρ_{F}gu 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 ρ_{L}gq(x) defines a separate load of potentially different density to that assumed for the infill.
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).
With the load term now entirely on the righthand side of the equation, we see it consists of two parts with the infill part, ρ_{F}gu (right hand side), dependent on the deflection u(x) for which we are solving. However, a fixed load ρ_{L}gq(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 ρ_{F}gq(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
An alternative development of the equation for solution is by division through by (ρ_{M}−ρ_{F})g, giving
Here we see immediately that in contrast to Eq. (4), a direct, noniterative 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
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 flexuredependent 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.
For the most part, finitedifference methods have been applied to solve the flexure equation numerically in both the onedimensional beamtype situation and for twodimensional, thin elastic sheets. For the simplest case of constant flexural rigidity, D, the lefthand side of Eq. (3), for instance, is discretised as (see Fig. 4)
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 twodimensional 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, onedimensional 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 onedimensional 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
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 twodimensional sheetlike problem as well (see passing from Eq. 3 to 7 of Van Wees and Cloetingh, 1994). The resulting finitedifference 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 “halfstation” method (Cyrus and Fulton, 1966, 1968) and avoids a product rule derivation prior to discretisation entirely. Instead, Eq. (2) is directly transformed into a finitedifference approximation, replacing the derivatives both within and outside the brackets with finitedifference 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 fourthorder 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 halfstation 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 halfstation and wholestation methods are proven (Cyrus and Fulton, 1966, 1968) to be secondorder accurate, finitedifference 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 wholestation 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 halfstation method is convergent (see Appendix C).
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 (socalled “wholestation” 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 wholestation and halfstation 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 wholestation 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 endloaded 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 halfstation method gives exactly the result expected, showing that it corresponds to an endloaded 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 wholestation method gives no meaningful result and is unable to simulate a “broken” plate.
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 rectangularshaped 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 halfstation 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 “raftlike” response. Again, for such a case, the wholestation 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 basementtransecting faults. McQueen and Beaumont (1989) initially created a simple force balance model wherein horizontal stress across a faultbounded 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 blocktilting 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 “selfloading” 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, selfloading 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^{12} 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.
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 wholestation 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 halfstation 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 wholestation scheme will ultimately reach a point at which it returns no solution at all, whilst the halfstation method is stable for any combination of loads and variations of h(x).
A more general mathematical issue concerns the positivity of any numerical solution to a fourthorder differential equation of this kind. As can be verified (see Appendix B), all discretisations of the flexure equation using secondorder finitedifference 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 ρ_{M}g, 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, $\frac{\mathrm{6}D}{\mathrm{d}{x}^{\mathrm{4}}}$ and ρ_{M}g, becomes of critical importance. This is due to the issue of roundoff, 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^{3}, 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 highestresolution, 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 TANDEMX data (DLR – German Aerospace Center, 2018) are relatively freely accessible to scientists on a 0.4 arcsec (∼10–12 m) grid. Hence, $\frac{\mathrm{6}D}{\mathrm{d}{x}^{\mathrm{4}}}$ will quite conceivably be of magnitude 10^{19} or more in “real” problems in geodynamics. The small term ρ_{M}g will always be of order 10^{4}. Under such conditions, roundoff will lead to an “effectively singular” matrix and the numerical problem will fail. One workaround is relatively easily available. The use of quadrupleprecision (128 bit) representation of floatingpoint 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 “realworld” dimensions and parameters will be correctly dealt with by the numerical algorithm.
Besides the general issue described above, the additional term for the platewide 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.
The halfstation 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 wholestation method (applying the product rule first), by contrast, is unable to be relied upon to do so. Perhaps the most worrying aspect of the wholestation 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 wholestation 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 wholestation discretisation would perform safely. However, it is also known that the halfstation method works with smoothly varying functions of D, and in general, the errors associated with the halfstation method are always smaller than those of the wholestation 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 finitedifference discretisation of the flexure equation should be carried out using the halfstation method.
The wider application of the halfstation 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 onedimensional 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 finitedifference 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 highresolution 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 wholestation (product rule) derivation of the numerical scheme and should be revised.
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 fourthorder 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 floatingpoint 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 onedimensional problems in this paper, it is nevertheless clear that everything shown here extends to twodimensional, thin elastic plate formulations as well. To this end, we give the twodimensional halfstation formulation and discretisation of the problem (Appendix B). We will discuss twodimensional solutions using this scheme in forthcoming papers.
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.
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 platewide horizontal stress. The term ρ_{M}⋅g (left hand side) represents a restoring force due to displaced mantle. On the righthand side of the equation, ρ_{L}gq(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). ρ_{F}gu 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
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.
B1 Halfstation discretisation
We apply secondorder finitedifference operators simultaneously for both second derivatives inside and outside the brackets in Eq. (A1), which is something referred to as the halfstation method (Cyrus and Fulton, 1968).
Hence, if
where dx is the grid spacing, and i is the node number, then
We discretise the whole term in brackets first on a grid $i=\mathrm{1},\mathrm{\dots},N$:
Then, substituting the terms in brackets and advancing the indices gives us
Collecting terms, we obtain
Discretising the remaining parts of the equation then gives
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
which is a nonlinear series of equations in u. We reformulate this as a recursive matrix fixedpoint 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 Wholestation discretisation
The wholestation discretisation begins from the result of applying the product rule to Eq. (B1), giving us
Discretisation involves applying secondorder finitedifference schemes directly to all derivatives.
Hence,
Adding the remaining terms from Eq. (B8) and explicitly writing to show the relationship to Eq. (B9), we have the following.
The same procedure as for the halfstation discretisation is employed to solve these equations. The iteration could be made more efficient.
B3 Halfstation, twodimensional discretisation
As has been discussed, all twodimensional, 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 (wholestation). It is equally possible to apply a halfstation derivation, however. We start from Eq. (3) of Van Wees and Cloetingh (1994), and assuming that ν is constant, we can write
where ${u}_{xx}={\displaystyle \frac{{\partial}^{\mathrm{2}}u(x,y)}{\partial {x}^{\mathrm{2}}}}$ and so on. Defining partial finitedifference operators on a twodimensional grid (i,j),
and assuming $\mathit{\lambda}:=\mathrm{d}x=\mathrm{d}y$, Eq. (B11) can be discretised term by term in a way similar to the halfstation method applied to the onedimensional form. We have used maxima to derive the solution.
Here we give a general proof of the limitations of the wholestation discretisation as well as a comparison of a single case of the application of a 1D halfstation problem to an (approximately) equivalent, 2D, finiteelement, numerical solution.
C1 Proof of convergence
Consider the simplified beam equation:
for $x\in \mathrm{\Omega}=(\mathrm{0},L)$, and with the boundary conditions
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}>0 such that D(x)≥D_{0} a.e.,

α∈L^{2}(Ω), 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 finitedifference operators.
When numerically solving Eq. (C1), two different methods can be used.

One is the halfstation 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 wholestation 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^{2} when the function D is twice differentiable.
C2 Wholestation approach and regularity of D
Let us show that this formula makes sense only if D is at least twice differentiable. The wholestation discretisation can be written as
Now consider the academic case L=1 with two different simple functions D:

$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}$;

$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 $\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 $\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}$.
Conclusion. In both cases it is impossible to verify the numerical accuracy of the wholestation discretisation. Moreover, this scheme leads to incorrect numerical solutions when D is not sufficiently regular.
Remark. In the halfstation 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 wholestation and the halfstation methods give the same results.
C3 2D numerical comparison
In Fig. C1 we show a case in which only the halfstation 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 halfstation method using 10 001 nodes over the model length with changes in thickness instantaneous over 1 node. We simulate the same problem with a finiteelement solution of the equations of elastic equilibrium in a 2D elastic beam, assuming an isotropic, elastic material, using the FENICS finiteelement 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 selfloaded 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 finiteelement solution is not exactly the same as the finitedifference 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^{′′}).
The codes used in the preparation of this paper are available from the GitHub repository at https://github.com/davidhindle/flexure1dhs (Hindle, 2021a) and on Zenodo (https://doi.org/10.5281/zenodo.4643989; Hindle, 2021b).
No data sets were used in this article.
DH had the idea to use the halfstation method of discretisation (at the time without realising he was doing so and by complete accident). OB pointed out the implications of the discretisation, discovered the earlier work on it, and also explained the reasons behind the instability of the numerical method when roundoff becomes an issue. OB modified the code to work with 128 bit, quadrupleprecision, floatingpoint representation. OB developed a number of benchmark codes to check results against original papers. OB wrote the maxima script to get the correct halfstation derivation of the twodimensional problem. The main codes used in the paper were written by DH. The main text, figures, and experiments were all prepared by DH.
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
David Hindle thanks Hans Petter Langtangen posthumously for his enormous contributions to the process of learning and understanding mathematics for natural scientists. David Hindle only discovered Hans Petter's recognition of the halfstation method amongst the many, rich web archives of course material Hans Petter left to us all. This is a testament to the enduring influence of an evidently much loved, talented, and inspiring individual, as well as his wonderful scientific career, which was cut far too short. Both authors thank the editor, Tara Gerya, and the reviewers, Chris Beaumont, Stefan Schmalholz, and one anonymous referee, for their helpful comments and patience with the subject matter.
This openaccess publication was funded by the University of Göttingen.
This paper was edited by Taras Gerya and reviewed by Stefan Markus Schmalholz, Christopher Beaumont, and one anonymous referee.
Adams, R. A. and Fournier, J. J.: Sobolev spaces, Academic Press, ISBN 0120441438, 2003. a, b
Audet, P. and Mareschal, J.C.: Variations in elastic thickness in the Canadian Shield, Earth Planet. Sci. Lett., 226, 17–31, https://doi.org/10.1016/j.epsl.2004.07.035, 2004. a
Axelsson, O.: Iterative Solution Methods, Cambridge University Press, https://doi.org/10.1017/CBO9780511624100, 1994. a
Barrell, J.: The strength of the Earth's crust, J. Geol., 22, 655–683, https://doi.org/10.1086/622145, 1914. a, b
Beaumont, C.: Foreland basins, Geophys. J. Int., 65, 291–329, https://doi.org/10.1111/j.1365246X.1981.tb02715.x, 1981. a, b
Bleyer, J.: Numerical Tours of Computational Mechanics with FEniCS, Zenodo [code], https://doi.org/10.5281/zenodo.1287832, 2018. a
Bodine, J., Steckler, M., and Watts, A.: Observations of flexure and the rheology of the oceanic lithosphere, J. Geophys. Res.Sol. Ea., 86, 3695–3707, https://doi.org/10.1029/JB086iB05p03695, 1981. a
Braun, J., Deschamps, F., Rouby, D., and Dauteuil, O.: Flexure of the lithosphere and the geodynamical evolution of noncylindrical rifted passive margins: Results from a numerical model incorporating variable elastic thickness, surface processes and 3D thermal subsidence, Tectonophysics, 604, 72–82, https://doi.org/10.1016/j.tecto.2012.09.033, 2013. a
Burkhard, M. and Sommaruga, A.: Evolution of the western Swiss Molasse basin: structural relations with the Alps and the Jura belt, Geological Society, London, Special Publications, 134, 279–298, https://doi.org/10.1144/GSL.SP.1998.134.01.13, 1998. a
Burov, E. and Watts, A.: The longterm strength of continental lithosphere: “jelly sandwich” or “crème brûlée”?, GSA today, 16, 4, https://doi.org/10.1130/10525173(2006)016<4:TLTSOC>2.0.CO;2, 2006. a
Cyrus, N. J. and Fulton, R. E.: An accuracy study of finite difference methods in structural analysis, NASA technical reports, technical memorandum, Fourth Conference on Electronic Computation, American Society of Civil Engineers, ntrs.nasa.gov:19680018134, 1966. a, b, c, d, e
Cyrus, N. J. and Fulton, R. E.: Accuracy of finite difference methods, NASA Technical Note TN D4372, 1968. a, b, c, d, e, f
DeCelles, P. G. and Giles, K. A.: Foreland basin systems, Basin Res., 8, 105–123, https://doi.org/10.1046/j.13652117.1996.01491.x, 1996. a
DLR – German Aerospace Center: TanDEMX – Digital Elevation Model (DEM) – Global, 12 m, https://gdk.gdide.org/geonetwork/srv/api/records/5eecdf4cde57462499e960086b032aea, (last access: 4 July 2021), 2018. a
Earth Resources Observation And Science (EROS) Center: Shuttle Radar Topography Mission (SRTM) 1 ArcSecond Global, https://doi.org/10.5066/F7PR7TFT, 2017. a
Egan, S.: The flexural isostatic response of the lithosphere to extensional tectonics, Tectonophysics, 202, 291–308, https://doi.org/10.1016/00401951(92)90115M, 1992. a
Fisher, O.: Results of a Transcontinental series of gravity measurements, Nature, 52, 433–435, https://doi.org/10.1038/052433a0, 1895. a
Gilbert, G.: The strength of the earth's crust, B. Geol. Soc. Am., 1, 23–27, 1889. a
Gilbert, G. K.: Lake Bonneville, vol. 1, United States Geological Survey, https://doi.org/10.3133/m1, 1890. a
Govers, R., Meijer, P., and Krijgsman, W.: Regional isostatic response to Messinian Salinity Crisis events, Tectonophysics, 463, 109–129, https://doi.org/10.1016/j.tecto.2008.09.026, 2009. a
Gunn, R.: A quantitative study of isobaric equilibrium and gravity anomalies in the Hawaiian Islands, J. Frankl. Inst., 236, 373–390, https://doi.org/10.1016/s00160032(43)902753, 1943a. a
Gunn, R.: A quantitative evaluation of the influence of the lithosphere on the anomalies of gravity, J. Frankl. Inst., 236, 47–66, https://doi.org/10.1016/S00160032(43)911986, 1943b. a, b, c, d, e, f, g
Gunn, R.: Quantitative aspects of juxtaposed ocean deeps, mountain chains and volcanic ranges, Geophysics, 12, 238–255, https://doi.org/10.1190/1.1437321, 1947. a
Hagen, E. S., Shuster, M. W., and Furlong, K. P.: Tectonic loading and subsidence of intermontane basins: Wyoming foreland province, Geology, 13, 585–588, https://doi.org/10.1130/00917613(1985)13<585:TLASOI>2.0.CO;2, 1985. a, b
Hindle, D.: davidhindle / flexure1dhs, GitHub repository [code], https://github.com/davidhindle/flexure1dhs (last access: 30 March 2021), 2021. a
Hindle, D.: davidhindle/flexure1dhs: flexure1dhs_v1.0beta (v1.1), Zenodo [code], https://doi.org/10.5281/zenodo.4643989, 2021. a
Hindle, D. and Kley, J.: The Subhercynian Basin: an example of an intraplate foreland basin due to a broken plate, Solid Earth, 12, 2425–2438, https://doi.org/10.5194/se1224252021, 2021. a
Horton, B. K. and DeCelles, P. G.: The modern foreland basin system adjacent to the Central Andes, Geology, 25, 895–898, https://doi.org/10.1130/00917613(1997)025<0895:TMFBSA>2.3.CO;2, 1997. a
Kaban, M., Chen, B., Tesauro, M., Petrunin, A., El Khrepy, S., and AlArifi, N.: Reconsidering effective elastic thickness estimates by incorporating the effect of sediments: A case study for Europe, Geophys. Res. Lett., 45, 9523–9532, https://doi.org/10.1029/2018GL079732, 2018. a, b
Karner, G. and Watts, A.: Gravity anomalies and flexure of the lithosphere at mountain ranges, J. Geophys. Res.Sol. Ea., 88, 10449–10477, https://doi.org/10.1029/JB088iB12p10449, 1983. a
Langtangen, H. P.: Finite difference methods for wave motion, http://hplgit.github.io/nummethodsforPDEs/doc/pub/wave/pdf/wave4print.pdf, last access: 3 November 2016. a
McKenzie, D. and Fairhead, D.: Estimates of the effective elastic thickness of the continental lithosphere from Bouguer and free air gravity anomalies, J. Geophys. Res.Sol. Ea., 102, 27523–27552, https://doi.org/10.1029/97JB02481, 1997. a
McQueen, H. and Beaumont, C.: Mechanical models of tilted block basins, Origin and Evolution of Sedimentary Basins and Their Mineral Resources, edited by: Price, R. A., AGU/IUGG Monograph, 48, 65–71, 1989. a, b, c, d
NASA/METI/AIST/Japan Spacesystems And U.S./Japan ASTER Science Team: ASTER Global Digital Elevation Model, https://doi.org/10.5067/ASTER/ASTGTM.002, 2009. a
PérezGussinyé, M. and Watts, A.: The longterm strength of Europe and its implications for plateforming processes, Nature, 436, 381–384, https://doi.org/10.1038/nature03854, 2005. a
Sebben, S. and Baliga, B. R.: Some extensions of tridiagonal and pentadiagonal matrix algorithms, Numerical Heat Transfer, Part B Fundamentals, 28, 323–351, https://doi.org/10.1080/10407799508928837, 1995. a
Stewart, J. and Watts, A.: Gravity anomalies and spatial variations of flexural rigidity at mountain ranges, J. Geophys. Res.Sol. Ea., 102, 5327–5352, https://doi.org/10.1029/96JB03664, 1997. a, b, c, d, e
Thor, R.: Mapping the thickness of the Martian elastic lithosphere using maximum likelihood estimation, http://resolver.tudelft.nl/uuid:ca2a860fa78e4517a829eb03a7254f45, last access: 28 June 2016. a
Van der Meulen, M., Buiter, S., Meulenkamp, J., and Wortel, M.: An early Pliocene uplift of the central Apenninic foredeep and its geodynamic significance, Tectonics, 19, 300–313, https://doi.org/10.1029/1999TC900064, 2000. a, b
Van Wees, J. and Cloetingh, S.: A finitedifference technique to incorporate spatial variations in rigidity and planar faults into 3D models for lithospheric flexure, Geophys. J. Int., 117, 179–195, https://doi.org/10.1111/j.1365246X.1994.tb03311.x, 1994. a, b, c, d, e, f
Walcott, R.: Flexural rigidity, thickness, and viscosity of the lithosphere, J. Geophys. Res., 75, 3941–3954, https://doi.org/10.1029/JB075i020p03941, 1970a. a, b, c
Walcott, R.: Flexure of the lithosphere at Hawaii, Tectonophysics, 9, 435–446, https://doi.org/10.1016/00401951(70)900569, 1970b. a
Walcott, R.: Isostatic response to loading of the crust in Canada, Can. J. Earth Sci., 7, 716–727, https://doi.org/10.1139/e70070, 1970c. a
Walcott, R.: Late Quaternary vertical movements in eastern North America: Quantitative evidence of glacioisostatic rebound, Rev. Geophys., 10, 849–884, https://doi.org/10.1029/RG010i004p00849, 1972. a
Watts, A.: The effective elastic thickness of the lithosphere and the evolution of foreland basins, Basin Res., 4, 169–178, https://doi.org/10.1111/j.13652117.1992.tb00043.x, 1992. a, b
Watts, A. B.: Isostasy and Flexure of the Lithosphere, Cambridge University Press, 2001. a, b
Wessel, P., Luis, J., Uieda, L., Scharroo, R., Wobbe, F., Smith, W., and Tian, D.: The generic mapping tools version 6, Geochem. Geophys. Geosyst., 20, 5556–5564, https://doi.org/10.1029/2019GC008515, 2019. a, b, c, d
Wickert, A. D.: Opensource modular solutions for flexural isostasy: gFlex v1.0, Geosci. Model Dev., 9, 997–1017, https://doi.org/10.5194/gmd99972016, 2016. a, b, c
 Abstract
 Introduction
 Numerical flexure solutions
 Comparison of the different numerical schemes
 Other numerical issues
 Implications
 Conclusions
 Appendix A: General aspects of the equation when solved numerically
 Appendix B: Discretisation schemes
 Appendix C: Convergence, numerical comparisons
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Numerical flexure solutions
 Comparison of the different numerical schemes
 Other numerical issues
 Implications
 Conclusions
 Appendix A: General aspects of the equation when solved numerically
 Appendix B: Discretisation schemes
 Appendix C: Convergence, numerical comparisons
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References