The Antarctic Crust and Upper Mantle: A Flexible 3D Model and Software Framework for Interdisciplinary Research

Stål, Tobias; Reading, Anya M.; Halpin, Jacqueline A.; Phipps, Steven J.; Whittaker, Joanne M.

doi:10.3389/feart.2020.577502

TECHNOLOGY AND CODE article

Front. Earth Sci., 27 November 2020

Sec. Solid Earth Geophysics

Volume 8 - 2020 | https://doi.org/10.3389/feart.2020.577502

The Antarctic Crust and Upper Mantle: A Flexible 3D Model and Software Framework for Interdisciplinary Research

Tobias Stål^1,2^*

Anya M. Reading^1,2

Jacqueline A. Halpin²

Steven J. Phipps²

Joanne M. Whittaker²

¹School of Natural Sciences (Earth Sciences/Physics), University of Tasmania, Hobart, TAS, Australia
²Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, TAS, Australia

Interdisciplinary research concerning solid Earth–cryosphere interaction and feedbacks requires a working model of the Antarctic crust and upper mantle. Active areas of interest include the effect of the heterogeneous Earth structure on glacial isostatic adjustment, the distribution of geothermal heat, and the history of erosion and deposition. In response to this research need, we construct an adaptable and updatable 3D grid model in a software framework to contain and process solid Earth data. The computational framework, based on an open source software package agrid, allows different data sources to be combined and jointly analyzed. The grid model is populated with crustal properties from geological observations and geochronology results, where such data exist, and published segmentation from geophysical data in the interior where direct observations are absent. The grid also contains 3D geophysical data such as wave speed and derived temperature from seismic tomographic models, and 2D datasets such as gravity anomalies, surface elevation, subglacial temperature, and ice sheet boundaries. We demonstrate the usage of the framework by computing new estimates of subglacial steady-state heat flow in a continental scale model for east Antarctica and a regional scale model for the Wilkes Basin in Victoria Land. We hope that the 3D model and framework will be used widely across the solid Earth and cryosphere research communities.

1. Introduction

Past, present, and future changes in the mass of the Antarctic ice sheets have a direct impact on global sea level (e.g., King et al., 2012; Shepherd et al., 2012; Golledge et al., 2015; Ritz et al., 2015; DeConto and Pollard, 2016; Golledge et al., 2019). During the 21st century and beyond, the projected rise in sea level in response to anthropogenic climate change is expected to have enormous social and economic consequences (e.g., Kulp and Strauss, 2019; Oppenheimer et al., 2019). Constraining the likely response of ice sheets to global climate change is therefore a high priority. The mechanisms controlling the extent and thickness of the cryosphere involve interaction with the atmosphere (e.g., Frieler et al., 2015; DeConto and Pollard, 2016; Lenaerts et al., 2016), the ocean (e.g., DeConto and Pollard, 2016; Dinniman et al., 2016; Rintoul et al., 2016), and the crust and mantle beneath, which is the focus of this contribution. Examples of solid Earth–cryosphere interaction include the impact of the heterogeneous Earth structure on glacial isostatic adjustment (e.g., Whitehouse, 2018), the amount and distribution of geothermal heat (e.g., Pattyn, 2010), and the history of erosion and deposition over geological time (e.g., Paxman et al., 2018). The continental crust is a highly heterogeneous layer usually characterized by a combination of geological observations, geochronological results, tectonic plate reconstructions, and geophysical surveys to obtain an overall picture of the composition, age, evolution, and 3D architecture of its constituent units. A sharp change in seismic wave speed, the Mohorovičić discontinuity (Moho), defines the boundary between the crust and the mantle beneath (Christensen, 1988; An et al., 2015a). The upper mantle provides a rigid and tectonically mobile component, which together with the crust forms the continental lithosphere. A deeper seismic discontinuity, the lithosphere–asthenosphere boundary (LAB), indicates the transition to a ductile mantle as a result of increasing temperature and pressure with depth (Artemieva, 2011). Many aspects of the Earth’s crust and mantle have significant spatial variability that impacts overlying ice sheets; hence, access to solid Earth research results has gained importance to the interdisciplinary research community (Whitehouse et al., 2019).

1.1. Geology, Geochronology, and Geochemistry

Our understanding of the Antarctic crust is restricted by the ice cover that leaves only 0.18% of the rocks exposed (Burton-Johnson et al., 2016), with access further limited by logistical difficulties. Early field campaigns enabled geological investigations to map out crustal domains along the Antarctic coast and Transantarctic Mountains (Ravich et al., 1965; Craddock, 1970; Adie and Adie, 1977; Tingey et al., 1991). Those interpretations are, to a large extent, still valid, although more recent field geological studies have expanded the number of outcrops visited. Geochronology and geochemistry have added insights to refine our understanding by constraining event chronologies, derive likely tectonic environments, and, in conjunction with geophysics, also allows geological correlation (regional and local studies include, e.g., (geographically) clockwise around the Antarctic continent: Halpin et al., 2005, 2012; Corvino et al., 2008; Williams et al., 2018; Daczko et al., 2018; Tucker et al., 2017; Morrissey et al., 2017; Maritati et al., 2019; Di Vincenzo et al., 2007; Goodge et al., 1992; Siddoway et al., 2004; Yakymchuk et al., 2015; Burton-Johnson and Riley, 2015; Will et al., 2009; Jacobs et al., 1998; Marschall et al., 2010).

