Estimating the lithospheric flexure of a plate with non-uniform flexural rigidity: a quantitative modelling approach

Lithospheric deformation is a fundamental process in plate tectonics. It is, therefore, critical to determine how the lithosphere responds to geological loads to better understand tectonic processes. The lithosphere can be modelled as the flexure of a thin, elastic plate over long-term (> 105 yr) geological timescales. The partial differential equation for the flexure of an orthotropic plate is used indirectly to calculate theoretical admittance and coherence, which are then compared against the observed admittance and coherence to invert for the non-uniform flexural rigidity (or effective elastic thickness, Te) of the plate. However, the process for accurately recovering variable lithospheric flexure remains unresolved, as the classical lithospheric model may overestimate the deflection of the plate. Here we adopt the classic lithospheric model with applied external and internal loads at the surface and Moho, respectively, and assume that the compensation material is denser than the mantle material beneath the Moho. The lithospheric flexure errors are derived mainly from the Te and Moho recovery errors in this lithospheric model. Synthetic modelling is then performed to analyse the influence of the Te and Moho errors. The analysis of synthetic modelling shows that: (1) the Te error-induced flexure errors exhibit a rippling pattern, and the rippling pattern is broader in high Te regions; (2) the Moho error-induced flexure errors mainly occur in the low Te regions, and applying Airy isostasy theory in low Te regions may still greatly overestimate the lithospheric deformation amplitude; and (3) the lithospheric flexure errors are dominated by the Te and Moho errors in the high and low Te regions, respectively.


Introduction
Our understanding of the mechanics of the Earth's crust and lithosphere have been greatly advanced since the 1970s, via numerous modelling, including experimental examinations of the behaviour of physical bodies (e.g., Lewis 1970, 1972;Lewis and Dorman 1970a, b) and derivations of the mathematics that describe lithospheric processes (e.g., McNutt 1980;Forsyth 1985;Burov and Diament 1992;Banks et al. 2001;Braitenberg et al. 2002;McKenzie 2003;Kirby and Swain 2009;Zhang et al. 2018a, b, c;Zhang et al. 2019a). These studies followed the principle that the lithosphere can be idealised if the main structural features are considered, thereby allowing simplified models to capture the key structures and geodynamics of the lithosphere. However, no information on the internal structure of the lithosphere has been used in these investigations, and its mechanical properties, such as the flexure of the tectonic plates, may be incorrectly estimated (Ribe 1982).
The lithospheric strength of tectonic plates reflects their resistance to vertical deformation in response to geological loads over long-term (> 10 5 yr) geological timescales (Watts and Burov 2003). This assumption allows the lithosphere to be modelled as the flexure of Open Access Full list of author information is available at the end of the article a thin, elastic plate via the partial differential equation for the flexure of an orthotropic plate (Timoshenko and Woinowsky-Krieger 1959). This equation is usually used indirectly to calculate theoretical admittance and coherence, which are then compared against observed ones to invert for the non-uniform flexural rigidity (or effective elastic thickness, T e ) of the plate (Kirby 2014). Numerous studies have used T e to analyse lithospheric rheology and deformation (Pérez-Gussinyé et al. 2009;Chen et al. 2015;Ji et al. 2017Ji et al. , 2020Lu et al. 2020). However, this idealised term does not refer to an existing thickness or physical layer within the Earth, but instead corresponds to the thickness of an ideal elastic plate that undergoes the same deformation as the lithosphere under the same loads (Watts 2001). As T e is a periphrastical parameter for understanding lithospheric rheology and deformation, directly modelling lithospheric flexure should provide key insights into the nature of tectonic evolution and dynamics.
Here we use the classic lithospheric model with applied external and internal loads at the surface and Moho, respectively, and assume that the compensation material is denser than the mantle material beneath the Moho (Fig. 1). The density contrast between the compensation material and mantle material beneath the Moho is set to 200 kg/m 3 , based on the density contrast between the uppermost and basal lithospheric mantle (Kaban et al. 2016). Notably, in some special areas, loads may be taken from specific geological attributes instead of topographies. In the subduction zones, for instance, slab pull may be considered as the main loads and usually be handled as boundary loading (e.g., Zhang et al. 2014Zhang et al. , 2021, applying the corresponding plate flexure models has made a great deal achievements (e.g., Zhang et al. 2018a, b, c;2019a, b). Taking the topographies, however, as the loads has a widespread application in estimating plate flexure and should have a wider prospect. The model used in this study is still an idealised model that considers the vertical density variation in the lithosphere. Furthermore, we obtain T e variations via the fan wavelet coherence method Swain 2004, 2008) and the Moho topography via a constrained gravity inversion method (Wu et al. 2017). Based on the finite-difference (FD) method shown by Kirby and Swain (2008), we reformat the equations to invert the lithospheric flexure of a two-layer lithospheric model. The errors associated with the lithospheric flexure estimation in such a lithospheric model are mainly derived from the T e and Moho errors. Synthetic modelling is then performed to analyse the influence of the T e and Moho errors on the lithospheric flexure estimation.

