Open-source modular solutions for flexural isostasy: gFlex v1.0

Isostasy is one of the oldest and most widely applied concepts in the geosciences, but the geoscientific community lacks a coherent, easy-to-use tool to simulate flexure of a realistic (i.e., laterally heterogeneous) lithosphere under an arbitrary set of surface loads. Such a model is needed for studies of mountain building, sedimentary basin formation, glaciation, sea-level change, and other tectonic, geodynamic, and surface processes. Here I present gFlex (for GNU flexure), an open-source model that can produce analytical and finite difference solutions for lithospheric flexure in one (profile) and two (map view) dimensions. To simulate the flexural isostatic response to an imposed load, it can be used by itself or within GRASS GIS for better integration with field data. gFlex is also a component with the Community Surface Dynamics Modeling System (CSDMS) and Landlab modeling frameworks for coupling with a wide range of Earth-surfacerelated models, and can be coupled to additional models within Python scripts. As an example of this in-script coupling, I simulate the effects of spatially variable lithospheric thickness on a modeled Iceland ice cap. Finite difference solutions in gFlex can use any of five types of boundary conditions: 0-displacement, 0-slope (i.e., clamped); 0-slope, 0shear; 0-moment, 0-shear (i.e., broken plate); mirror symmetry; and periodic. Typical calculations with gFlex require 1 s to∼ 1 min on a personal laptop computer. These characteristics – multiple ways to run the model, multiple solution methods, multiple boundary conditions, and short compute time – make gFlex an effective tool for flexural isostatic modeling across the geosciences.