Interpretations of Antarctic geology are often contextualized in a tectonic reconstruction framework (Du Toit, 1937; Whittaker et al., 2013b; Matthews et al., 2016; Williams et al., 2019) and can hence be guided by data from continents that were adjoined in Gondwana, especially Australia, India, and Africa (e.g., Yoshida et al., 1992; Fitzsimons, 2000; Aitken et al., 2014; Daczko et al., 2018). Blocks of once continuous Archean cratons and orogenic belts are split between east Antarctica and Africa, India, and Australia. West Antarctica mostly consists of younger Phanerozoic crust (Siddoway, 2008; Boger, 2011; Artemieva and Thybo, 2020; Jordan et al., 2020). Archean and Paleoproterozoic crust is mainly cratonic, Proterozoic crust is formed by the reworked orogens of Nuna and Rodinia, and more recently, Phanerozoic crust has been added by Gondwanan and Cenozoic accretions and volcanism. Extensive reviews have drawn wellfounded interpretations for coastal regions (e.g., Boger, 2011; Harley et al., 2013; Jordan et al., 2020), but due to the lack of data, geological and tectonic maps of the ice covered interior rely significantly on extrapolation. An ongoing challenge is to access and incorporate the large amount of often inconsistent geological, geochronological, and geochemical studies. Initiatives such as the GeoMAP project (Cox et al., 2018) and compilations of rock sample data (e.g., Gard et al., 2019) aim to facilitate geological studies of Antarctica, using the broad range of published data.

1.2. Geophysics

Significant emphasis is placed on geophysical methods, particularly for East Antarctica, to infer geological information about ice-covered regions from remotely observed physical properties. Geophysical data are acquired from ground measurements, airborne instruments and satellites (Fowler, 1990).

Seismic measurements are sparse in Antarctica, and are often clustered according to the given regional study (e.g., Winberry and Anandakrishnan, 2004; Reading, 2006; Hansen et al., 2010; Hansen et al., 2016; Heeszel et al., 2016; Shen et al., 2018). Data from Antarctic deployments and global databases are used to generate continental scale seismic models (An et al., 2015a; Lloyd et al., 2020). Airborne geophysics coverage is variable across the continent. Large international campaigns, such as ICECAP (e.g., Roberts et al., 2011; Young et al., 2011; Aitken et al., 2014; Graham et al., 2017), acquire data over multiple summer seasons enabling extensive spatial coverage. Multiple datasets, including high resolution magnetic and gravity anomalies, surface elevation and ice penetrating radar are usually acquired simultaneously (e.g., Aitken et al., 2014; Robert et al., 2017) and Antarctic research has been accelerated by carefully curated compilations of such data (Fretwell et al., 2012; Morlighem et al., 2019). Notable regional airborne campaigns include targets such as Gamburtsev Subglacial Mountains (Ferraccioli et al., 2011), the South Pole satellite polar gap (Sneeuw and van Gelderen, 1997; Forsberg et al., 2017), Dronning Maud Land (Jacobs et al., 2015; Ruppel et al., 2018) and Transantarctic Mountains (Goodge and Finn, 2010). Magnetic data has been compiled as continental scale maps (ADMAP and ADMAP2, von Frese et al., 2007; Golynsky et al., 2013; Golynsky et al., 2017). Global satellite gravity surveys such as GOCE and GRACE are of particular importance in Antarctica due to the consistent cover of long-wavelengths anomalies (Visser, 1999; Pail et al., 2011; Förste et al., 2013). Continuous satellite measurements facilitate the identification of changes over time, such as mass loss (Velicogna, 2009; King et al., 2012), and changes in altimetry of the glacial surface from, e.g., CryoSAT-2 altimetry (Slater et al., 2018).

Modeling studies that are particularly important in the Antarctic context include making use of the curvature of gravity field (Ebbing et al., 2018), finding the elastic crustal thickness (Chen et al., 2018), comparison of models of, e.g., Moho depth from various approaches (Baranov et al., 2018; Pappa et al., 2019) and integrating density, compositional and thermal models (Haeger et al., 2019). Interpretation of magnetic anomalies combined with other datasets can support delineation of crustal domains (Goodge and Finn, 2010; Aitken et al., 2014; Ruppel et al., 2018; Paxman et al., 2019), and are also used to infer depth to the Curie temperature isotherm (Maule et al., 2005; Martos et al., 2017).

1.3. Solid Earth-Cryosphere Interactions

Mapping tectonic domains from geological data provides a first order segmentation of the lithosphere for 3D glacial isostatic adjustment models (Kaufmann and Wolf, 1999; Nield et al., 2018). Crustal heat production can to some extent be estimated from geochemistry (Hasterok and Webb, 2017) and geochronology (Jaupart and Mareschal, 2013). Likewise, mass transport by glacial exhumation and deposition is informed by geological and geochronological observations. From ground, airborne and satellite data, modeling exercises, and from comparisons with other continents, it is becoming increasingly apparent that we should expect large spatial variations in the subglacial physical properties of the crust and upper mantle in the Antarctic interior. This heterogeneity impacts solid Earth-cryosphere interaction on regional and local scales.

1.3.1. Glacial Isostatic Adjustment

Glacial isostatic adjustment (GIA) is the response of the viscous mantle and rigid lithosphere to changes in ice load (e.g., Whitehouse, 2018). As ice sheets melt, mass is transferred from the continent to the ocean, and the continental crust rebounds in response to the resulting buoyancy force. Lateral variations in lithospheric thickness and the viscosity of the deforming Earth’s mantle impact the rate and nature of this rebound (e.g., Kaufmann and Wolf, 1999; Nield et al., 2014; Nield et al., 2018). The crustal movement is measured by GPS time series (e.g., Martín-Español et al., 2016), and past uplift can be reconstructed from geomorphological observations by dating raised beaches, glacial erratics and sediments (White et al., 2010; MacKintosh et al., 2011). The observed elevation does not, in general, represent isostatic equilibrium as the Antarctic lithosphere is at present adjusting in response to changes in ice load and global sea level (Peltier, 2004; Whitehouse et al., 2012; Gunter et al., 2014; Whitehouse, 2018).