Methodology
The key to calculating lithospheric flexure via the FD method is to obtain accurate spatial distributions of the Moho and T e . The Moho depth is usually recovered via gravity inversion (Wu et al. 2017). There are currently two main methods for mapping the spatial variations in T e based on the coherence between the Bouguer gravity anomaly (BGA) and topography/bathymetry. One is the multitaper method (McKenzie and Fairhead 1997;Pérez-Gussinyé et al. 2004), and the other is fan wavelet coherence method Swain 2004, 2008). The characteristic flexural wavelength is correlated with the strength of the lithosphere, and the flexural amplitude varies with both T e and the magnitude of the load (Watts 2001 Fig. 1 Conceptual mechanical model of the lithosphere and asthenosphere, the compensation depth is set to be at the Lithosphereasthenosphere interface (LAI). ρ c , ρ m and ρ F are density values (Table 1), and �ρ 1 = ρ c − ρ f , �ρ 2 = ρ m − ρ c , �ρ 3 = ρ F − ρ m are interface density contrasts, ρ f represents either the air density on land or seawater density in the ocean (revised from Banks et al. 2001) on a stronger lithosphere can be truncated when the selected window size is too small, resulting in an artificially low T e . Conversely, T e can be blended across multiple features, with a potential bias toward the largest feature in the window, when the selected window size is too large (Pérez-Gussinyé et al. 2004Kalnins and Watts 2009). The fan wavelet coherence method can alleviate these issues by allowing the T e features to be located at different wavelengths, and is employed to estimate T e in this study.

Gravity inversion for Moho topography
Joint inversions that involve gravity measurements and other geophysical constraints are generally employed to reduce the non-uniqueness of the best-fit models (Chappell and Kusznir 2008;Bai et al. 2014;Ji et al. 2018). Here we perform a three-dimensional gravity inversion that is constrained by seismic observations (Wu et al. 2017) to recover an accurate Moho topography. The Moho topography that corresponds to the complete Bouguer gravity anomaly can be expressed as h(r) = h e (r) + �h(r) , where h e (r) is the initial interface. The difference between the Moho and initial interface can then be obtained by substituting this difference into the gravity inversion equation (Oldenburg 1974): where F [] is the Fourier transform, g is the complete Bouguer gravity anomaly, ς is Newton's gravitational constant, ρ is the density contrast across the interface, k is the wave vector of the transformed function, r is the entire x-y plane and z 0 is the reference depth of the density interface, n is the iterations, whose value depends on the convergence rate of Eq. (1) and the error threshold.
We largely follow the procedure outlined in Wu et al. (2017), except for the implementation of Eq. (1) to determine the interface difference during each iteration of the inversion, to ensure that both the Moho depth and Bouguer gravity anomaly satisfy the inversion equation (Oldenburg 1974) after the Moho topography is recovered. The detailed process is as follows: (1) linearly fit the Moho depth values and the corresponding Bouguer gravity anomalies with the constraint points; (2) then, use the fitting parameters to form the initial Moho grid h e (r) with the observed Bouguer gravity anomaly grid; (3) perform the iterative algorithm of Eq. (1); (4) the final Moho topography is obtained iteratively until the interface difference between adjacent iteration steps meets the error threshold, which is set to be 0.001 km in this study.

Thin plate flexure
For quick reference, we provide the key equations for estimating thin plate flexure, which can be easily calculated via the approach of Kirby and Swain (2008). We consider an elastic slab model that involves only the initial external (at the surface) and internal (at the Moho) loads, which are represented as ρ c − ρ f gh i and [(ρ m − ρ c )gw i ] , respectively. As the plate model is a linear system, the deflection can be divided into four components. The deflection of amplitude w T at the Moho and h T at the surface caused by the initial external load, w B at the Moho and h B at the surface caused by the initial internal load. Then, the final Moho topography after flexure is and the final surface topography after flexure is The final deflection amplitude v can then be expressed as The relations of initial loads and deflection amplitude in the partial differential equation for flexure of an orthotropic plate are given by Timoshenko and Woinowsky-Krieger (1959): where the bending and twisting moments are and where D is the flexural rigidity, σ is Poisson's ratio. Furthermore, substituting the initial loads with Eq. (4), and simplifying Eq. (5) with the Laplace operator that gives (Timoshenko & Woinowsky-Krieger 1959;Stark et al. 2003): where all variables are functions of the (x, y) position, � = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the Laplace operator. The minus sign on the right side of Eq. (6) demonstrates the plate acts as resistance to the topographic deformation. When D reaches zero, Eq. (6) follows the Airy isostatic theory. The solution to Eq. (6) via the FD method is shown in Appendix A.

A larger Winkler foundation
A foundation that acts with a force proportional to the deflection at every point is known as a Winkler foundation (e.g., Huang and Thambiratnam 2001;Mofid and Noroozi 2009). An elastic foundation is usually expressed as ρ m gv in linear elastic plate theory. However, the deflection predicted via linear theory is always too high; this non-linearity decreases with increasing plate thickness, such that plates with T e > 20 km are only slightly nonlinear (Ribe 1982). This overestimated deflection should be suppressed in ocean settings since T e is generally low in these settings. Here we use a two-layer model ( Fig. 1), whereby the internal load is confined to the Moho, such that the entire lithosphere is modelled with a lithospheric elastic foundation with higher density ( ρ F ) than that of the uppermost layer of the mantle ( ρ m ). Equation (6) is then changed to The main advantages of this model are: 1. The overestimated deflection of the lithosphere predicted via linear theory can be suppressed; 2. The density of the mantle usually increases with depth, so the model fits real data better in the framework of lithospheric-scale deflection; 3. The main inversion function (Eq. (6)) in such a model can be converted to Eq. (7), which can be considered as a forward equation with a small Winkler foundation, thereby ensuring that stable calculation results are always obtainable.

Fan wavelet coherency method
The wavelet method for calculating the coherence between the topography and gravity data was first proposed by Stark et al. (2003), who employed a Gaussian window function as the transform kernel to estimate T e in southern Africa. The continuous wavelet transform offers the ability to obtain the local phase and amplitude of the topography and gravity anomalies, whereas the Fourier transform requires wavenumber averaging of the spectral estimates. The Gaussian wavelet fails to reproduce the Fourier power spectrum exactly (Kirby 2005). Two-dimensional Morlet wavelets, which are Gaussian-modulated complex exponential functions, can reproduce the Fourier spectrum (Kirby and Swain 2004). The spatial resolution of T e can be improved by adjusting the central wavenumber of the Morlet wavelet (Kirby and Swain 2011). Only the square of the real component of the wavelet coherence needs to be employed to estimate T e Swain 2008, 2009). The coherency, which is formed by obtaining auto-and cross-spectra at a series of azimuths over the 0°-180° range in a 'fan' and then averaging them (Kirby and Swain 2009), is computed as the real and positive coherence, which suffers from the problem of information loss. The squared real coherency can be expressed as a function of scale s and position x: where B sxθ and H sxθ represent the complex wavelet coefficients of the Bouguer gravity anomaly and topography, respectively, the * indicates complex conjugation, and the angular brackets denote azimuthal averaging performed over the entire 180° azimuthal range. The continuous wavelet transform (CWT) is performed through the following equations: Fourier transform of the 2D Morlet wavelet, scale s is also the dilation parameter, the "central wavenumber" k 0 = π 2/ ln 2 ≈ 5.336 is used in this study. Note that v only represents the frequency domain coordinate parameter in Eq. (9), but elsewhere in this article is the flexure.
Flexural components in Eqs. (2) and (3) can be replaced with equations including initial internal and external loads w i and h i . As these two equations are linear, they form the so-called "load-deconvolution equations" in the wavenumber domain with the assumption that T e is uniform (Forsyth 1985;Kirby and Swain 2008). Following the formulation in Kirby and Swain (2011) and making the equations available for a two-layer model that The parameters are where G, H mean the observed Bouguer gravity anomaly and surface topography, | k| = k 0 /s is the equivalent Fourier wavenumber, Z m is the Moho depth, Z F is the LAI depth, ς , E and σ are constant values (Table 1).
After Eq. (10) is solved, initial loads and flexural components can be obtained. Then, the predicted wavelet coherency is where T is for surface load components, and B for subsurface load components. With the advantage of averaging being performed over the wavelet azimuth in the wavelet method, Eq. (12)  mean the azimuthal averaged wavelet coefficients of initial surface and subsurface loads. By considering each local wavelet spectrum to be independent of adjacent spectra (the "decoupling" assumption), the wavelet coherency can be applied to approximately invert for local T e (Stark et al. 2003;Kirby and Swain 2004). Then the final T e distributions are obtained by finding the T e that provides the minimum misfit between the observed and predicted coherencies at each grid point, using Brent's method of 1-D minimization (e.g., Press et al. 1992).