Theory to describe deflections of the lithosphere under loads has evolved significantly over the past 160 years (Watts, 2001). The development of this theory started with simple approximations of perfect buoyant compensation of loads by a lithosphere with no strength overlying a mantle of known density (Airy, 1855;Pratt, 1855). These approximations allowed surveyors to explain the observed lack of significant gravity anomalies around large mountain belts (cf., Göttl et al., 2009). While this theory, called isostasy, revolutionized the way topography was viewed on the Earth, more realistic solutions for isostatic deflections of the surface of Earth take into account the bending, or flexure, of a lithospheric plate of nonzero but finite strength. This strength may be defined as the elastic thickness, the effective thickness of a flawless plate of the equivalent strength, or as the flexural rigidity that is characteristic of a plate of a given thickness (see Eq. A7). By bending over distances of several tens to hundreds of kilometers, the lithosphere low pass filters a discontinuous surface loading field into a smoothed solid-Earth response.
Even though the early geological theories of Pratt (1855) and Airy (1855) focused on simple buoyancy, the differential equation basis for solving lithospheric bending already existed at that time. Bernoulli (1789) and Germain (1826, and earlier work) developed the first differential-equation-based theories for plate bending. Lagrange (1828) reviewed the prize that Germain won in 1811 for her work on elastic plate flexure, and, on realizing an error in the lumping of terms due to Germain's incorporation of an incorrect formula by Euler (1764), corrected it and produced the first complete flexure equation (see reviews by Todhunter and Pearson, 1886;Ventsel et al., 2002). Around the same time, Cauchy (1828) and Poisson (1828) better connected the theory of elasticity to plate bending problems. These works predated Kirchhoff (1850), who developed the classical or Kirchhoff-Love plate theory that remains in use today (Ventsel et al., 2002). While many further advances have been made (e.g., Love, 1888;Timoshenko et al., 1959) especially for structural and aeronautical engineering, it is the Kirchhoff-Love plate theory that has been used most widely for geological applications (e.g., van Wees et al., 1994). Comer (1983) tested classical Kirchhoff plate theory, which is a thin-plate theory that simplifies the plate geometry and therefore the mathematics required to solve for it, against a thick-plate theory of lithospheric flexure. While this thick-plate theory relaxes several approximations, its solutions are very similar to those for thin-plate flexure (Comer, 1983).
In the first half of the twentieth century, Vening Meinesz (1931,1941,1950) and Gunn (1943) applied analytical solutions of the plate theory of Kirchhoff (1850) to geological problems. They employed analytical solutions that relate the curvature of the bending moment of a plate of uniform elastic properties to an imposed surface point load, line load, or si-nusoidal load. These load solutions could be used to compute flexural response to any arbitrary sum of individual loads in either the spatial or spectral domain, due to the linear nature of the biharmonic flexure equation (Eqs. 1 and 2), and may be combined with a variety of boundary conditions (Watts, 2001).
Computational advances allowed discretized models to replace purely analytical solutions. These models fall into one of several categories. Many take advantage of the linear nature of the flexure equation for constant elastic thickness to superimpose analytical solutions of point loads (in the spatial domain) or sinusoidal loads (in the wavenumber domain) in order to produce the flexural response to an arbitrary load (Comer, 1983;Royden and Karner, 1984). Other models produce numerical solutions to the thin plate flexure equation by solving the local derivatives in plate displacement with numerical (mostly finite difference) methods (e.g., Bodine et al., 1981;van Wees et al., 1994;Stewart and Watts, 1997;Pelletier, 2004;Govers et al., 2009;Sacek et al., 2009;Wickert, 2012;Braun et al., 2013). Models in this latter category allow for variations in the elastic thickness of the plate, a factor of growing importance as variations in elastic thickness through space and time are increasingly recognized, measured, and computed (e.g., Watts and Zhong, 2000;Watts, 2001;Van der Lee, 2002;Flück, 2003;Pérez-Gussinyé and Watts, 2005;Tassara et al., 2007;Pérez-Gussinyé et al., 2007Tesauro et al., 2009;Kirby andSwain, 2009, 2011;Lowry and Pérez-Gussinyé, 2011;Tesauro et al., 2012bTesauro et al., , a, 2013Braun et al., 2013;Kirby, 2014). In spite of these efforts, the community currently lacks a robust, easy-to-use, generalized tool for flexural isostatic solutions that can be used by modelers and data-driven scientists alike.
Here I introduce a broadly implementable open-source package of solutions to flexural isostasy. This package, called gFlex (for GNU flexure), advances and makes more accessible an earlier model, generically called flexure (Wickert, 2012). gFlex has been released under the GNU General Public License (GPL) version 3 and is made available to the public at the University of Minnesota Earthsurface GitHub organizational repository, at https://github. com/umn-earth-surface/gFlex, and through the Python Package Index (PyPI). This allows for rapid collaborative editing of the source code and easy automated installation. It is written in Python (e.g., Rossum et al., 2012) for easy interoperability with a range of other programming languages, models, and geographic information systems (GIS) packages, and to take advantage of the numerical packages for Python that allow for much more rapid matrix solutions than would be typical with a more basic interpreted language (Jones et al., 2001;Davis, 2004;Oliphant, 2007;van der Walt et al., 2011). See Section 5 for further information on obtaining and running gFlex.
gFlex can solve plate flexure in two major ways (Fig. 2). First, it can produce analytical solutions to flexural isostasy generated by superposition of local solutions to point loads in the spatial domain (i.e., as a sum of Green's functions) (e.g., Royden and Karner, 1984). These use biharmonic equation for plate flexure with uniform elastic properties (Eqs. 1 and 2) (Bodine et al., 1981). Second, it can compute finite difference solutions for both constant and arbitrarily varying lithospheric elastic thickness structures. These solutions follow the work of van Wees et al. (1994), and hence Braun et al. (2013), except that gFlex does not incorporate terms for end loads but does include a wider range of implementable boundary conditions (Table 1). gFlex can be run as a standalone program with an input file, as a component of the indevelopment Landlab landscape modeling framework Tucker et al., 2013Tucker et al., , 2015 and by extension as a component within the Community Surface Dynamics Modeling System (CSDMS) (Syvitski et al., 2011;Overeem et al., 2013), or as a pair of add-ons to the Geographical Resources Analysis Support System (GRASS) Geographic Information System (GIS) (Neteler et al., 2012). The GRASS GIS implementation is particularly important, as it provides pre-built and standardized command-line and graphical interfaces and the ability to directly pull inputs from and compare solutions against field data in their native coordinate systems.

Methods and model development
Two solution types for flexural isostasy are provided in gFlex, and these are formulated for both one-dimensional (line load, assumed to extend infinitely in an orientation orthogonal to the line along which the equation is solved) and two-dimensional (point load) cases. The derivation that forms the basis for both of these is provided in Appendix Sect. A, and similar approaches to this derivation may be found in the work of Timoshenko et al. (1959) and Turcotte and Schubert (2002). The analytical and finite difference approaches are compared and shown to approximate each other well in Fig. 3. Table 1. Boundary conditions. Names provided here are the same as those used in the model, and b.c. stands for boundary condition. The first five can be selected for numerical solutions. The final one, NoOutsideLoads, is the outcome of superposition of analytical solutions, which allows the entire space to respond to local loads as if the 0-deflection boundaries were infinitely far away. In this notation, the subscript b indicates the boundary, generically. Where 0 and n are included as subscripts, i.e., for the mirror and periodic boundary conditions, these indicate boundaries at the first and last node of the model domain along a particular axis. Subscript x, which is a stand-in for x or y, is a variable distance to indicate the symmetry across a mirror boundary. Each of these boundary conditions requires a corresponding boundary condition for flexural rigidity.

Name
Mathematical Description Rigidity b.c.
dx 3 = 0 broken plate with a free cantilever end

Superposition of analytical solutions
The first solution type takes advantage of the linear nature of the analytical solution for flexure of a plate of constant thickness and elastic properties when subjected to a point or line load. These solutions may be superposed (i.e., summed) in space to compute the full flexural response. The second approach is to solve the equation for lithospheric flexure as a matrix equation by employing a finite difference scheme. This employs a sparse matrix elimination solver (e.g., Davis, 2004). The primary gFlex finite difference solution follows the approach of van Wees et al. (1994) to permit computations with steep gradients in flexural rigidity (Appendix Sect. A2), but gFlex also offers the discretization of Govers et al. (2009). The analytical solution imposes the assumption that scalar flexural rigidity, D, is uniform. This leads to biharmonic expressions for plate bending in one and two dimensions, respectively: Here, w is vertical deflection of the plate, q is the applied surface load, and ρ = ρ m − ρ f is the density of the mantle minus the density of the infilling material; see Fig. A1 for a diagrammatic description of all variables. The ρ term represents the feedback by which flexural subsidence can lead a depression to be filled by material, which leads to additional flexural subsidence. This can occur, for example, in a system that is fully underwater (e.g., an underwater volcano load) or one in which the depression is completely filled with sediments. If this infilling material is not uniform in density and/or spatial extent -for example, due to onlap or offlap of water along a shoreline -then one may solve this feedback instead via iteration, by solving for ρ f = ρ air ≈ 0 and adding water (or another load) to regions that match certain conditions after every cycle of the iteration.
The above equations are linearizable, and therefore can be solved by superposition of analytical solutions. In gFlex, this is done in the spatial domain on both structured grids and as a response to an arbitrarily placed set of point loads. Spectral solutions are possible (Stephenson, 1984;Stephenson and Lambeck, 1985) and efficient using fast Fourier transform algorithms (cf. Welch, 1967), but have not been implemented. The one-and two-dimensional solutions for lithospheric flexure take the form of an exponentially damped sinusoid. In one dimension, this is represented by the following expression: Here, the i subscript indicates that this is the response to a line load at a single x position, x i . α 1-D is the one-dimensional flexural parameter, defined by Vening Meinesz (1931) (following Hertz, 1884): The significance of the flexural parameter is that the flexural wavelength, λ α is related to the flexural parameter as λ α = 2π α. The distance from a point load to the first flexural bulge (forebulge) that it creates around its local depression, for example, is a flexural half-wavelength, π α. This nature of plate bending as an exponentially decaying periodic function can be seen most easily in the one-dimensional analytical (constant T e ) solution in Eq. (3). Brotchie and Silvester (1969) derived that the exponentially damped sinusoid due to a point load in two dimensions should be expressed by kei (Abramowitz and Stegun, 1972), which is the zeroth-order Kelvin function that satisfies the equation ker(r) + ikei(r) = K 0 re πi/4 , where K 0 is the zeroth-order modified Bessel function of the second kind. This function was defined by Lord Kelvin to solve for electrical current density in a circular wire with an applied oscillating (alternating) current (Barron and Barron, 2012, Appendix 5), and its solution has been broadly applied to the two-dimensional bending of a plate (e.g., Timoshenko et al., 1959;Lambeck, 1981;McNutt and Menard, 1982;Watts, 2001).
The subscripts i, j indicate that this is the flexural response to a single point load at the x and y positions x i and y j . The two-dimensional flexural parameter, α 2-D , contains D instead of 4D in the numerator because it does not need to include implicit loads and deflections along the y orientation that are required in the one-dimensional line-load plate bending case.
Lithospheric flexure calculated by superposition of analytical solutions can be represented as a simple sum across all line loads q l or point loads q p : For a given elastic thickness, each flexural response to a line or point load is similar in shape, but different in amplitude. Therefore, I optimize solution speed by pre-calculating the flexural response to a unit load in the center of a template array. This pre-calculated unit deflection array has twice the linear dimensions of the solution array, and is subsampled and re-scaled to compute the distributed response to each cell in the grid that contains a load. This technique works for all for rectilinear grids with uniform x and y grid spacing, though the x and y grid spacing do not have to be equal to one another. A similar optimization is possible for onedimensional solutions, but these are so rapid that this has not been found to be necessary. Within gFlex, this solution type is termed SAS, which stands for superposition of analytical solutions.
The analytical solution response to point or line loads can also be computed for a scattered set of loads and a scattered (and not necessarily the same) set of points at which the flexural response is calculated. This solution type is termed SAS_NG, which stands for, superposition of analytical solutions: no grid. Because it lacks the grid uniformity that permits the a solution template to be used, its computational time is not optimized in this way (Sect. 2.5).

Finite difference solutions
Finite difference solutions in one and two dimensions employ Eqs. (A19) and (A20), respectively. For these solutions, dx and dy may differ from one another, but each must be constant. First, for the one-dimensional solution, the expansion of Eq. (A19) is The two-dimensional solution is based on an expansion of Eq. (A20) (van Wees et al., 1994): These equations are discretized using a second-order accurate centered finite difference approximation (Fornberg, 1988, Table 1). Finite difference solutions in two dimensions may also be generated following the solution and discretization of Govers et al. (2009), which produces solutions for a more limited range of flexural rigidity variations.
The finite difference solution is computed as a linear matrix equation, where A is a sparse matrix of operators from a linear decomposition of Eqs. (A19) or (A20), W is a vector of deflections (typically unknown), and Q is a vector of imposed loads (typically known). It is solved directly by using the sparse LU (lower upper) factorization package UMFPACK (unsymmetric-pattern multifrontal package) (Davis, 2004) or, at the user's choice, iteratively with one of the many solvers that are available with the SciPy (Scientific Python) package (Jones et al., 2001).

Boundary conditions
gFlex supports a number of boundary conditions, and these are summarized in Table 1 and schematically drawn in Fig. 4. The finite difference (sparse matrix) numerical solutions can freely define any combination of no-displacement-and-no-slope (0Displacement0Slope), nobending-moment-and-no-shear (0Moment0Shear), no-slopeand-no-shear (0Slope0Shear), and mirror boundaries. Periodic boundaries may be mixed with any combination of the aforementioned boundary conditions, with the requirement that they exist on both sides of the deflection array, as having (for example) deflections at the west end of the array sensitive to loading and deflections to the east, but with those on the east not, in turn, sensitive to the west, being nonsensical. Superposition of analytical solutions naturally produce a 0displacement boundary at infinite distance from each point load (NoOutsideLoads). This can be seen by solving Eqs. (3) and (5) as x → ∞ and y → ∞. Each of these boundary conditions can be related to geological processes or locations that one may wish to model (Fig. 5). The 0Displacement0Slope (or clamped) boundary condition ( Fig. 4a) may be used to approximate a NoOutsideLoads case for the finite difference solutions (Fig. 3). When placed one flexural wavelength away from a point or line load, the surface displacement should, for a plate of constant elastic thickness, be ∼ 0.2 % of that at the point of maximum deflection, which is negligible compared to most sources of geological error. It is conceivable that a difference in elastic thickness in a continuous plate may exist that is so great that the thicker plate can be approximated to not bend; a 0Displace-ment0Slope boundary condition may also be used to simulate this, though one must debate whether to do this or to compute the flexural response across a plate with a prescribed elastic thickness variability.
The 0Moment0Shear boundary condition (Fig. 4b) means that the edge of the plate is completely free to flex, like the cantilevered end of diving board. This is appropriate for places in which the elastic thickness of the lithosphere goes to zero. Such broken-plate boundary conditions have been used in analytical solutions to simulate flexure of the lithosphere beneath the Hawaiian volcanoes, where heating significantly weakens the lithosphere (Wessel, 1993). This approximates the (zero-dimensional) single point discontinuity of a hotspot as a (one-dimensional) line boundary condition. A broken-plate solution has also been used for zones beneath mountain ranges where sufficient deformation may weaken the lithosphere (Stewart and Watts, 1997), and may be best suited for continental rift zones (Burov et al., 1994), as these closely approximate a linear discontinuity in an otherwise thick lithosphere. In all three of these cases, the lithosphere should lose strength as it approaches the boundary condition. For this reason, 0Moment0Shear is implemented only for the finite difference solution, which allows for spatial variations in elastic thickness.
The 0Slope0Shear boundary condition (Fig. 4c) may be considered to be a flat clamp on the boundary of the plate that may be freely moved upwards or downwards. While it may require creative thought to uncover a geological process that holds a plate edge flat but allows it to move freely in the vertical, this boundary condition can also be used at an appropriate distance away from the load(s) to approximate a NoOutsideLoads boundary for a finite difference so-lution, though typically the 0Displacement0Slope boundary provides a closer match.
The periodic boundary condition (Fig. 4d) wraps one side of the model around to the other side such that they form an infinite loop. To visualize this, one may imagine taking a paper map and taping either the east and west sides together or the north and south sides together, such that the flexure induced by loads on one edge is continuous with load-induced flexure on the opposite edge. Elastic thickness and loads both wrap around this boundary, making it possible to, if one is not careful, create sudden jumps in elastic thickness at the edge of the model. This takes somewhat longer to solve (Fig. 6c), but can be useful to compute a flexural response to the load of a long mountain belt by modeling just a limited region perpendicular to the strike of the range crest and allowing this slice to infinitely repeat in the range-crest-parallel orientation; at the limit of a very narrow slice of model space, this approaches the one-dimensional line-load solution. If a future model of lithospheric flexure relaxes the current assumption in gFlex that dx and dy may be different but must be constant in space, the periodic boundary condition should enable a finite difference flexural model to be employed on a closed surface, such as a sphere, enabling full global modeling. This is, to the best knowledge of the author, the first time that a periodic boundary condition has been implemented for lithospheric flexure.
The mirror boundary condition (Fig. 4e) reflects the elastic thickness and load structure across a plane of symmetry at the boundary. This may be used to speed a solution where a plane of mirror symmetry may be implied, which is important for large grids or where gFlex is used as part of a coupled set of numerical models (e.g., through CSDMS: Syvitski et al., 2011;Overeem et al., 2013;Peckham et al., 2013;Tucker et al., 2015). Example usage cases include topographic unloading by erosion of a symmetrical mountain range (Fig. 5c, and 5d), isostatic adjustment under a symmetrical ice cap, and emplacement of a volcanic load. The latter two cases often have fully radial symmetry, and therefore may be placed at the corner of the solution array with mirror boundary conditions on both adjacent sides to further limit the needed computational area. This is also to the best knowledge of the author the first application of a mirror boundary condition to modeling of lithospheric flexure, which is surprising considering its potential utility.
The names of the boundary conditions are based on their effects on deflections, w, but solutions also require boundary conditions to be placed upon the flexural rigidity, D; these are listed in Table 1. For the 0Displace-ment0Slope, 0Slope0Shear, and 0Moment0Shear deflection boundary conditions, a 0-curvature flexural rigidity boundary condition has been chosen. This allows for near-boundary gradients in flexural rigidity to be assumed to continue outside the computational domain. As noted above, mirror and periodic boundary conditions are applied to the rigidity field as well. For the analytical solutions, the approximation is an infinite plate of constant elastic thickness.
In two-dimensional solutions, boundary conditions meet at corners. Where a boundary condition meets another of the same boundary conditions at the corner, the two generate a continuous boundary condition that includes the corner of the array. This is always the case for the analytical solutions with implicit NoOutsideLoads boundary conditions. Where mirror or periodic boundary conditions meet themselves at corners, these produce doubly reflecting or doubly periodic boundaries; if every boundary is mirror or periodic (necessary in the latter case as periodic boundary conditions must always exist as pairs on opposite sides), these generate an infinite tessellated plane of loads and elastic thicknesses. Some boundary conditions in gFlex can work harmoniously with others. Periodic and mirror boundary conditions propagate 0Moment0Shear, 0Slope0Shear, and 0Displacement0Slope boundary conditions that exist orthogonally to them. Where mirror and periodic boundary conditions intersect at a corner, the periodic boundary condition will propagate the mirror boundary to ±∞. Those boundary conditions that do not reflect or repeat the effects of the other boundary conditions do not share the corners equally: in gFlex, 0Displacement0Slope boundary conditions dictate all corners where they meet other boundary conditions, forcing them to remain fixed at 0; physically, this means that the clamp of the 0Displace-ment0Slope boundary condition continues through the edges of the perpendicular boundaries. 0Moment0Shear boundary conditions were chosen to control the corners where they meet 0Slope0Shear boundary conditions, as the 0Mo-ment0Shear boundary condition has been recognized in geological work (e.g., Wessel, 1993;Burov et al., 1994;Stewart  The gridded superposition of analytical solutions (SAS) scales with the total number of grid cells times the number of cells with loads, as this is the total number of computations that must be made. (c) Finite difference solutions are computed with sparse matrices with dimensions equal proportional to the grid dimensions, squared, and therefore scale with the number of total grid cells. All of the solution time relationships are close to linear except for the two-dimensional finite difference solutions, due to the added complexity of their finite difference stencil. Many fits are to a subset of the data to avoid those solutions that are so rapid that the amount of time required for the non-solver portions of the code becomes significant. All marker symbols are semi-transparent, meaning that darker symbols than those that appear in the legend imply additional data points underneath. and Watts, 1997), while the 0Slope0Shear boundary condition has not.

Discontinuities and limit as T e → 0
Two notable issues inherent to the finite difference solutions and the treatment of a continuous plate become apparent as T e → 0. The first is that a region of T e = 0 must have a width of at least three cells to produce the expected local isostatic equilibrium; this is a result of numerical diffusion in the central difference discretization provided in Eqs. (9) and (10). The second is that because any region of 0 elastic thickness will enter isostatic equilibrium with its local loads and not be affected by nonlocal effects; if this region lies along the edge, it will ignore all boundary conditions. If a T e = 0 region along a boundary is ≥ 3 cells wide, it imposes a 0Dis-placement0Slope boundary condition on the interior cells; smaller regions of T e = 0 will allow some information on the ultimate boundary condition to leak through via numerical diffusion.
These issues are important to note, but unlikely to be important in most cases of gFlex. First, discontinuous transitions to zones of T e = 0 may also be modeled by segmenting the inputs into multiple arrays, running gFlex for each array with a 0Moment0Shear (broken-plate) boundary condition applied to the model domain edges representing the discontinuities, and then recombining the outputs into a continuous displacement field. Second, and more importantly, the conditions for broken-plate or T e = 0 solutions to be required are rare on Earth. Elastic thickness of 0 implies that there is no elastic lithosphere, and a broken-plate solution implies that there is no shear between adjacent lithospheric blocks. These conditions are most likely to be met in rift zones, though even these have some nonzero thickness of brittle crust. End loads, which are not currently included in gFlex, could be used in combination with a 0Moment0Shear boundary condition to better parameterize faults (e.g., van Wees et al., 1994) and expand the utility of gFlex. However, the typical case for which gFlex is designed involves glacial-isostatic adjustment, large-scale water loads, sedimentary basin development, large-scale erosional unloading, and other processes that extend across a swath of heterogeneous lithosphere that may contain many faults. In these cases, it has been found to be sufficient to simply characterize a variable field of finite elastic thickness across the domain, where elastic thickness falls around fault zones (e.g., Manríquez et al., 2013).

Model benchmarking
A set of tests was performed to measure the speed at which gFlex computes solutions. In these tests, an elastic plate that is 1000 km long (one-dimensional and two-dimensional) and 1000 km wide (two-dimensional) is subjected to a square load at its center that ranges from 100 km to the full 1000 km on each side. This load places a normal stress of 9 702 000 Pa on the surface, which is equal to 300 m of mantle material (3300 kg m −3 ). In these scenarios, there is no assumed infilling material (ρ f = 0). gFlex computed solutions for uniform rectilinear grids of increasing size using gridded and ungridded superposition of analytical solutions (SAS and SAS_NG, respectively) and finite difference (FD) methods. All boundary conditions (Table 1 and Fig. 4) were tested, though not in combination. The finite difference solutions include scenarios with both constant (25 km) and variable (10-40 km) effective elastic thickness, with the latter varying sinusoidally over a wavelength of 500 km such that the plate contains two full T e cycles. In the two-dimensional case, T e varies in both dimensions to produce a smoothed checkerboard pattern of elastic thickness. Finite difference solutions reported employ the direct solver UMFPACK (Davis, 2004), as it has been better tested in gFlex than the iterative solution methods and is therefore the default solver. Fig. 6 displays computation time for all of the benchmarking tests, and Fig. 7 is a comparison of the SAS_NG, SAS, and FD solution techniques for the case in which every point at which the solution is calculated also contains a nonzero load. These solution times do not account for file input or output or graphics generation. They do include the initialization time for the solution steps of gFlex; therefore, a number of the power-law fits to solution time do not include the times calculated with the smallest arrays, for which initialization time is a significant fraction of the total model runtime.
The factors that determine computation time are the solution method and the inclusion of periodic boundary conditions. While the SAS_NG method scales the best with increasing grid size, it is so much slower than the other methods for standard model-run grid sizes that it will not often exceed their speed. The finite difference method is the fastest if every cell contains a load, but can become slower than the an-Number of grid cells Figure 7. Comparison between solution methods where every cell in the domain contains a load. The ungridded superposition of analytical solutions (SAS_NG) scales best but in these tests is the slowest. It can, however, be faster when fewer cells contain loads. Some fits are to a subset of the data to avoid those solutions that are so rapid that the amount of time required for the non-solver portions of the code becomes significant. All marker symbols are semitransparent, meaning that darker symbols than those that appear in the legend imply additional data points underneath.
alytical methods if only a few cells contain loads, as analytical methods must make one set of calculations across the grid per load. Standard runtimes are between a fraction of a second and a few minutes on a personal laptop computer (Dell XPS 13 Developer Edition running Ubuntu 14.10) (Figs. 6 and 7).

Model interfaces and coupling
Some users of gFlex may want to run a single calculation, while others may want to produce many solutions as part of a larger-scale numerical modeling exercise, such as an inversion for elastic thickness or coupling with another model. Therefore, four different methods to use gFlex have been prepared: 1. standalone, with input files; 2. as part of a Python script; 3. driven by GRASS GIS (Neteler et al., 2012) to simplify integration of geospatially registered data with the lithospheric flexure model; 4. as a component for the CSDMS framework (Syvitski et al., 2011;Overeem et al., 2013;Peckham et al., 2013), including its tight integration into Landlab, a CSDMS-led Python-based Earth-surface modeling framework that is currently being developed Tucker et al., 2013Tucker et al., , 2015. GRASS GIS integration is also possible for model coupling using Python, including efforts that use the Landlab framework.

Standalone with input files
Some users may want to employ gFlex as a single calculation, for example to calculate the flexural response to a set of loads generated by a sedimentary deposit that was measured in the field. The user prepares an input file of model settings, an input ASCII grid of loads, and, should the elastic thickness be nonuniform, an input ASCII grid of lithospheric elastic thicknesses. Outputs from this mode of running gFlex include an ASCII grid of surface deflections and a set of plots of surface deflections and loads.

As part of a Python script
gFlex may also be imported as a Python module to be run either as a standalone simulation or as a component in a multi-model integration effort. This allows it, for example, to be a part of a flexural backstripping toolchain or a model of glacial-isostatic adjustment. Backstripping calculations may be performed by simply removing the sedimentary load (Roberts, 1998), or, in the case of a foreland basin, by inverting for the mountain belt loading history and lithospheric elastic thickness that would be required to produce the basin (Ballato et al., 2016). A programmatic approach is also useful for scenarios in which material infills a depression, but not over the whole domain and/or not with uniform density. While the flexure equations require that ρ f be constant, a more flexible way to solve for the effect of infilling material is to compute flexural response with ρ f = ρ air ≈ 0, add loads based on some set of rules, and then re-calculate flexure iteratively until convergence is achieved. This can occur in regions with a complex set of sedimentary deposits (see also Watts et al., 1982;Watts, 2001) and/or to be used for seawater loading across a shoreline (see also Mitrovica and Milne, 2003).

Driven by GRASS GIS
gFlex is also prepared for integration with the open-source geospatial software GRASS GIS (Neteler et al., 2012) as two add-ons or extensions named r.flexure and v.flexure, which are raster and vector operations, respectively. As GRASS GIS is a map-based application, r.flexure and v.flexure employ two-dimensional solutions (both analytical and finite difference), though future extensions of these modules to compute flexure from line loads along chosen onedimensional profiles would be possible. r.flexure can use the finite difference or SAS solution methods, whereas v.flexure exclusively uses the SAS_NG solution method to take advantage of its ability to produce solutions for an arbitrary scatter of point loads. Advantages of GRASS GIS include 1. full integration within a geospatially registered environment, meaning that data can be used directly as model inputs, and that model outputs may be compared against data; 2. a documented and standardized command-line interface; 3. a pre-built and standardized graphical user interface (GUI).
The graphical user interface is incorporated into the GRASS GIS wxPython GUI (Landa, 2008;Neteler et al., 2012), and this is particularly helpful for researchers, who are not as accustomed to command-line interaction with computers to use gFlex with their data. For computer modelers, the GRASS GIS coupling may be used to support broader model coupling and data-model integration efforts (see, for example, Srinivasan and Arnold, 1994).

Modeling frameworks
CSDMS (broadly) and Landlab (in particular) both include methods for integrating modular blocks of code as part of their respective efforts towards the community-wide goal to make modeling of Earth systems less time intensive and more streamlined (Voinov et al., 2010;Syvitski et al., 2011;Overeem et al., 2013;Peckham et al., 2013;Hobley et al., 2013;Tucker et al., 2013Tucker et al., , 2015. gFlex is included as a modular component of the still-in-development Landlab Earthsurface modeling framework Tucker et al., 2013Tucker et al., , 2015. Landlab integration provides wrapping with the CSDMS Basic Model Interface (BMI) and Component Model Interface (CMI) using the CSDMS Standard Name construction conventions (Peckham et al., 2013). The standard interfaces provided by both of these modeling frameworks will streamline model coupling that uses gFlex and help to prevent duplication of effort in building plate bending models. Furthermore, the inclusion of gFlex in Landlab will allow numerous Earth-surface systems to be modeled more precisely (Fig. 1).

Application example: Iceland
As a first example to utilize both the ability of gFlex to generate solutions with variable effective elastic thickness and its incorporation into GRASS GIS, gFlex is used along with a simple and efficient GIS-enabled glacier and ice cap model modified from the work of Colgan et al. (2015) to model a hypothetical expansion of the Iceland Ice Cap. While the importance of flexural isostasy in ice dynamics modeling has long been well-known (cf. Cuffey and Paterson, 2010), the author knows of no dynamic ice model that runs with a variable elastic thickness lithosphere, making this possibly the first such exercise. Earth's crust at Iceland has been built by Flexural isostasy with a constant 3.7 km elastic thickness (c) (following Hubbard, 2006) reduces ice cap extent and causes some interior ice thickening when compared to the case without flexure (b), as the ice cap conforms to the bowl-shaped depression that it creates. Deformation in the case with variable elastic thickness (f) is focused along the ridge and extends farther on the southwestern side that has greater elastic thickness, and modifies the topography of western Iceland (low elastic thickness) to produce spatially variable ice thickness changes (e, h).
the unique intersection of the Iceland hotspot and the Mid-Atlantic Ridge, which together produce a weak lithosphere with spatially variable elastic thickness, resulting in shortwavelength variability in the solid-Earth response to loading.
Here I test the two-way coupling between ice dynamics and solid-Earth deformation and the differences in steady-state ice caps that are produced in a modest climate change and ice cap extent scenario. This coupled ice dynamics and flexural isostatic model of Iceland requires four input components: the elastic thickness structure around Iceland, the modern topography of Iceland, the modern surface temperature field of Iceland, and modern precipitation rates across Iceland. The ice cap model used here (cf. Colgan et al., 2015) employs a shallow-ice approximation with basal sliding as a linear function of driving stress, which is intentionally much simpler than the modeling approach (Hubbard, 2006) that Hubbard et al. (2006) used to model the Last Glacial Maximum (LGM) Iceland Ice Cap. This is because the goal here is to show schematically the importance of including lateral variations in elastic thickness on the reconstructed thickness of an ice cap for a given paleoclimate, with less emphasis on actually reproducing any particular extent of the Iceland Ice Cap.
The elastic thickness structure under Iceland, in this schematic example, is related to the age of the oceanic crust. Calmant et al. (1990) related elastic thickness to the age of the lithosphere with the simple equation that results from the square-root time dependence of lithospheric cooling via thermal conduction (cf. Stein and Stein, 1992): Here, T e is given in kilometers and the age of the lithosphere, t, is given in millions of years. As continental material also exists within the computational window, the elastic thickness map of Tesauro et al. (2012a, b) is used for all subaerial landmasses. Across the continental shelves, the oceanic-crustbased map and the map from Tesauro et al. (2012a, b) are blended using spline interpolation within GRASS GIS (Neteler et al., 2012). The regional age of oceanic crust is provided by Müller et al. (2008), but their map indicates that even crust at the ridge in Iceland has an age of 6-7 Ma, resulting in a greater computed effective elastic thickness than would be expected based on the presence of the ridge or from heat flow data (e.g., Flóvenz and Saemundsson, 1993). While the structure of Iceland is certainly more complicated than the simpler parts of the ridge due to the effects of the hotspot and its tectonic environs (e.g., Watts and Zhong, 2000;Foulger, 2006), the assumption here is that the lithospheric effective elastic thickness structure due to the ridge is as if young crust continued along the Mid-Atlantic Ridge through all of Iceland, and the elastic thickness map (Fig. 8i) was modified to approximate this for the sake of this example.
The underlying digital elevation model, GEBCO_08 (British Oceanographic Data Centre , BaODC), includes the modern ice caps on Iceland, but these are already flexurally compensated and are small compared to the ice cap modeled here. While their removal would improve reconstructed ice discharge, they are ignored due to the schematic nature of this modeling effort.
Modern temperature and precipitation fields are from the Monthly NOAA-CIRES (National Oceanographic and Atmospheric Administration-Cooperative Institute for Research in Environmental Sciences) 20th Century Reanalysis (V2) by Compo et al. (2006Compo et al. ( , 2011 (for further background on their methods, see Whitaker and Hamill, 2002). These provide twentieth century mean conditions on a 2 • × 2 • latitude-longitude grid (temperature) or a 94 × 192 Gaussian grid (precipitation). These were cast as point data and interpolated using splines in GRASS GIS (Neteler et al., 2012). Prior to this spline interpolation, temperature was projected to sea level using the mean cell elevation (with lapse rate of 4.7 • C km −1 ), following Anderson et al. (2014) for ice caps; after interpolation, the resultant temperatures were then interpolated up to their respective surface elevations using the same 4.7 • C km −1 lapse rate. Although not all of the Icelandic surfaces are covered in ice at present, this rule was prescribed uniformly for the sake of a schematic model.
Three experiments were run: one with no flexure, one with flexure using a constant elastic thickness of 3.7 km (following Hubbard, 2006, and assuming E = 65 GPa), and one in which the full spatially variable flexure was used. In each of these runs, temperature was reduced from its present value by 5 • C and ice expanded to cover an area approximately equal to the currently subaerially exposed continent, approximately consistent with the previous modeling results of Hubbard et al. (2006) and with a temperature change that is much less than the LGM drop of 10-13 • C that was predicted to cause ice to spread onto the continental shelves as well . Mass balance was simulated by a positive degree-day melt model. June, July, and August temperatures were used to compute ablation, with a melt factor of 6 mm d −1 K −1 . Precipitation was held constant and all precipitation was assumed to contribute to positive mass bal-ance. Each scenario was run for 4000 years to reach full glacial and isostatic equilibrium, with isostatic equilibrium being assumed to occur instantaneously to facilitate more rapid computation of the equilibrium solution.
The results in Fig. 8 summarize the experiments. Figure 8c and f show the modeled flexural isostatic deformation and Fig. 8b and e show the deviation from the case with no isostasy; each of these pairs is for constant and variable elastic thickness, respectively. Figure 8h shows that with variable elastic thickness (Fig. 8i), ice thickness variability is concentrated where lithospheric elastic thickness is low.
The example of isostatic response to ice advance in Iceland is just one possibility of a feedback between an Earth-surface (or other geological) process and flexural deformation. Further such scenarios involving, for example, orogenesis and foreland basin formation (in settings such as that studied by Ballato and Strecker, 2014), rifting (Braun et al., 2013), and river delta morphologic evolution (Kim et al., 2006), will improve our understanding of the dynamic interactions between Earth's surface and subsurface (e.g., Braun et al., 2013).

Code availability
gFlex is available from the University of Minnesota Earth-surface GitHub repository at https://github.com/ umn-earth-surface/gFlex. It runs on Linux, Windows, and Mac computers running Python 2.(X ≥ 7).Y. It may be downloaded as an archive that is a snapshot of the state of the code, or cloned into an updatable copy of the software on the computer of an end-user. Version 1.0, described in this paper, is stored at https://github.com/umn-earth-surface/gFlex/ releases/tag/v1.0. gFlex is also stored on the Python Package Index (PyPI) at https://pypi.python.org/pypi/gFlex for easy automated download and installed with the commandline tool pip. gFlex documentation is available in its file README.md that is displayed at the main GitHub repository page, and some additional information is presented at the gFlex CSDMS Wiki page at http://csdms.colorado.edu/ wiki/Model:GFlex.
Interfaces to GRASS GIS and Landlab are available from their respective repositories. The GRASS GIS interface works with GRASS GIS 7.X and can be downloaded and installed automatically with the g.extension tool within GRASS (Neteler et al., 2012) or be downloaded through the subversion repository at http://trac.osgeo.org/grass/browser/ grass-addons/grass7. The Landlab interface is located in the Landlab GitHub repository at https://github.com/landlab/ landlab/tree/master/landlab/components/gflex. Both require a locally installed copy of gFlex to run.

Conclusions
gFlex is a new, open-source, easy-to-use model to compute isostatic deflections of Earth's lithosphere with uniform or nonuniform flexural rigidity due to arbitrarily distributed surface loads. It can be run as a standalone model through a configuration file, a Python module, a component in the Landlab and CSDMS community modeling frameworks, or via one of two GRASS GIS add-ons for a direct link to geospatial data. Its open-source code base may be updated and improved by the community, it may be easily installed using automated tools, and it is poised to be coupled with other models in efforts to understand interactions between multiple components of the Earth system. These attributes all embody my primary aim in creating gFlex: to provide an accessible set of flexural isostatic solutions for work across the geosciences by field scientists and modelers alike. Plates resist bending (i.e., flexure) through fiber stresses that develop in response to loading-induced deformation. In this appendix, the background of the theory is provided by an abridged derivation of plate flexure, which provides the background for the assumptions and solution methods employed in both the analytical and finite difference one-dimensional and two-dimensional solutions. Components of the theoretical background are also relevant for the various boundary condition options introduced in the main text.
A derivation of flexural response to a load can be subdivided into two components. The first is the bending moment, which describes the internally generated torques that resist bending. The second is the relationship between the bending moment and the imposed load.

A1 Bending moments
The bending moment of a plate, M, is the resistance of the plate to bending. This resistance exists because when a plate of > 0 thickness is bent, layers within the plate on the inside of the curve are placed under compression and layers within the plate on the outside of the curve are placed under tension. These fiber stresses are denoted σ x x in the along-plate coordinate system (x , z ) depicted in Fig. A1, and cause each infinitesimal layer of the plate to act like a spring that resists plate bending.
Classical (Kirchhoff-Love) plate theory is derived using an approximation of cylindrical bending (cf. Love, 1888).
Over short distances, the bent plate is assumed to follow the arc of a circle (Fig. A1). Arc length, s, is the product of the radius of curvature, r c , and the angle over which the arc is defined, θ.
The layer halfway between the top and the bottom of a homogeneous plate experiences no net extension or shortening during bending. This midpoint layer is therefore taken to be the reference radius of curvature, r c = r 0 , of a plate that extends from r = −T e /2 to r = T e /2, where T e is the effective elastic thickness of the plate. Flexural isostatic deflections are small compared to the length scale over which they occur, meaning that r 0 T e /2, and therefore approximates the true radius of curvature regardless of through-plate position z . To calculate the range of arc lengths, s, that exist above and below the reference layer at r 0 , one can note that Eq. (A1) describes a linear relationship between arc length and radius of curvature. Therefore, it is possible to use the definition of strain and Eq. (A1) to define the fiber strain in each layer as a function of its distance from the midpoint. The normal strain along the x orientation, ε x x , is given by Here, z is defined to be zero at r 0 and θ is held constant and therefore cancels out. Sign conventions are unimportant due to the symmetry of the problem above and below the midpoint layer (Fig. A1). As radius of curvature decreases, curvature increases: The long horizontal length-scales involved in flexure problems result in a small slope, (dz/dx) 1, and this squared becomes so much smaller than 1 that the denominator on the right-hand side ≈ 1. This small slope also allows the small angle approximation to be used, meaning that (x , y ) ≈ (x, y). As noted above, sign convention is unimportant -half of the plate experiences tension and the other half experiences compression -so taking the absolute value, which is included in the definition in order to maintain a positive radius of curvature, becomes unnecessary. Substituting w for z and r 0 for r c , removing the slope-related term in the numerator of Eq. (A3), and combining it with Eq. (A2), provides the magnitude of fiber strains as a function of distance from the midplane in the plate and curvature of the plate: Equation (A4) becomes important in the final step to define the bending moment because it relates fiber strains directly to deflections that can be measured and/or modeled. The bending moment itself, M, balances the torques generated by plate flexure. It therefore describes the resistance of the plate to bending, and is defined as the sum through the thickness of the plate of all fiber stresses σ x x times their respective lever arms z (cf. Turcotte and Schubert, 2002).
It is possible to rewrite this in terms of strain instead of stress via an elastic constitutive relationship (Hooke's Law), σ xx = Eε xx + νσ yy . Here, E is Young's modulus, which is a generalized spring constant that typically ranges between 10 10 and 10 11 Pa for rock (Turcotte and Schubert, 2002, p. 106), and is 65 GPa by default in gFlex. ν is Poisson's ratio, which describes how much material tends to extend (or shorten) in one direction when shortened (or extended) in another, and is commonly taken to be 0.25 for the lithosphere. An analagous equation, σ yy = Eε yy + νσ xx , exists in the y orientation. The stress required for these additional strains reduces the strain in a given orientation by a factor of 1/(1 − ν 2 ): In the one-dimensional case, ε yy = 0. Both E and ν lie outside of the integral because they are assumed constant over z . It is possible to solve for the bending moment in one dimension by using Eq. (A4) to replace ε xx in Eq. (A6). As the orientation of the curvature (d 2 w/dx 2 ) is orthogonal to the direction of integration, the integral is simple to solve and results in the solution for the bending moment: The terms to the left of the derivative define the scalar flexural rigidity, D: ET 3 e 12(1 − ν 2 ) .
As D is the key parameter that controls flexural response, and is a function of T e , E, and ν, gFlex contains the additional simplifying assumption that E and ν are uniform constants. This permits variations in scalar flexural rigidity to map to variations in effective elastic thickness via Eq. (A8). It prevents overparameterization in gFlex, and implicitly states the assumption that changes in the effective elastic thickness of the lithosphere, cubed, are more significant than changes in Poisson's ratio, squared, or Young's modulus.
To generalize the bending moment of a plate that is loaded in two dimensions, one can start by writing a vector of curvatures, κ (cf. van Wees et al., 1994): The first term in κ, ε xx /z , scales with the x-directed normal strain that is also part of the one-dimensional solution. ε yy /z is the equivalent for y-oriented fiber normal strains, and γ xy /z is the fiber shear strain term that accounts for torsion.
The flexural rigidity must also be defined in three dimensions, and is defined here following linear elasticity: Using Eqs. (A10) and (A11), one can define the bending moment as Solving Eq. (A12) for only the upper (left) terms allows one to recover Eq. (A7), the one-dimensional case.

A2 Force and torque balance
A static lithospheric plate that experiences the downward force of an imposed load, qdx, responds with differential shear forces, V and V +dV , that develop in response to bending. Further vertical normal stresses that influence plate flexure are generated by the sum of the buoyant restoring force of displaced mantle, ρ m g(−w), and additional driving forces by any surface loads that fill the flexural depression, ρ f gw. Summed together, these form the additional term − ρgw, where ρ = (ρ m − ρ f ) The total vertical force balance is therefore These shear forces generate torques that must be balanced in turn by the bending moments. Here I explicitly ignore end loads because they are not part of the numerical solution in gFlex, which was designed with surface loads in mind, though they are straightforward to include (see van Wees et al., 1994;Turcotte and Schubert, 2002;Braun et al., 2013). The resultant torque (τ ) balance is: After noting that dV V , Eq. (A15) simplifies to This can be rearranged to define the shear force as the negative slope of the bending moment, which in turn is proportional to the curvature of the deflection (Eq. A7): This observation is key to defining the 0Moment0Shear and 0Slope0Shear boundary conditions (Table 1 and Fig. 4). Equations (A14) and (A17) can be combined by substituting V in Eq. (A14) to relate the bending moment and deflection to the imposed load, q.
To solve the two-dimensional case, one can follow van Wees et al. (1994) in using the differential operators for curvature inκ (Eq. A9) to generalize the one-dimensional flexural solution. This is then combined with the infill and buoyancy term. Written in compact form, the two-dimensional flexural isostatic equation is: (A20)