1.3.2. Subglacial Geothermal Heat

Geothermal heat flow, often termed ‘heat flux’ in ice sheet modeling studies, is a necessary boundary condition in many ice sheet models (e.g., Winkelmann et al., 2011). Heat at the base of slow flowing ice sheets can cause melting that impacts ice flow speed and can reduce the stability of the ice sheet. It can also affect the ice viscosity and hence affect internal deformation (e.g., Matsuoka et al., 2012; Petrunin et al., 2013; Pattyn et al., 2016). Heat is generated in the interior of the Earth and reaches the surface due to the temperature gradient. This is regulated by the thermal conductivity of the crust and mantle. Heat flow is known to be highly variable on continental, regional and local scales (Cull, 1982; Beardsmore and Cull, 2001; McLaren et al., 2003; Ramirez et al., 2016; Begeman et al., 2017; Jordan et al., 2018; Pollett et al., 2019). At plate margins and locations such as extensional basins, heat flow through convection or advection, by moving fluids and/or magma at depth, may be dominant.

Several different approaches are in current use to estimate the subglacial heat flow from modeled temperature gradients (Discussed by Lösing et al. (2020) and Burton-Johnson et al. (2020)). Magnetic derived heat flow maps are produced from either equivalent source magnetic dipole models (Maule et al., 2005) or magnetic spectral analysis from high resolution airborne data (Martos et al., 2017). Both methods are used to estimate a depth to the Curie temperature isotherm. Another approach uses seismic wave speed as an indirect measure of temperature at depth. Temperature is the main controlling factor of lateral variations in seismic wave speed in the upper mantle (Goes et al., 2000; Cammarano et al., 2003; Shapiro and Ritzwoller, 2004; An and Shi, 2007). An et al. (2015a) presented a surface wave tomography model constrained by receiver functions. From the wave speed, upper mantle temperatures are inferred and thermal gradients to the surface estimated (An et al., 2015b). Both the magnetic and seismic approaches have limitations due to their underlying assumptions, accuracy and resolution. A significant challenge when estimating subglacial heat flow is the need to account for the unconstrained lateral variations in heat production and thermal conductivity in the crust. Heat production varies over a large range for different rock types (Carson et al., 2014; Jaupart et al., 2016; Hasterok and Webb, 2017), and including geological knowledge in regional studies is of great value (e.g., McLaren et al., 2003; Burton-Johnson et al., 2017; Burton-Johnson et al., 2020). Direct measurements of the subglacial heat flow are very sparse in Antarctica (e.g., Fisher et al., 2015; Begeman et al., 2017), and some studies derive subglacial conditions from measurements within the ice (discussed by e.g., Mony et al., 2020). Heat anomalies are also known from radar images of the ice sheet (e.g., Schroeder et al., 2014; Jordan et al., 2018), the presence of subglacial lakes (Pattyn et al., 2016) and by inversion of ice sheet models (Pattyn, 2010).

1.3.3. Erosion and Deposition

The subglacial topography of Antarctica is the result of its tectonic evolution overprinted by cycles of erosion, exhumation and redeposition of sediment by rivers and glaciers. Topography can influence ice sheet dynamics through parameters such as direction of slope (e.g., Greenbaum et al., 2015), and fine-scale roughness (Goff et al., 2014; Graham et al., 2017). Subglacial topography is constrained by ice penetrating radar, gravity and seismic data. With data compilations such as Bedmap2 and BedMachine (Fretwell et al., 2012; Morlighem et al., 2019), a substantial part of the Antarctic subglacial landscape is revealed, but in many areas there are still large uncertainties (Fretwell et al., 2012; Graham et al., 2017). Glaciers are efficient in eroding and forming the landscape (Koppes and Montgomery, 2009; Cowton et al., 2012; Morlighem et al., 2019). Large amounts of sediment have been transported from Antarctica to the continental shelf and continental slopes (Whittaker et al., 2013a; Sauermilch et al., 2019), but in some areas the erosion has been very limited due to cold-based ice sheets that tend to preserve the existing topography (Jamieson and Sugden, 2008; Wilson et al., 2012; Paxman et al., 2018).

Understanding of the subglacial landscape evolution by erosion and deposition calls for an interdisciplinary approach, whereby ice sheet development, geophysical data and geological data are combined to constrain Antarctica’s past and present landscape, and isostasy (Jamieson and Sugden, 2008; Jamieson et al., 2010; Mackintosh et al., 2014; Paxman et al., 2016; Paxman et al., 2018; Paxman et al., 2019).

1.4. Motivation for the 3D Grid Model

Reproducible models of the Antarctic crust and upper mantle are needed to progress interdisciplinary studies such as those relating to GIA, heat flow and topography. A better understanding of the solid Earth is achieved by combining multiple data sources (Begg et al., 2009; Pappa et al., 2019; Stål et al., 2019). Populating models with current data presents a challenge, especially given the present rate of new data releases that have the potential to improve existing results. Lateral variations of crustal properties are often absent from large scale geophysical studies. One successful attempt to facilitate data access is the Quantarctica project that links data to users via a GIS application (Roth et al., 2017). Quantarctica allows users to directly visualise and compare datasets of a different nature. However, GIS might not be the first choice for multidimensional data processing, and a scripted framework is desirable for geophysical modeling and analysis.