Synthetic modelling
Synthetic modelling is commonly used to test the accuracy of the wavelet method (Stark et al. 2003;Swain 2004, 2008). Random fractal initial loads and T e models are used as inputs to the flexural equation, which outputs the final surface topography (h) and Moho topography (w). We follow the same approach in our synthetic modelling. We employ Eq. (6) in the inversion since it is a thorough rearrangement of the flexural equation (Eq. (11)) in Kirby and Swain (2009) where the inputs are the final topography (h) and a random fractal initial internal load (w i ).
The equivalent topography, T e and Moho topography, are required to estimate lithospheric flexure. The equivalent topography can be specifically calculated using a global elevation model, sediment thickness model and sediment thickness-density conversion formula (12) 1/2 , (e.g., Sclater & Christie 1980;Ji et al. 2020), such that this random fractal initial load (h i ) and flexure (v) derived topography (h) can be used in the synthetic modelling. The influences of both the recovered T e errors and inverted Moho errors on the FD-derived lithospheric flexure are investigated to determine how these inputs shape the model results.

Model generation
We follow the synthetic modelling procedure of Macario et al. (1995), and use the mid-point displacement method (Peitgen and Saupe 1988) to generate two-dimensional fractal Brownian surfaces that represent the initial external (surface, h i ) and internal (subsurface, w i ) loads. In addition, the T e model is also obtained by applying the mid-point displacement method, as Kirby and Swain (2008) noted: "fractal T e distributions are likely to be closer to the real Earth's T e structure." A fractal dimension of 2.5 is used for the surface (Fig. 2a) and subsurface loads (Fig. 2b), as well as the T e model (Fig. 2c). The grids consist of 255 × 255 points at a 20-km spacing and are scaled to possess a realistic amplitude. Furthermore, the T e grid is low-pass filtered with a cutoff wavelength of 150 km. We choose a central wavenumber of 5.336 in the subsequent fan wavelet coherence inversion since we place greater importance on the T e accuracy than its spatial resolution. We do not set a constant F (loading ratio), as in Kirby and Swain (2008). Although this approach may go against Forsyth's fourth assumption about the proportionality of the power spectral densities of the initial loads (Forsyth 1985;Simons and Olhede 2013), the feasibility of the T e inversion with multiple varying F interfaces via the fan wavelet coherence method has already been tested by Kirby and Swain (2009). The topography (Fig. 2d), Moho (Fig. 2e) and lithospheric flexure models (Fig. 2f ) are obtained by applying the FD method to Eq. (7). Before applying the FD method, all the data sets were mirror extended to be 9 times of their original size to reduce the computation errors caused by boundary condition, and the extended boundaries were removed after all the calculation processes were completed. The Moho interface is then used to obtain the forward Bouguer gravity anomaly (Fig. 3), with the relevant parameters listed in Table 1. The Moho inherits high-frequency information from the initial subsurface. However, the modelled Bouguer gravity anomaly may suffer from signal degradation due to upward continuation of the signal since the observation surface is much higher than the Moho interface, resulting in attenuation of this highfrequency information.

