# Seismological Research Letters

- © 2016 by the Seismological Society of America

## ABSTRACT

GeoTess is a model parameterization and software support library that manages the construction, population, storage, and interrogation of data stored in 2D and 3D Earth models. The software is available in Java and C++, with a C interface to the C++ library. The software has been tested on Linux, Mac, Sun, and PC platforms. It is open source and is available online (see Data and Resources).

## INTRODUCTION

The use of multidimensional global Earth models in both the seismological research community and in the operational systems of global monitoring agencies has become increasingly ubiquitous in recent years and is destined to become even more prevalent in the future. This trend is driven by the significantly enhanced fidelity of these models as compared to 1D models, resulting in the promise of more accurate predictions of seismic observations and hence to more accurate and precise estimates of seismic event locations and magnitudes. Use of multidimensional global models has become feasible due to the dramatic increase in recent years of the computational capabilities of readily available commercial off‐the‐shelf computational resources.

Efforts to evaluate multidimensional models and incorporate them into operational systems are complicated by the unique formats and software support systems that are implemented by each development team (Reiter *et al.*, 2011). This makes evaluation and implementation of 3D models by monitoring agencies difficult and inefficient and reduces the probability that research products developed by the seismological research community will ever find their way into operational systems.

GeoTess is an attempt to address this issue by providing a common model format designed to facilitate the implementation and sharing of 3D Earth models. GeoTess is composed of a model parameterization for multidimensional Earth models and extensible software that implements the construction, population, storage, and interrogation of data stored in the model. GeoTess is not limited to any particular type of data; in GeoTess, the data are just a collection of values associated with each node in the grid.

GeoTess uses triangular tessellations to represent the geographic geometry and topology of the model grid. Although software to support a parameterization based on a regular latitude–longitude grid would be much more straightforward to develop, such grids suffer from severe unintended variability in cell areas, with cell areas approaching zero near the poles (Fig. 1). Software that supports triangular tessellations, on the other hand, is somewhat more complicated to develop but results in grids with much more uniform cell size and ∼25% fewer vertices. Triangular tessellations have been used in a number of Earth science applications (Wang and Dahlen, 1995; Wang *et al.*, 1998; Simmons *et al.*, 2011; Ballard *et al.*, 2009; Myers *et al.*, 2010; Pasyanos *et al.*, 2014).

## CAPABILITIES

GeoTess software provides tools to perform the following functions:

Construct the discretized grid of a model by defining the positions of all the nodes (geometry) and the connectivity of the nodes with one another (topology or connectivity) to form nonoverlapping triangles that completely span the unit sphere. Construction of grids with variable resolution in both geographic and radial dimensions is supported.

Query a model grid for information about the nodes, cells, or tessellations.

Associate data structures with the nodes.

Query the model for the data associated with a specified node.

Modify the data associated with a node.

Rapidly find arbitrary positions within the grid hierarchy (point searching).

Given an arbitrary position within the model space that does not correspond to a node of the model geometry, calculate a set of interpolation coefficients that can be applied to data at nearby nodes of the model geometry to obtain an estimate of the data value at the point of interest. GeoTess currently implements linear and natural‐neighbor (Sibson, 1981) interpolation algorithms and may implement additional methods in the future.

Calculate weighting coefficients for a specified ray path through a model (useful in tomographic inversion applications).

Read models from, and write them to, files in custom ASCII and binary formats.

Grid construction is implemented with an application called GeoTessBuilder, which is written in Java and can be run on any computer platform that supports that language. The remaining functions described above are implemented in a software library called GeoTess, which is available in Java and C++, with a C interface to the C++ library. The libraries have been compiled and tested on SunOS, Linux, Mac OSX, and Windows operating systems and can be downloaded from the Sandia National Laboratories website (see Data and Resources). The website also has links to several example codes that illustrate the use of some of the most common GeoTess functions.

## MODEL COMPONENTS

A GeoTess model is comprised of the following elements.