In this contribution we present a flexible 3D grid model of the Antarctic crust and upper mantle. We populate the grid with datasets that have been used in univariate studies to constrain lithospheric rheology, heat flow and erosion and uplift: e.g., seismic wave speed, thermal properties, subglacial topography, geology and crustal segmentation models (Table 1). As a computational framework, we use agrid, an open software environment for storing, analysing and modeling multivariate and multidimensional data with functionality to visualize and export the results (Stål and Reading, 2020). agrid depends on well documented Python packages such as numpy (Harris et al., 2020), scipy (Jones et al., 2015), xarray (Hoyer and Hamman, 2017), dask (Rocklin, 2015), and rasterio (Gillies, 2013). Computations using numpy are as fast and memory efficient as compiled code (Van Der Walt et al., 2011), and chunk parallelization is made possible using dask arrays.

TABLE 1

TABLE 1. Datasets used to populate the grid, in alphabetic order.

The 3D grid model and computational framework are intended for a wide range of applications, and are designed to be updated as additional data become available. Thus, we make constraints and related uncertainty from geology, geochronology and geophysics available in a form that is usable by researchers in geoscience, glaciology and ice sheet modeling. Through this contribution, we aim to facilitate interdisciplinary studies on the interaction between the solid Earth and cryosphere of Antarctica.

2. Data

Our model and framework includes numerous geological and geophysical datasets, together with the source reference, as listed in Table 1. We limit the spatial extent of the grid to the present coastline and ice shelf grounding line (Mouginot et al., 2017). Some processing, such as resampling and interpolation, is applied when the data are imported. Data in global projections are first reprojected, then interpolated to avoid artifacts and distortion when interpolating across the South Pole and anti-meridian line. Some of the datasets included in this contribution certainly contain spatial distortion due to reprojection. This distortion typically has its origin when published results are stored to a global grid. We do not aim to correct those artifacts in this contribution, as this would modify the datasets and require further discussion. Instead, we include the datasets as they are published.

Uncertainty information relating to each parameter is included where available (E.g., Martos et al., 2017). Those provided uncertainty values might not capture the total range of uncertainty that arise from necessary assumptions and resolution. Refined analysis of datasets and uncertainty can be achieved in the framework. However, this is beyond this contribution.

All data are also associated with provenance information and metadata that links the original source. Metadata are stored with the dataset in the grid. The agrid package (Stål and Reading, 2020) contains methods to access the data directly from the original sources, open online repositories and through Quantarctica (Roth et al., 2017). Links to web addresses, current at the time of writing, are provided in the Supplementary Material. In the case that a link becomes outdated, error handling is provided. There is no limitation to the number of datasets that can be included in a model. The datasets listed here are included to produce the test cases for appraisal of the framework.

3. Methods and Results

In this section we outline the methods used to construct the 3D grid and illustrate the functionality of the computational framework through usage examples. All computations in this study are performed using the Python package agrid (Stål and Reading, 2020). Use of agrid facilitates easy programming and compact scripts, with the underlying software being tailored to computations that use data, and metadata, held in the 3D grid. The figures in this study are generated using only a few lines of high level code, and functions provided with agrid. Where applicable, we utilize perceptually linear color representation (Crameri and Shephard, 2019; Morse et al., 2019).

3.1. Populating the 3D Grid

To populate the 3D model, the datasets listed in Table 1 are imported. Datasets are re-sampled and interpolated to the defined extent, resolution, projection and cell sizes. Here we use bi-linear interpolation, but other refined techniques are available. Data imported from polygon vectors are rasterized and attributes saved to the grid using a map function. Observations at point locations, such as geochronological data (compiled by Gard et al., 2019), are binned to the containing grid cells. Datasets are projected to WGS 84/Antarctic Polar Stereographic (EPSG:3031), with very limited distortion in continental Antarctica. The total grid extent is set to 6,200 $\times$ 6,200 km with a horizontal resolution of 20 $\times$ 20 km (Figures 1–3). The extent and resolution of the grid can easily be modified and multiple resolutions can be used simultaneously. Using the same code, but with smaller extent and higher resolution, the Wilkes Subglacial Basin is shown as a grid with 2 $\times$ 2 km cells (Figures 4, 5C,D). The choice of values for depth sections can also be easily modified and is illustrated in Figure 1.

FIGURE 1

FIGURE 1. Oblique view of data held in the 3D grid model and illustration of plotting functionality. The model space is delineated by the Antarctic coastline and ice shelf grounding line (Mouginot et al., 2017). Depth sections are set to; 0, 8, 16, 25, 30, 35, 40, 45, 50, 75, 100, 150, 200, and 300 km. (A) Temperature in the crust and upper mantle derived from shear wavespeed by merging models AN1-Tc and AN1-Ts (An et al., 2015b) interpolated to defined grid. (B) Heat production in the crust from a simplified exponential function of depth, average production from age (Jaupart and Mareschal, 2013), segmentation by Schaeffer and Lebedev (2015) and crustal thickness from An et al. (2015a).

FIGURE 2

FIGURE 2. Examples of simple calculated outputs and visualisation in map view. Color representation is optimised for visibility. (A) Minimum subglacial heat flow from three studies (Maule et al., 2005; An et al., 2015b; Martos et al., 2017) using provided uncertainty ranges. (B) Maximum heat flow from the same three studies using provided uncertainty ranges. (C) Disagreement as standard deviation of the spread of the three studies. (D) Surface elevation with adjusted isostasy for ice removed. Calculated from Fretwell et al. (2012) and assuming constant density of ice 916.7 kg/m³, the crustal and mantle densities from Afonso et al. (2019). Moho from An et al. (2015a) and LAB from An et al. (2015b). A simple smoothing represents the rigidity of the lithosphere, as described in text.

FIGURE 3