Effect of T e errors
The fan wavelet coherence method is employed to invert for T e (Fig. 4a) using the generated surface topography, Moho topography and forward-modelled Bouguer anomaly grids, and to calculate its associated errors (Fig. 4b).
Note that, all the grids of errors are calculated by subtracting the forward-modelled grids from the inverted grids. This T e inversion recovers well the main T e features in the T e model. The errors are randomly distributed across the entire map, but the trend of the local errors is consistent with the T e gradient belts. This consistency is due mainly to the resolution of the transition wavelets, as well as two other adverse factors: (1) the decoupling assumption and (2) the assumption of uniform T e in the loading equations (Kirby and Swain 2008). We then used the generated topography, Moho topography and inverted T e grids to estimate (Eq. (7)) the lithospheric flexure (Fig. 5a) and the associated errors (Fig. 5b). The long-wavelength component of the two flexure maps is basically the same (Figs. 2f and 5a), and the amplitudes of the flexure errors (Fig. 5b) are up to one order of magnitude smaller than the T e errors (Fig. 4b).
The T e error-induced flexure errors exhibit a rippling pattern due to the deflection of a higher T e plate possessing a lower amplitude and slower amplitude decay rate (see Fig. 6), although such a rippling pattern is truncated and/ or deformed by the attitudes of T e and total load. Furthermore, the range of the rippling pattern is affected by the T e values and load amplitude, although it has a higher correlation with the T e values in this synthetic modelling study, as the rippling pattern is obviously broader in the high T e regions.
To further explore the relation between the flexure errors and the T e errors, we simulated another data sets with a new T e model, who possesses an opposite numerical distribution pattern against the first T e model, and is notated as T e r . More exactly, T e r = maximum(T e ) + minimum(T e )-T e (Fig. 7a). With the same procedure described above, we obtain the maps of lithospheric flexure with T e r (Fig. 7b), inverted T e r grid (Fig. 8a), T e r errors due to the inversion (Fig. 8b),  (Fig. 9a) and lithospheric flexure errors for the model inversion using the inverted T e r grid (Fig. 9b). The high T e region is concentrated in the central part of the T e r grid (Fig. 7a), this enhances the resistance of the central area to the applied loads and make the lithospheric flexure in the corresponding areas smoother (Fig. 8b). The inverted T e r grid still reflects well the long wavelength part of the T e r model, with the errors due to the inversion (Fig. 8b) much more concentrated at the central area. However, the consistency of the T e r errors and the T e r gradient belts is obvious as before. Figure 9 is dominated by the lithospheric flexure errors of rippling pattern, as significant errors are always surrounded by a series of errors in opposite sign. Though the difference between T e and T e r grids is great, the response of lithospheric flexure is quite similar.