*Layers:*A set of layers (Fig. 2). Each layer spans the entire 2D geographic extent of the model. The boundaries at the top and bottom of a layer may have topography. Within each layer, model data values are continuous, both geographically and radially. Model data values may be discontinuous across layer boundaries. Layers may have zero thickness at some or all geographic locations. An important limitation of the parameterization used by GeoTess is that layer boundaries may not fold back on themselves; that is, any radial line emanating from the center of the Earth must intersect each layer boundary exactly one time.*Tesselations:*A set of multilevel tessellations (Fig. 3). Each layer will be associated with one multilevel tessellation, but many layers may be associated with the same multilevel tessellation; that is, there is a many‐to‐one relationship between layers and multilevel tessellations. By associating layers that are deep in the Earth with low‐resolution multilevel tessellations and layers at shallower levels in the Earth with higher‐resolution multilevel tessellations, the resolution of the model can be varied radially as necessary to achieve more appropriate sampling (Fig. 2).*Topology:*The topology of each multilevel tessellation will consist of a set of levels (Fig. 3), with each level consisting of a set of triangles that spans the surface of a unit sphere, without gaps or overlaps. The triangles on a given tessellation level are obtained by subdivision of the triangles on the previous tessellation level, with the first tessellation level being either a tetrahedron, octahedron, icosahedron, or a tetrahexahedron. Each multilevel tessellation may have variable resolution in the geographic dimensions by only subdividing a subset of the triangles on the previous tessellation level.*Geometry:*The geometry of each multilevel tessellation will consist of a set of vertices that defines the positions of the corners of the triangles. If a model is comprised of more than one multilevel tessellation, they will share common vertices, to the extent possible.*Data arrays:*Each data array is a 1D array of data values that may be of type double, float, long, int, short, or byte. All the data arrays in the model must be of the same type and must have the same number of elements.*Profiles:*Each profile is composed of a list of monotonically increasing radii and a list of data arrays. Each profile is associated with a single vertex and a single layer in the model. The first and last radii in a profile define the bottom and top of the associated layer at the geographic position of the vertex. Several types of profiles are supported:*N*‐Point profiles consist of two or more radii and an equal number of data arrays, with one data array associated with each radius.Constant value profiles consist of two radii and a single data array that defines the data values for the entire radial span of the profile.

Thin profiles consist of a single radius and a single data array. They have zero thickness; that is, the radius of the bottom and top of the profile are equal.

Empty profiles consist of two radii but no data arrays.

Surface profiles consist of only a single data array. They have no radius values. These are used to support 2D models.

Empty surface profiles contain no information at all. These are used to support 2D models.

The data values within a profile are continuous.

7.

*A 2D array of Profiles with nVertices × nLayers elements:*The first index refers to one of the vertices of the model geometry, and the second index refers to one of the layers of the model. For a given vertex index, the 1D array of profiles contains a profile for each layer of the model, stored in the order of increasing radius. The last radius of each profile in a 1D profile array must be equal to the first radius of the next profile in the same 1D profile array. Although the data values within a single profile are continuous, data values may be discontinuous across the profile (i.e., layer) boundaries.8.

*Radial interpolators that interpolate data values in profiles:*These include linear interpolators, cubic spline interpolators, and potentially others.9.

*2D interpolators that interpolate values in the two geographic dimensions:*These include a linear interpolator that interpolates values within a single triangle of the 2D tessellations and the natural neighbor interpolator that provides continuous spatial derivatives of the data values (Sibson, 1981).10.

*(2+1)D interpolators that combine 1D and 2D interpolators to interpolate data in 3D:*They first use a 1D interpolator to interpolate values at a specified radius in a neighborhood of profile arrays, and then apply a 2D interpolator to those values to find an interpolated value at the desired 3D location.

## GRID GENERATION

Grid generation is accomplished with a software application called GeoTessBuilder, which operates in two fundamental modes. In the first, GeoTessBuilder will construct a new 2D geographic grid from scratch. The user first specifies a default grid resolution (the edge length of the triangles, in degrees). If nothing further is specified, then a uniform tessellation with that geographic resolution will be constructed. The user may also specify a collection of points, paths, and polygons, or any combination thereof. GeoTessBuilder will refine the grid in and around those points, paths, and polygons, to whatever grid resolution the user requests (Fig. 4). The points, paths, and polygons can be specified in ASCII text files or can be accessed from KML files generated by Google Earth. Specification of the radii and data values along the profiles of the model is performed separately by user‐generated software (examples are provided with the code).

In the second mode of operation, GeoTessBuilder will construct a new GeoTess model by refining an existing model. The user supplies an existing model and a list of grid node indexes. GeoTessBuilder will add additional grid nodes around the specified grid nodes, effectively doubling the spatial resolution around the supplied model. Refinement is performed in both the geographic and radial dimensions. Alternatively, the user can specify a threshold value of one of the data attributes stored in the model. GeoTessbuilder will add additional nodes around existing nodes where the value of the specified data attribute exceed the specified threshold value.