FIGURE 3. New maps generated to show the methodology of using data held in the 3D grid model. (A) Segmentation from seismic tomography (Schaeffer and Lebedev, 2015). (B) Schematic geological age map (Tingey et al., 1991). (C) Actual geochronology compiled by Gard et al. (2019). The dataset is clipped by mapped rock outcrops from Burton-Johnson and Riley (2015) to mitigate errors. (D) Geological age estimated from a combination of the previous three datasets, with Gard et al. (2019) as preferred and indicated with shading in a strong tone, Tingey et al. (1991) as midtone, and Schaeffer and Lebedev (2015) in faint tone. Continental crustal age, and geochronological data are divided into three classes (Janse, 1984) and as discussed in text: Archean (purple), Proterozoic (green) and Phanerozoic (brown). Suggested oceanic crust in Schaeffer and Lebedev (2015) is shown in blue. White indicates no data (B,C).

FIGURE 4

FIGURE 4. New maps generated for Wilkes Basin showing data held in the 3D grid model and calculated outputs at higher resolution. (A) Present subglacial topography, Bedmap2 (Fretwell et al., 2012). (B) Subglacial topography, with ice removed and isostasy corrected using crustal thickness from An et al. (2015a) and lithospheric thickness from An et al. (2015b). Crustal and upper mantle densities from Afonso et al. (2019). Bedrock elevation and surface elevation from Bedmap2 (Fretwell et al., 2012). (C) Estimates of crustal age. Cells with geological observations in strong tone (Burton-Johnson and Riley, 2015; Gard et al., 2019), schematic geology (Tingey et al., 1991) in mid tone, and segmentation from Schaeffer and Lebedev (2015) in light tone. Continental crustal age is classified into three classes, Archean (purple), Proterozoic (green), Phanerozoic (brown), together with oceanic crust (blue). Methods are discussed in the text. (D) Crustal thickness from An et al. (2015a).

FIGURE 5

FIGURE 5. New maps generated by combining constraints from geophysical and geological data held in the 3D grid model: a new steady-state heat flow model, as discussed in text. (A) Heat flow map of East Antarctica, AqSS.ea (B) Uncertainty, as defined by the datasets used, excluding lateral uncertainties. (C) Heat flow map of Wilkes Basin, AqSS.wsb, (D) Uncertainty.

3.2. Computational Framework: Usage Examples

The agility of our 3D framework allows the rapid generation of maps or other outputs. Such products may be used to support research discussion or as numerical inputs for other studies (e.g., boundary conditions for ice sheet models).

3.2.1. Temperature in the Lithosphere and Heat Production in the Crust

Illustrating basic computation and oblique 3D visualisation using agrid and Antarctic datasets, Figure 1A shows lithospheric temperatures combined from AN-Ts and AN1-Tc (An et al., 2015b), interpolated to fit the grid. Figure 1B displays a first-order estimate of crustal heat production as a combination of crustal thickness (An et al., 2015a), segmentation (Schaeffer and Lebedev, 2015), heat production estimate from crustal age (Jaupart et al., 2016) and decreasing heat production as an exponential function of depth:

A = - A_{0} \times e^{- z / d_{Moho}} (1)

where A is the value of heat production in W/m³, $A_{0}$ is the average heat production, given the age of the crust, at that location, and $z / d_{Moho}$ is the fraction of depth to Moho, at the location.

3.2.2. Calculated Outputs Based on Multiple Geophysical Datasets

Illustrating further examples of computation and visualisation in map view, Figure 2 shows constraints from multiple heat flow models, and adjusted surface elevation based on multiple datasets. Minimum heat flow (Figure 2A) and maximum heat flow (Figure 2B) are the lowest and highest values at each grid cell in any of Maule et al. (2005), An et al. (2015b), and Martos et al. (2017), including provided uncertainty. Figure 2C shows the standard deviation as a measure of disagreement between the heat flow maps from aforementioned studies. Areas are readily seen where ice sheet modellers should be particularly careful when using the geothermal heat contribution as a boundary condition. The property maps shown in Figures 2A–C could therefore be useful for sensitivity studies of the impact of geothermal heat on the ice sheet at a continental scale.

Isostatic models are used to understand how the Antarctic crust and upper mantle interact with the cryosphere (e.g., O’Donnell and Nyblade, 2014). Figures 2D, 4B show bedrock elevation for isostatically relaxed ice-free conditions. Such computations are easy to perform in our framework, for example, using the simplified formula:

{DEM}_{iso} = {DEM}_{sg} + \frac{({DEM}_{s} - {DEM}_{sg}) \times ρ_{ice} \times D_{LAB}}{ρ_{crust} \times D_{Moho} + ρ_{mantle} \times (D_{LAB} - D_{Moho})} (2)

where ${DEM}_{iso}$ is the adjusted elevation model, ${DEM}_{sg}$ is the Bedmap2 subglacial elevation, ${DEM}_{s}$ is the surface elevation (Fretwell et al., 2012), $ρ_{ice}$ is the density of ice, assumed to be constant (916.7 kg/m³), and $ρ_{crust}$ and $ρ_{mantle}$ are applied from average crustal and lithospheric density in Afonso et al. (2019) reference model. We apply a 2D Gaussian kernel, with standard deviation of 60 km to include a simple constant model for the rigidity of the lithosphere. Figure 2D shows the elevation if the present ice mass were to be removed and the lithosphere regained its isostatic buoyancy. For ice sheet reconstructions of the past, or predictions of the future, the isostatic response of the solid Earth must be considered, as the coastline and ice shelf grounding lines are not static. Using our 3D model and framework, research tasks, such as testing alternative reconstructed ice masses, and recalculating the isostatic correction, are as straightforward as importing the modeled map of ice thickness.

3.2.3. Mapping Crustal Age by Merging Geological and Geophysical Datasets