Effect of Moho errors
Bouguer gravity anomaly reflects the long-wavelength Moho information. It follows that the short-wavelength component of the flexure-derived Moho will be the main source of error (Fig. 10b) when we use the Bouguer anomaly to invert for the Moho topography (Fig. 10a). Although the maximum Moho error is 6 km, the Moho inversion yields reliable results as the root-mean-square error is only 0.91 km. Furthermore, the even distribution of the Moho errors provides insights into how the Moho inversion errors impact the flexure errors.
We also used the generated topography, T e model and recovered (inverted) Moho topography to estimate (Eq. (7)) the lithospheric flexure (Fig. 11a) and the associated errors (Fig. 11b). A comparison of the Moho and flexure results yields the following findings. (1) The Moho inversion errors are distributed evenly across the entire region, whereas the main flexure errors occur in the low T e regions. (2) The Moho errors tend to induce flexure errors with similar magnitudes but opposite signs (see the white rectangles in Figs. 10b and 11b). (3) The plate still acts as a low-pass filter to the Moho disturbance, even in the low T e regions, with an increase in the intensity of filtering, where there is an increase in the T e values.
Although the flexure error in the low T e regions, e.g., oceanic lithosphere, is sensitive to the Moho disturbance, this method is still suitable for investigations of the oceanic lithosphere. The reference depth of the oceanic Moho is shallow, which may allow us to Easting (km) Fig. 6 Cross sections of the lithospheric flexure of uniform T e plates to a cuboid loading. The centre of the cuboid is placed at (2500, 2500) km, with height of 1 km, both length and width of 80 km, density of ρ c . The grey solid line is for a T e = 25 km plate, the grey dotted line is for a T e = 5 km plate and the black solid line is calculated by the grey solid curve minus the grey dotted curve retain more short-wavelength information, such that more accurate inversion results can be obtained. In addition, the FD method-derived flexure is suited for Airy isostasy theory in regions, where T e reaches zero. However, we have to mention that, directly apply Airy isostasy theory in some low T e regions may greatly overestimate the lithospheric deformation amplitude. For example, take the extreme point whose coordinate is (2.52, 2.54)*10 3 km in the lower white rectangle in Fig. 10b, the maximum Moho error reaches 4.258 km and we may expect a corresponding lithospheric flexure error of −(ρ m − ρ c )/(ρ F − ρ m ) * 4.258 = −9.5805 km with the Airy isostasy theory (the coefficient −(ρ m − ρ c )/(ρ F − ρ m ) comes from Eq. (7) by setting D to be zero). However, this extreme point in Fig. 10b gets a value of − 2.520 km, the amplitude is only 26.3% to the expected value from the Airy isostasy theory. This huge overestimation clarifies the significance of applying the FD method even in low T e regions when researching the vertical lithospheric/crustal deformation.