## SALSA3D MODEL

We offer two examples of models based on the GeoTess parameterization. The first is SALSA3D, a global *P*‐velocity model of the Earth mantle (S. Ballard *et al.*, unpublished manuscript, 2016; see Data and Resources). The SALSA3D model is deployed on a grid that was generated with an adaptive grid refinement algorithm. The starting model consisted of a triangular tessellation with uniform 8° edge lengths (Fig. 5a). After an iterative tomographic inversion procedure had converged, the model resolution matrix was computed (Fig. 5b). All grid nodes where the diagonal of the model resolution matrix was greater than 0.25 were refined to effectively double the spatial resolution of the model in the vicinity of the nodes with high model resolution. This procedure (tomography followed by the model resolution calculation and adaptive grid refinement) was repeated four times, at the end of which the grid had spatial resolution as low as 1° (Fig. 5c) and the model resolution was everywhere less than 0.25 (Fig. 5d). This procedure resulted in a model with approximately an order of magnitude fewer nodes than would have been the case if the spatial resolution had been 1° everywhere.

## LibCorr3D

LibCorr3D is a software library that extends GeoTess and is used to manage station‐phase specific predictions of seismic travel time and the uncertainties of those predictions. For every station in a network and every phase of interest, a separate GeoTess model is generated that stores two quantities: the predicted travel time and the uncertainty of the prediction. Seismic event location software can use these models by looking up predictions and model uncertainties, rather than computing them on the fly. An alternative strategy would be that LibCorr3D would store corrections to values computed with a 1D model, such as ak135 (Kennett *et al.*, 1995), rather than the predictions themselves, and the locator would add the corrections to the values it calculated with the 1D model.

Figure 6 illustrates the contents of a typical LibCorr3D model for station GERES, located in Germany. The difference in travel time computed using SALSA3D and ak135 is illustrated in the Figure 6a. Figure 6b illustrates the uncertainty of the travel‐time predictions calculated using SALSA3D (Ballard *et al.*, unpublished manuscript, 2016). Travel‐time differences and uncertainties extend out to 95° from the station, because that is the distance at which seismic rays begin to interact with the core–mantle boundary. Although the surface values of predicted travel time and uncertainty have been emphasized, it is important to note that these values are stored in a 3D grid. Figure 6c is the depth of the bottom of the LibCorr3D models. Over most of the globe, the model is only one node thick (blue in Fig. 6c), but the model extends deeper than any recorded earthquake where warranted by recorded seismicity. Standard 1D distance‐dependent travel‐time uncertainty is shown in Figure 6d for comparison with the path‐dependent travel‐time uncertainties computed with SALSA3D.

This strategy has implications for computational performance. A typical model for a single station phase, stored at a resolution of 0.5° horizontally and 30 km radially, will require roughly 50 MB of memory once it is loaded into computer memory. For a network of ∼60 stations, a single phase, and two attributes (travel time and travel‐time uncertainty), this will require ∼6 GB of RAM to store, which is a modest amount of memory on a modern computer.

## SUMMARY

GeoTess is a model parameterization and software support library that manages the construction, population, storage, and interrogation of data stored in 3D Earth models. It is our hope that nuclear‐explosion monitoring agencies will familiarize themselves with GeoTess and implement the ability to accept 3D models in GeoTess format. If this happens, and members of the research community develop the ability to deliver data products in this format, it would greatly facilitate the sharing of 3D Earth model information between model developers and monitoring agencies and would likely increase the probability that research products are properly evaluated and ultimately implemented by monitoring agencies.

## DATA AND RESOURCES

The GeoTess software and examples are available at www.sandia.gov/geotess (last accessed September 2015). The unpublished manuscript “SALSA3D–A tomographic model of compressional wave slowness in the Earth’s mantle for improved travel time prediction and travel time prediction uncertainty” by S. Ballard, J. R. Hipp, M. L. Begnaud, C. J. Young, A. V. Encarnacao, E. P. Chael, and W. S. Phillips submitted to BSSA.

## ACKNOWLEDGMENTS

We thank Scott Phillips and David Yang for trying out GeoTess in their research. We are also grateful to David Yang, Megan Slinkard, Stephen Arrowsmith, and an anomymous reviewer for reviewing this article. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the U.S. Department of Energy’s National Nuclear Security Administration under Contract Number DE‐AC04‐94AL85000.