Mapping crustal age provides an illustration of merging geological and geophysical sources, addressing the challenge of combining categorical and numerical data types. We utilize geochronological measurements compiled by Gard et al. (2019). The number of samples (Supplementary Material), mode, average value and standard deviation are calculated and binned to each cell. The legacy schematic geology map from Tingey et al. (1991) is used for reference and to guide moderate extrapolation of geology. Age estimates expressed in geological time are converted to age in years (Stål, 2020). Where no geological observations or extrapolation are available, we use crustal segmentation informed by seismic tomography. Most global regionalization studies often exclude or oversimplify Antarctica due to the limited available data (e.g., Jordan, 1981; Artemieva and Mooney, 2001; Artemieva, 2006; Artemieva, 2009). We implement one of the few continental scale segmentation models that covers Antarctica, the k-means clustering of surface-wave dispersion from Schaeffer and Lebedev (2015), which makes use of methods by Lekić et al. (2010) and data first presented by Schaeffer and Lebedev (2013). Examples of the standardised content reduced to three age classes and oceanic crust are shown on a continental (Figure 3) and regional scale (Figure 4C). The shading tone indicates the source, and hence, the robustness of the constraint. Direct observations (Gard et al., 2019) are strong in tone, schematic geological domains (Tingey et al., 1991) are shown in midtone and geophysical regionalisation (Schaeffer and Lebedev, 2015) is shaded in a faint tone. Combining data of different types is straightforward in concept, but challenging in practice, and the new framework shows that this can be achieved in a repeatable manner.

3.2.4. Calculated Outputs at Higher Resolution

Illustrating the functionality of the 3D model and framework at a regional scale, Figure 4 shows data held in the 3D grid and calculated outputs for the Wilkes Subglacial Basin. Figure 4A is a representation of the Bedmap2 dataset (Fretwell et al., 2012). Figure 4B shows the same simplified isostatic correction as Figure 2D in higher resolution. Figure 4C shows the combined model of crustal stabilisation age, using same methods as for Figure 3D, again at higher resolution, for the Wilkes Basin.

3.2.5. AqSS, a Steady-State Heat Flow Model

We further illustrate the functionality of the computational framework through generating a Steady-state heat flow model, AqSS, which combines geophysical and geological data. steady-state models can be reduced to two components that are identified as sources of geothermal heat: heat from the Earth’s core and mantle, reaching the crust as heat flow through the Moho, $q_{m}$ , and a commonly larger component, heat generated within the crust.

q_{g} = q_{m} + d_{m} \times A_{c} (3)

where $q_{g}$ is the subglacial heat flow, $q_{m}$ is the heat flow at the Moho, $d_{m}$ is the crustal thickness (Fretwell et al., 2012; An et al., 2015a) and $A_{c}$ is an average heat production within the crust.

From studies in different geological settings and methods, the mantle component has been constrained to $q_{m} = 14 \pm \sim 3 {mWm}^{- 2}$ (Guillou et al., 1994; Rudnick and Nyblade, 1999; Roy and Rao, 2003; Jaupart et al., 2016), which is within the uncertainty of most Antarctic heat flow estimates (e.g., Martos et al., 2017).

Uncertainty for AqSS is calculated from the uncertainty provided with each dataset, assuming they are independent.

σ_{q} = \sqrt{σ_{q m}^{2} + {(d_{m} \times A \sqrt{{\frac{σ_{A c}}{A}}^{2} + σ_{d_{m}}^{2}})}^{2}} (4)

where $σ_{q}$ is the absolute heat flow uncertainty, $σ_{q m}$ is the absolute uncertainty of heat flow into the crust, 3 mWm⁻² (reviewed by Jaupart et al., 2016). The relative uncertainty of crustal thickness ( $σ_{d_{m}}$ ) is set to 15%, A is the absolute mean heat production and $σ_{A c}$ is half of the range of heat production as suggested by Jaupart and Mareschal (2013) and listed in Table 2.

TABLE 2

TABLE 2. Thermal properties assigned to crustal domains.

By assuming steady-state conditions throughout East Antarctica and applying a constant contribution from the mantle (Mareschal and Jaupart, 2004), we avoid invoking any assumptions regarding temperatures in the lower crust or upper mantle. The larger part of the total heat flow is heterogeneous and originates from the crust (e.g., Jaupart et al., 2016; Burton-Johnson et al., 2017). To assign crustal heat production (A), we use the geological observations and crustal segmentation, as described in the previous section. We divide the crust into three classes according to stabilisation age: Archean-Paleoproterozoic, Meso-Neoproterozoic and Phanerozoic (Janse, 1984; Begg et al., 2009; Jaupart and Mareschal, 2013; Jaupart et al., 2016). For each class, an average heat production range is applied from Jaupart and Mareschal (2013). Crustal thickness is constrained from seismology (An et al., 2015a) and shown in Figure 4D. Details of the classification are given in Table 2.

We use the segmentation in Figures 3D, 4C to calculate new heat flow maps based on geophysical and geological input data using the methods described in the previous section. The resulting steady-state heat flow and associated uncertainties for the approach used, are shown in Figure 5. This provides an illustration of the further ability to compute output based on data of different types. Figure 5A shows our new mapped heat flow estimate, AqSS.ea, at continental scale. Figure 5C shows a regional equivalent for the Wilkes Subglacial Basin, AqSS.wsb, as an illustration of working at higher resolution. Calculated uncertainties are shown in Figure 5B, for East Antarctica, and Figure 5D for Wilkes Subglacial Basin.

3.2.6. Appraisal of the Steady-State Heat Flow Model, AqSS, and Previous Models

