Toolbox for Analysis of Flexural Isostasy (TAFI)—A MATLAB toolbox for modeling flexural deformation of the lithosphere

In this paper we present the Toolbox for Analysis of Flexural Isostasy (TAFI). TAFI is a suite of tools for the MATLAB computing platform that supports forward modeling of flexural subsidence and uplift of the lithosphere in response to tectonic loading. TAFI also calculates the gravity anomaly due to flexure of the lithosphere, and allows for easy comparison of modeled deflection and gravity curves with observations to facilitate a rapid interactive search for the best-fitting flexural model.

Vertical tectonic loads emplaced on the lithosphere may be in the form of static surface loads such as topography, static subsurface loads due to lateral variations in lithosphere or asthenosphere density, or dynamic subsurface loads created by tractions at the base of the lithosphere imposed by flow in the underlying mantle. These tectonic loads result in vertical deformation, typically involving subsidence in the vicinity of the load (assuming the load represents an addition of weight to the lithosphere) and uplift in surrounding areas (the peripheral “bulge”) (Fig. 1). The magnitude and wavelength of the flexural deflection depends on the magnitude of the load and the strength of the lithosphere, with greater strength resulting in more broadly distributed (larger-wavelength) deformation. Flexure results in variations in the regional gravity field due to deformation of density interfaces within the lithosphere. Such density interfaces may include the crust-mantle boundary, the interface between the sedimentary basin fill and underlying rocks, the seafloor in marine areas, and other significant intra-crustal interfaces. The relationships between flexural deformation and the gravity field are dependent upon the partitioning of the load between surface (topography) and subsurface (buried density anomalies) components (Karner and Watts, 1983; Forsyth, 1985; Macario et al., 1995; McKenzie, 2003).

The flexural behavior of the Earth’s lithosphere can be approximated with a model in which the lithosphere is treated as a thin elastic plate overlying an inviscid fluid asthenosphere (Watts, 2001). Analytical solutions describing the flexural bending of homogeneous infinite and semi-infinite thin elastic plates subjected to point, line, or harmonic loads have been widely used in flexural isostatic analysis (Hertz, 1884; Nadai, 1963; Walcott, 1970). These elastic plate models have been shown to match the long-wavelength (>25 km) bathymetric features and the gravity field near most major oceanic tectonic features, including ocean trenches, ridges, island chains, transform faults, and hotspot swells (Turcotte, 1979; Watts et al., 1980; Bodine et al., 1981; Dahlen, 1981). On continents, elastic plate models have also been applied to analyze isostasy in continental rifts, mountain belts, and sedimentary basins (e.g., Karner and Watts, 1983; Watts et al., 1982; Chase and Wallace, 1988; Forsyth, 1985; McNutt et al., 1988; Flemings and Jordan, 1989; Weissel and Karner, 1989; Stern and ten Brink, 1989; ten Brink et al., 1997). Elastic plate models have also been used to assess glacial isostatic rebound (e.g., Stern et al., 2005) and to analyze isostasy on scales spanning continents and ocean basins (Banks et al., 1977; Watts and Burov, 2003; McKenzie and Bowin, 1976; Watts, 1978; McNutt and Menard, 1982; Lowry and Smith, 1994; Djomani et al., 1995; McKenzie and Fairhead, 1997; Petit and Ebinger, 2000; Flück et al., 2003; Crosby and McKenzie, 2009).

TAFI provides an easy-to-use interactive tool to model isostatic deformation of the lithosphere as flexural bending of an elastic plate. TAFI consists of a suite of MATLAB scripts and functions, coded as m-files, which facilitate rapid interactive forward modeling of flexural deformation of a thin elastic lithosphere overlying an inviscid asthenosphere. Only intrinsic MATLAB functions are used, providing portability across operating systems without need of special libraries or packages other than MATLAB itself. TAFI computes analytical Green’s functions that represent the flexural response of a thin elastic plate subjected to vertical impulsive (line or point) loads. The flexural response for non-impulsive loads can be obtained by convolving the Green’s functions with discretized functions representing loads that are distributed either in a radially symmetric pattern (simulating, for example, a seamount), along linear profiles (simulating loading by a thrust belt or seamount chain), or arbitrarily on a horizontal (x-y) plane (topography).

The Green’s function for a point load represents a two-dimensional (2-D) radial cross section in the r-z plane (r is the radial coordinate position) across a radially symmetric flexural basin, normalized by the magnitude of the load. In TAFI, these are referred to as axisymmetric plate models, with the basin shape described along a profile measured radially away from the point load. The point load Green’s function can be convolved with a spatially distributed load to compute the flexural response for axisymmetric loads of finite radius and for loads that are distributed arbitrarily in the horizontal (x and y) directions. In the first instance, the user supplies a load discretized as a function of radial distance, with TAFI returning a flexural profile across a similarly axisymmetric basin. In the second instance, the user supplies a load discretized on an x-y grid. In this case, TAFI returns a grid representing the three-dimensional (3-D; x, y, and z) shape of the basin. Accordingly, TAFI refers to such models as 3-D plate models.