Integrated effects of the T e and Moho errors
Neither T e nor the Moho can be recovered with absolute accuracy using real data sets. We therefore investigated the integrated effects of the T e and Moho errors on the estimated lithospheric flexure errors to explicate the sources of the flexure errors. We used the generated topography, recovered T e and inverted Moho topography to invert (Eq. (7)) for the lithospheric flexure (Fig. 12a) and calculate their associated errors (Fig. 12b). The longwavelength component of the lithospheric flexure is still well recovered, indicating that the FD method is stable when an integrated disturbance is applied. The lithospheric flexure errors are dominated by the T e errors in the high T e regions (Figs. 5b, 11b and 12b). Correspondingly, the flexure errors are dominated by the Moho errors in the low T e regions, although the amplitudes of the flexure errors have been somehow changed in the low T e regions (Figs. 11b and 12b).

Conclusions
Here, we employed the classic lithospheric model with an applied external load at the surface and internal load at the Moho, and assumed that the compensation material was denser than the mantle material beneath the Moho, in an attempt to estimate lithospheric flexure. The resultant lithospheric model is closer to the actual lithospheric structure, with effective suppression of the overestimated lithospheric deflection predicted via linear theory, and conversion of the main lithospheric  (6)) into a forward equation with a small Winkler foundation (Eq. (7)). This lithospheric model, which based on T e estimates from the fan wavelet coherence method and gravity-derived Moho topography, yields a finite-difference solution of the partial differential equation that recalibrates the lithospheric flexure with high accuracy. Our main conclusions are as follows: 1. The T e error-induced flexure errors exhibit a rippling pattern, and the rippling pattern is broader in high T e regions.
2. The Moho error-induced flexure errors mainly occur in the low T e regions, and applying Airy isostasy theory in low T e regions may still greatly overestimate the lithospheric deformation amplitude. 3. The flexure errors in the low and high T e regions are dominated by the Moho and T e errors, respectively.

Appendix A: Finite-difference method
The FD solution in Eq. (5) is mainly from Chen (2013), but we use the final loads instead of the initial loads on the right side of the equation. Equation (5) can be further simplified to The differential operations in Eq. (A1) are approximated using central difference operators. A schematic diagram of the central difference grid (Ghali et al. 1989) is shown in Fig. 13, whereby a square grid is adopted, such that dx = dy . The corresponding central difference operators for the first-order, second-order and Laplace differential operations are shown in Fig. 14. The approximation to the differential operations is then undertaken with the central difference operators, followed by formula expansion and combination. The following coefficients of the indicated v i,j terms for an isotropic plate can then be obtained: where A = 1−σ 8 D i+1,j+1 − D i−1,j+1 − D i+1,j−1 + D i−1,j−1 , dx = dy = , and all of the coefficients must be divided by 4 . For Eq. (7), only the first term in Eq. (A2) should be changed to The data sets of the study area can be meshed into M × N grids, and M × N deflection values need to be calculated. Linear algebraic equations that incorporate the coefficients in Eq. (A3) can be used to calculate the deflection values as follows:  where − → v is the deflection vector (a column vector with M × N rows), R is the difference coefficient matrix (a matrix with M × N columns and M × N rows), and − → h and − → w are the final surface and subsurface topography vectors (column vectors with M × N rows), respectively. We employed periodic boundary conditions during the calculations, and 24 kinds of boundary conditions that should be considered (Eq. A2 and Fig. 13). Furthermore, compression of the sparse matrix R is recommended, since it contains very few nonzero elements.
14 Schematic diagram of the node operation for the central difference operators. a First-order differential operations. b Second-order differential operations. c Second-order mixed partial differential operations. d Laplace differential operations. Null nodes indicate zero coefficients