Our final set of functionality examples illustrate using the framework to appraise alternate models for a given parameter. Figure 6A compares AqSS, minimum and maximum values, with earlier published models and calculated heat flow from borehole measurements in western part of Australia (compiled by Hasterok, 2019). The Australian dataset includes transient and shallow processes, that are not captured in AqSS nor some of the other geophysically derived estimates.

FIGURE 6

FIGURE 6. Comparison of heat flow models enabled by the 3D model and framework. (A) Distribution of heat flow values. For East Antarctica, distributions from Fox Maule et al. (2005), An et al. (2015b), Martos et al. (2017), example heat flow values for AqSS (derived from values mapped in Figures 5A,B) as minimum estimate (blue line) and maximum estimate (red line); for Australia, distribution of actual measurements in southern and western Australia compiled by Hasterok (2019). (B) Heat flow model from seismic data (An et al., 2015b) minus AqSS. (C) Heat flow model from magnetic data (Martos et al., 2017) minus AqSS. Subtracting AqSS, which is a steady-state heat flow model, from published maps of total heat flow indicates non-steady-state contributions to total heat flow.

Figures 6B,C show examples of comparing two observation-derived datasets with a constructed reference model to inform the discussion of lithospheric properties. We show An et al. (2015b) and Martos et al. (2017) heat flow maps minus steady-state heat flow from AqSS. These two alternative results are effectively the additional heat flow likely generated from neotectonic and other non steady-state processes, such as recent rifting, volcanism and orogenesis.

3.2.7. Variation of Thermal Gradients With Depth

Figure 7 illustrates an example of extracting the variation of a property with depth. We show thermal gradients from locations in West and East Antarctica as a Gaussian kernel density estimate (KDE), including seismic-derived temperatures (An et al., 2015b) and magnetic-derived Curie temperature depth, including uncertainty bounds (Martos et al., 2017). The KDE is calculated over the depth dimension for East and West Antarctica separately. We also include uncertainties when defining the kernel size. In West Antarctica, the example is from Lake Whillans, the location of one of few direct measurements of heat flow in Antarctica (Fisher et al., 2015). In East Antarctica, the example is from Dome C. The location maps, showing West and East Antarctica, are obtained by importing a polygon vector to use as a factor (inset in Figures 7A,B).

FIGURE 7

FIGURE 7. Illustration of framework capability to extract depth profiles for model comparison. Thermal model of the lithosphere, populated with data from Antarctic heat flow models for West and East Antarctica reduced to kernel density estimations (KDE). Temperatures derived from seismic data, An et al. (2015b), in black contours showing highest concentration of thermal profiles. Depth to Curie temperature isotherm with uncertainty derived from magnetic data (Martos et al., 2017) in green contours. Surface and subglacial elevation from Fretwell et al. (2012) and subglacial temperature from van Liefferinge et al. (2018) in red at the surface. KDE Gaussian kernel for mantle temperatures set to 100°C/10 km, for Curie temperature isotherm 25°C/2 km and for surface 5°C/0.1 km. Plotted profiles in red show two examples locations of 1D temperature models using combined input. The subglacial heat flow is proportional to the gradient of temperature and the thermal conductivity in the upper crust. To facilitate KDE, only every fifth grid cell is computed. The figure is cropped at 250 km depth. Insets show sampled area.

The contours show the range of allowed values and how the two models, An et al. (2015b) and Martos et al. (2017), compare in depth section. The profile of temperature with depth varies over a large range for both example locations (Figure 7 red line), and when an average kernel is displayed (Figure 7 gray contours). This result, and the use of the 3D grid and framework in comparing models and sensitivity to different parameters, is further discussed below.

4. Discussion

We first outline the most significant limitations of the 3D model and framework, and then discuss aspects of our newly generated heat flow example, as an exemplar of how the research environment might be used.

4.1. Limitations

There is a trade-off between resolution and computational expense for any numerical model. Moreover, numerical stability is, in general, required for grid-based calculations. The continental scale model in 20 $\times$ 20 km grid, is presented as an example that is too coarse to contain and represent detailed observed geology and finer crustal geophysics. In terms of continental scale heat, the segmentation used to estimate the likely crustal heat production is not sufficient for ice sheet models that depend on heat transfer on a fine scale (van Liefferinge et al., 2018). The second provided example of the Wilkes Subglacial Basin in 2 $\times$ 2 km grid is more detailed in some areas, but includes interpolations from coarse data, and hence, the resolution appears finer than the data used. The open framework (Stål and Reading, 2020) facilitates a transparent workflow where the impact of, for example, model resolution can be tested.

The model functionality allows for the inclusion of uncertainty values matching each dataset. Therefore, the impact of the noted limitations can be mitigated. The model can be realized with a desired extent, resolution and data content to suit the needed outcome and stage of research. In this contribution, we include the uncertainties provided with the datasets. Those metrics may not cover the true uncertainty of the datasets, when resolution and artifacts from the methodology are considered. The strength of the framework is that the impact of such concerns can be understood as data coverage improves.

4.2. Insight From Examples

The heat flow estimate exemplifies how our multidimensional and multivariate grid may be used to combine input data of different types, and execute calculations across the grid. This provides, we hope, a constructive approach to reconcile the differences between published heat flow models for Antarctica (Figure 2C).