The Green’s function for a line load is a 2-D flexural profile (in the x-z plane) across a basin of infinite length, normalized by the magnitude of the load. In TAFI these are referred to as 2-D plate models, with the basin shape described along a profile measured perpendicular to the trend of the line load. TAFI computes the flexural deformation for spatially distributed 2-D loads (loads with finite width and infinite along-strike extent) by convolving the line load Green’s function with a load discretized along a profile trending perpendicular to strike (the x-direction). Greens functions for both infinite (unbroken) plates and semi-infinite (broken) plates are provided in TAFI. For 2-D models, TAFI also provides an analytical solution for harmonic loads, in which the load magnitude varies sinusoidally in the x-z plane and is of infinite extent in the y-direction.

TAFI computes the gravity field resulting from the flexed plate using Parker’s method (Parker, 1973). This allows the user to specify a layered density structure for the lithosphere, with the gravitational attraction of each flexed layer summed to compute the total gravity anomaly. Examples of density interfaces that contribute strongly to the gravity anomaly are the crust-mantle boundary, the interface between the sedimentary basin fill and underlying crust, and in marine environments, the seafloor.

The MATLAB geodynamic modeling functions provided in TAFI to calculate plate flexure and the corresponding gravity anomaly are similar to those available in other software packages. Examples include Lithflex1 and Lithflex2 (http://csdms.colorado.edu/wiki/Geodynamic_models; Slingerland et al., 1994); Flex2D (http://www.ux.uis.no/∼nestor/work/programs.html); MECAIR (http://aconcagua.geol.usu.edu/∼arlowry/code_release.html); gFlex, (https://github.com/awickert/gFlex; Wickert, 2016), and MATLAB routines by Fredrik Simons (http://geoweb.princeton.edu/people/simons/software.html). The primary feature that distinguishes TAFI from these other software packages is TAFI’s easy-to-use graphical user interface (GUI), which facilitates interactive flexural modeling through easy modification of model parameters and instantaneous visual comparison of model results to data. The GUI provides a powerful instructional tool, allowing users to interactively explore how changes in plate boundary conditions, plate strength, load position and magnitude, and superposition of loads affect flexural deformation. TAFI’s GUI is designed to serve as a template to which additional isostatic analysis tools, such as solutions for different plate rheologies (e.g., viscous, viscoelastic, or visco-plastic), may be added by users familiar with MATLAB coding.

The maximum deflection of the plate (w_max), amplitude of the peripheral uplift (w_b), distance from the load to the crest of the peripheral uplift (x_b), and the point of zero deflection (x₀) between the basin and peripheral uplift are attributes commonly used to compare flexural models to observations (Fig. 1). These attributes, which describe the shape and width of the basin, depend on the flexural parameter, α, and are obtained in closed form for the line and harmonic load models (Table 1). For the point load model, these values are determined by searching the calculated deflection curve for the minimum, maximum, and zero value nearest to the load.

Here, F(Δg) is the Fourier transform of the gravity anomaly, r is the position vector, k is the wavenumber, G is the universal gravitational constant, Δρ is the density contrast across the interface, F(h) is the Fourier transform of the de-trended flexural relief on the interface (prescribed by the plate flexure), and z₀ is the mean depth of the interface. Inverse Fourier transforming F(Δg) provides the gravity field produced by topography on the buried interface. The gravitational effect of multiple density interfaces within the lithosphere is obtained by summing the effect of each interface.

TAFI is structured hierarchically, with an upper-level script that provides a GUI to lower-level functions that perform the geodynamic calculations, visualization, and input-output operations (Fig. 2; Table 2). The appearance of the GUI is defined in the file “TAFI.fig” (Fig. 3 ¹). This is a MATLAB encoded file that defines the elements available from within the GUI. TAFI.fig also associates the GUI elements with actions and variables used in TAFI for modeling and visualization. The actions invoked by user interaction with the GUI elements are prescribed in the elements’ respective callback functions, which are coded in the file “TAFI.m”. The callback functions call TAFI lower-level functions that perform computational, data input-output, or visualization tasks and pass variable values obtained from the GUI elements (e.g., the flexural rigidity) to these functions.

The geodynamic and gravity modeling functions are MATLAB implementations of Equations 4–7 and Equation 8, respectively. These functions can also be called from the command line, allowing access from outside of the GUI or from within user-written programs (Table 3). The MATLAB implementation of Equations 4–6 to obtain Green’s functions (by setting Q₀ = 1) for the flexural response of an elastic plate to impulse (line or point) loads is straightforward, as is the implementation of Equation 7 to describe harmonic loading. Flexure due to non-unit impulse loads is obtained by scaling the Green’s function by a user-specified load magnitude. For distributed 2-D loads, the flexural response is calculated by convolving the user-supplied discretized load function with the semi-infinite or continuous plate Green’s functions for line loads (Equations 4–5). For 3-D distributed loads, the flexural response is calculated by convolving the user-specified load matrix with the point load Green’s function (Equation 6). MATLAB’s intrinsic convolution functions require that the load and Green’s functions be discretized at similar uniform intervals. For 2-D models, the load is resampled to the Green’s function interval (specified by the user) prior to convolving. For 3-D loads, the Green’s function is calculated at the grid spacing used in the user-supplied load function, and no resampling is required.

To compute the gravity anomaly, TAFI requires the Fourier transform of the deflection in the wavenumber domain (Equation 8). Estimation of the gravity field in the wavenumber domain via the series summation in Equation 8 is then straightforward. The gravity field in the wavenumber domain is then inverse Fourier transformed to yield the gravity field as a function of position. Four terms of the series in Equation 8 are used, as empirical testing showed that this is sufficient to reproduce the calculated gravity anomaly to within 10^–9 mGal in realistic modeling scenarios.

TAFI’s installation, directory structure, use and properties of GUI elements, default behaviors, input-output file formats, and other operating instructions are described in Supplemental File 1². The TAFI m-files are in Supplemental File 2³. Here, we briefly describe TAFI’s functionality. All operations described below can be done outside of TAFI’s GUI by invoking TAFI’s functions from the command line as described in Supplemental File 1 [^{footnote 2}]. The GUI is designed to allow the user to vary model parameters and observe the results immediately using sliders, edit boxes, and dropdown menus, thus facilitating interactive fitting of flexural models to data.

Flexure and gravity modeling in TAFI requires passing model parameters as arguments to TAFI’s geodynamic modeling functions, which return calculated deflection and gravity fields. The model parameters are specified in the function invocation when using TAFI from the MATLAB command line. When using the TAFI GUI, the model parameters are specified with sliders, edit boxes, and dropdown menus (Fig. 3 [see ^{footnote 1}]). The flexure and gravity curves are displayed in the GUI plot window and are updated interactively when the model parameters are changed.

To build a flexure model in TAFI, the user must specify the plate type, load geometry, flexural rigidity, load magnitude, load position, and density contrast between the mantle and basin fill. The functional form of the flexural response of the elastic plate for 2-D and 3-D flexure models is prescribed by the plate geometry (infinite or semi-infinite plate) and the load geometry (selecting from line, point, 2-D harmonic, 2-D distributed, axisymmetric distributed, or 3-D distributed load geometries; Table 4). The remaining model parameters are defined by the user and passed to TAFI’s geodynamic functions. In the case of distributed loads (2-D or 3-D), load magnitudes are provided by the user as a MATLAB column vector (2-D) or matrix (3-D). For sinusoidal loads, the load wavelength (λ) instead of the load magnitude is passed to the corresponding geodynamic function. In the GUI, the user can vary these variables with sliders and edit boxes and see the change in the computed flexural curve displayed in real time. Using the TAFI GUI, the user can also interactively change the load position and scale the magnitude of distributed loads. TAFI uses default values for g (9.8 m/s²), Earth’s radius (R, 6370 km), E (8.5 × 10¹⁰ N/m²), and Poisson’s ratio (0.25). These values can be modified by users as described in Supplemental File 1 (^{footnote 2}).

The elastic plate is required to have at least two density interfaces, one at the top of the plate and another at the base of the plate. These density interfaces represent the boundaries between the basin fill and underlying crust, and between the crust and mantle (Fig. 1). The mantle and infill densities are required to calculate the flexural parameter (α) (Table 1), which is used to compute the flexural deflection of the lithosphere for non-harmonic loads (Equations 4–6). The density contrasts between these layers and the depth of interfaces are also used to calculate the gravity anomaly associated with flexure of the lithosphere (Equation 8). When using the TAFI GUI, the density interfaces are defined by specifying the basin fill, crust, and mantle densities. When invoking TAFI from the MATLAB command line, the interfaces are defined by specifying the density contrasts between the layers and their depths. Additional density interfaces, representing lithological boundaries within the lithosphere, can be specified when accessing TAFI from the command line. These additional density interfaces are used to compute the gravity anomaly, but are not used in the flexural deformation calculation.

TAFI can be used to explore plate flexural behavior without loading data. However, the primary purpose of TAFI is to facilitate forward modeling to determine the flexural parameters that best fit the observed vertical deflection and gravity data. To demonstrate the TAFI workflow, we model the flexural deformation of the lithosphere along a 2-D profile crossing the Aleutian Trench (Fig. 4). The Aleutian Trench extends ∼2900 km from the Gulf of Alaska (USA) to the Kamchatka Peninsula (Russia) and is the plate boundary where the Pacific plate is being subducted under the North American plate. A number of studies have shown that the seafloor bathymetry near the Aleutian Trench can be fit with a 2-D semi-infinite elastic plate model subjected to a line load (Caldwell et al., 1976; Levitt and Sandwell, 1995; Watts, 2001; Bry and White, 2007).

We chose a profile that crosses the trench ∼45 km to the west of the “Aleutian 8” profile of Levitt and Sandwell (1995) (Fig. 4). Bathymetric and free-air gravity data along this profile were downloaded from the U.S. National Geophysical Data Center (NGDC). We prepared the data for use in TAFI by removing the mean and regional tilt to correct for elevation changes not associated with flexure (primarily, changes in seafloor depth due to plate age). The corrected flexural and gravity profiles were saved in text (.txt) format (format is described in Supplemental File 1 [^{footnote 2}]).

Launching TAFI from the GUI shows the default infinite plate model. Dropdown menus are used to change the plate geometry and load shape from the default settings to semi-infinite and 2-D line load respectively. Following Levitt and Sandwell (1995), the TAFI model is then modified to change the infill, crust, and mantle densities and the depths to their associated interfaces to 1000 kg/m³, 2750 kg/m³, 3200 kg/m³, 7000 m, and 12,000 m respectively. TAFI’s plot window displays both the initial and modified flexure and gravity curves.

To constrain the flexure and gravity modeling of the Aleutian Trench by curve fitting in TAFI, the Aleutian bathymetry and gravity data (saved earlier) are imported into TAFI (Fig. 5) (data import process and file format are described in Supplemental File 1 [^{footnote 2}]). Curve fitting involves adjusting the modeled flexure and gravity curves to fit the imported data. This is done by revising the flexural rigidity, load magnitude, and load position through their respective sliders and edit boxes in TAFI’s GUI. The model curves are updated each time a parameter value is changed, with the new model curves added to the existing plot. The curve-fitting process is terminated when the modeled curves fit the data reasonably well (Fig. 6). The final flexural model parameters and attributes for the best-fitting Aleutian Trench TAFI model (Fig. 6) are found to lie within the range of values estimated by Levitt and Sandwell (1995) and Bry and White (2007) (Table 5).

The flexural responses for 2-D and 3-D line and point loads (Equations 4–6) and 2-D sinusoidal loads (Equation 7) were verified by comparing the results obtained from TAFI with solutions computed in Microsoft Excel and Generic Mapping Tools (Wessel et al., 2013), and by comparing the zero crossing, maximum deflection, and flexural bulge position and amplitude with analytical solutions where available (Table 1). Gravity solutions for the line, point, and harmonic loads were verified by comparing TAFI solutions with gravity models calculated from the flexural admittance function (Watts, 2001). TAFI’s solutions for distributed 3-D loads were verified by computing 3-D models with (1) a point load located in the center of the load grid, and (2) a line load crossing the load grid. Radial deflection and gravity profiles extending away from the point load and striking perpendicular to the line load were verified to match with TAFI’s comparable axisymmetric point load and 2-D line load solutions, respectively (Fig. 7).

The Toolbox for Analysis of Flexural Isostasy (TAFI) consists of MATLAB functions that compute the flexural response of an infinite or semi-infinite thin elastic plate overlying an inviscid substrate and subjected to line, point, and 2-D periodic loads. TAFI also provides functions to convolve the flexural Green’s functions with non-impulsive loads, allowing computation of deflection under spatially distributed 2-D or 3-D loads, and provides functions to compute the Bouguer or free-air gravity fields associated with plate flexure. The geodynamic functions in TAFI can be called from the MATLAB user interface or within user-written MATLAB code, or accessed through TAFI’s graphical user interface (GUI). The GUI facilitates interactive flexural modeling and comparison of the model to gravity observations and data constraining flexural subsidence and uplift. TAFI uses only MATLAB intrinsic functions and functions distributed as m-files with TAFI, and does not require additional toolboxes or libraries.

We would like to thank geophysics graduate and undergraduate students of the Department of Geosciences, Warner College of Natural Resources, Colorado State University, who reviewed the toolbox and provided constructive comments. We thank Dr. Andrew Wickert, University of Minnesota, whose reviews made this paper better. We also thank the science editor of Geosphere, Dr. Raymond M. Russo, for his helpful suggestions. This material is based upon work supported by the National Science Foundation under grants 1043700 and 1358664, and under Cooperative Agreement 0342484 through subawards administered and issued by the ANDRILL Science Management Office at the University of Nebraska–Lincoln, as part of the ANDRILL U.S. Science Support Program.