The comparison of the results from magnetic and seismic studies provides new insight into deep Earth properties since both approaches estimate temperature gradients, but using different methods. The differences in Curie temperature depths from seismic (An et al., 2015b) and magnetic (Martos et al., 2017) studies are larger in East Antarctica than in West Antarctica (Supplementary Material). These observations imply properties of the lithosphere such as fluid content and heterogeneous heat production that are not captured in the methods used. Compositional variations and presence of fluids impact the seismic wave speed and hence estimated temperatures (Hirth and Kohlstedt, 1996; Goes et al., 2000; Haeger et al., 2019). Magnetic models depend on a simplified crustal thermal and magnetic structure. As an example of a departure from the assumed case, shallow felsic intrusions can provide a large contribution to the surface heat flow, and this could be observed as a deeper Curie temperature isotherm because removal of radiogenic heat producing material facilitates cooling of the lower crust (Jaupart et al., 2016). Figure 7 highlights the large range and uncertainties involved in present heat flow estimates and also illustrates the much steeper thermal gradient in the crust compared to the upper mantle. We note that thermal conductivity generally decreases with temperature (Xu et al., 2004; McKenzie et al., 2005). However, geothermal heat is not lost rapidly through the crust, so crustal heat production must have a large influence on geothermal heat flow at the surface. New outputs such as Figure 7 show how the limitations in available evidence give rise to temperature changes with depth in the upper mantle that are, taken together, implausible. For example, a temperature decrease with depth is highly unlikely in stable lithosphere. The valuable studies that we have compared note their underlying assumptions and logical simplifications. Our new model and framework allows the implications of such simplifications to be better understood.

We have introduced a new conceptual heat flow model, AqSS, where we base the calculations on the energy balance of the lithosphere, rather than estimated temperature gradients. Our method represents a new approach in the Antarctic context and and uses a reduced number of assumptions. With negligible heat generated in the lithospheric mantle (An et al., 2015b; Jaupart et al., 2016; Martos et al., 2017), the Moho steady-state heat flux must be equal to the flux at the lithosphere-asthenosphere boundary. For old and stable crust, the mantle component of the heat can be reduced to a low and constant value in the range between 10–20 $m W m^{- 2}$ (Roy and Rao, 2003; Michaut et al., 2007; Jaupart et al., 2016), however, in more dynamic regions with thinner lithosphere, we need to include the non-steady-state contribution due to, e.g., tectonism (estimated from a geothermal gradient, but understanding the thermal properties in the crust as discussed above). AqSS provides us with an initial model that maps stable regions of the Antarctic interior. We then estimate the amount of transient (non steady-state) heat by subtracting the steady-state model from comprehensive models. This difference highlights dynamic regions in West Antarctica (Figures 6B,C). Including dynamic Earth processes ideally requires that not only crustal geology, but also hydrology, constraints from glaciology and the dynamic mantle are fully incorporated. Our framework, we hope, enables current and future progress toward that goal.

4.3. Use Cases for the 3D Model and Software Framework

The main use cases for Antarctic research, with an emphasis on interdisciplinary studies of the interaction of the solid Earth and cryosphere, are listed below:

(1) Computing results based on geophysical datasets. A broad range of datasets can be combined in the same frame and uncertainty bounds included, as illustrated in this contribution. The extensive toolboxes from, e.g., the Python ecosystem are available for modeling and analysis. Import, export and visualisation functions simplify the workflow. Supplementary Figure S4 shows the potential for experimentation in data visualisation.

(2) Combining geophysics and geological constraints, and making use of the merged result in ongoing calculations, as illustrated in this contribution. Constraints from glaciology could potentially be included in the same way, e.g., as a constraint on shallow processes to facilitate discussion of heat flow estimates for given regions.

(3) Appraisal of models. Comparisons between datasets, or calculated differences, can provide insights that are beyond the potential of the individual contributing studies, again, as we have illustrated in this contribution.

(4) Working with uncertainty and probabilistic methods. With the large uncertainties involved in Antarctic solid Earth research, probabilistic tools are essential to progress in the understanding of the Antarctic lithosphere. A productive way forward is to embrace the uncertainties and build probabilistic models (e.g., Stål et al., 2019). The computational framework that is presented here is well-suited to this task and provides an environment where data and associated uncertainties, probabilities and likelihoods can be processed.

(5) An enabling capability for the international research community. Building robust models of the Antarctic crust and upper mantle is a community effort, that will be refined incrementally with additional data. When a specific research product is desired, e.g., a reference heat flow map to include in ice sheet models, we can now draw constraints from multiple studies and/or easily test a range of alternative maps.

5. Conclusions

We present a new 3D grid model and framework: a computing environment tailored to interdisciplinary research. The software framework is easy to use, allows geophysical and geological data to be combined, and provides a virtual laboratory to develop and test, for example, solid Earth models. The model points directly to published data sources and the data contained can easily be updated.

This contribution aims to facilitate progress in Antarctic research concerning solid Earth-cryosphere interaction. Physical property maps and grids, of utility to studies of glacial isostatic adjustment, geothermal heat and the shaping of topography can be performed; bridging between the solid Earth and cryosphere research communities. The usage examples that we provide include a conceptually new steady-state heat flow map based on the energy balance of the lithosphere for comparison with maps based on modeled thermal gradient.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: Code and output products to be made available from Pangea: https://doi.pangaea.de/10.1594/PANGAEA.918549. Code includes updated URL to all used datasets. Code also published on Zenodo: https://doi.org/10.5281/zenodo.3775166.

Author Contributions

TS developed the software, built the 3D model, generated the examples and wrote the first draft text. AR guided the overarching research direction and advised on the geophysics. JH advised on the geology. SP advised on the interdisciplinary context. JW advised on the plate tectonics and basin geoscience. All authors contributed to revising the text.

Funding

This research was supported under Australian Research Council’s Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001) and the Center for Southern Hemisphere Oceans Research, a joint research center between QNLM and CSIRO.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This research is a contribution to the SCAR SERCE program.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.577502/full#supplementary-material

References

Adie, R. J., and Adie, R. J. (1977). Earth sciences: the geology of Antarctica: a review. Philos. Trans. R. Soc. Lond. B Biol. Sci. 279, 123–130. doi:10.1098/rstb.1977.0077