## Abstract

Global geophysical observations constrain all theories of terrestrial dynamics. We jointly interpret EGM2008 gravity, RET2014 topography and the Global Centroid Moment Tensor database from a structural point of view. We hypothesize that lateral variations of gravity and topography reflect the scale-dependent competence of rocks. We compare the spectral and spatial characteristics of the observed fields with structural predictions from the mechanics of differential grade-2 (DG-2) materials. The results indicate that these viscoelastic materials are a powerful tool for exploring dynamic processes in the Earth. We demonstrate that the known spectral range of Earth's gravity and topography can be explained by the folding, shear banding, faulting and differentiation of the crust, lithosphere and mantle. We show that the low-amplitude long-wavelength bias apparent in the disturbance field can be explained by perturbations to Earth's overall ellipsoidal shape, induced by internal slab loading of the mantle. We find by examining the directional isotropy of the data that the zonal energy in Earth's gravity disturbance is maximized about an axis coincident with the shape-perturbation minimum. The symmetry of tectonic features about this axis, extending from eastern Borneo to Brazil, and its coincidence with the equator suggest the coupling of current plate motions to true polar wander.

Gravity, topography and seismicity reflect the historically integrated response of our rocky planet to surface tractions and internal stresses. The purpose of this contribution is to interpret these global datasets from a structural point of view. Four questions come to mind: Is there evidence in global high-resolution gravity–topography spectra for competent layering in Earth's radial structure? How is the azimuthal anisotropy of Earth's gravity and topography (first identified by Bills & Lemoine 1995) related to the dynamics of global plate motions? To what extent do Earth's topography and near-surface density variations mask the dynamics of the deep interior? Are extreme variations in effective viscosity necessary to drive global plate motions?

We hypothesize that the observed variations of gravity and topography reflect the scale-dependent competence of rocks. This hypothesis was considered earlier by Gilbert (1889), Barrell (1914) and Watts (2001). Its evaluation was based on gravity measurements adjusted for, or predicted by, simple isostatic models like Airy or Pratt compensation or elastic plate flexure. While it is clear that topographical loads are isostatically compensated to some degree, the parameters needed for these models (e.g. density, layer thickness, age at loading, mechanical properties) are either poorly known or inconsistent with the observed fabrics of deformed rocks.

For example, the elastic plate model has been used extensively to interpret gravity and topography data (Watts 2001). This work demonstrates conclusively that the oceanic lithosphere possesses finite rigidity. However, significant limitations of this model are that it offers no explanation for localized deformation at any scale, and that it predicts smooth elastic distortions with wavelength *λ* to thickness *z* ratios, *ξ* = *λ*/*z*, of 25–30 or more.

These results have led some of the earliest proponents of plate tectonics to argue that the lithosphere is effectively only about 25 km thick. This is at odds with estimates of plate thickness based on thermal subsidence of the seafloor with age (Parsons & Sclater 1977; Stein & Stein 1992), as well as with recent broadband seismic observations of a shear-zone-like structure in the continental lithosphere at depths ranging from about 80 to 130 km (Cooper *et al.* 2017).

Subsequently, by combining the unprecedented high-resolution and global consistency of modern geophysical data with the mechanics of differential grade-2 (DG-2) materials (Rivlin & Ericksen 1955) (see the subsection ‘Diharmonic structural theory’ in Appendix A), we re-evaluate the scale-dependent competence hypothesis. Specifically, we exploit the fact that these non-linear viscoelastic materials exhibit a range of responses to loading, from homogeneous deformation to folding to shear banding to effective rigidity, with transitions in the range *ξ* and thermomechanical competence, the ratio of thermal *κ* to mechanical *χ* diffusivities (Patton & Watkinson 2010), to probe the depth extent of competent layering and shear localization in the Earth.

We demonstrate that the known spectral range of Earth's gravity and topography can be explained by the folding, shear banding, faulting and differentiation of the crust, lithosphere and mantle. We show that the low-amplitude long-wavelength bias apparent in the disturbance field can be explained by perturbations to Earth's overall ellipsoidal shape, induced by internal loading of the mantle by subducted slabs. We find by examining the directional isotropy of the data (Bills & Lemoine 1995) that the zonal energy in Earth's gravity disturbance is maximized about an axis coincident with the shape-perturbation minimum. The apparent symmetry of active tectonic features about this axis, extending from eastern Borneo to Brazil, and its coincidence with the equator suggest the coupling of global plate motions to true polar wander. We interpret crustal detachments, the subduction interface and mantle seismic discontinuities in light of our findings. Finally, we compare and contrast our recent work on the poloidal–toroidal problem for DG-2 materials with earlier work on Newtonian and power-law viscous fluids. Together, these results indicate that DG-2 materials are a powerful tool for exploring dynamic processes in the Earth.

## Gravity, topography and seismicity

Gravity pervades the physical environment. Revealed by the plumb bob and the carpenter's bubble, it defines the vertical and the horizontal. Consistent with Newton's mechanics (*c.* 1686), its spatial variations reflect both the shape of and distribution of mass within the Earth.

Experience shows that fluids at rest in a gravity field naturally seek an equipotential configuration, free of shear stresses. As such, mean sea level – the geoid – serves as an intuitive vertical datum for measurements of height and depth. The geoid is approximated to within 10^{−}^{5} of Earth's mean radius by an ellipsoid of revolution (Kaula 2000).

The equatorial radius of this reference ellipsoid is about 21 km longer than its polar radius. This flattened shape dominates the gravity field, such that gravity at the equator is about 5.6 gal less than at the poles (Fig. 1a). Gravity variations residual to Somigliana's normal field (*c.* 1929), which accounts for the ellipsoidal shape effect, are no more than 10^{−}^{3} of the mean field. Although small, these gravity disturbances (Fig. 1b) have amplitudes well above the noise of modern gravimetric techniques.

Prior to the age of artificial satellites, a geocentric reference system could not be realized (Hoffmann-Wellenhof & Moritz 2005). Although accurate regional reference systems existed, they all differed by unknown offsets of their coordinate origins from Earth's centre of mass. Consequently, a global comparison of heights was not strictly possible. Because satellites orbit the centre of mass and can be tracked accurately from the ground, they provide the observational basis for truly global reference systems. Today, for example, the global positioning system, fully operational since 1995, gives accurate positions in an Earth-centred Earth-fixed reference frame, without first being referenced to the geoid. Earth's geometry is measureable directly from orbit.

Earth's solid physical surface is not an equipotential of gravity. The height of the land and depth of the ocean comprise what we call topography. Its variations with respect to the geoid are no larger than 10^{−}^{3} of the mean planetary radius. Historically, for reasons of accessibility, the height of the land in near-coastal regions was known with reasonable accuracy. For the same reason, knowledge of the height of mountainous regions and depth of the ocean was relatively poor.

It was not until the mid-twentieth century that echo soundings from ships at sea revealed the presence of the globe-girdling system of mountain ranges beneath – the oceanic ridges. Combined with magnetometer surveys of the seafloor and the reversal chronology of the geomagnetic field, their morphology could be understood as the creation of new crust by magmatism in an axial rift and its subsequent lateral motion away from the ridge axis. This is the ‘seafloor spreading hypothesis’ (Dietz 1961; Hess 1962).

Seismic activity has been known since antiquity (Udías 1999). However, a precise intrumental record of related motions on a global scale did not exist prior to the establishment of the worldwide standardized seismograph network in 1962. Although initially intended to monitor nuclear weapons testing, it clearly had broader scientific value; both as a means of measuring the elastic structure of the planet (Dziewonski & Anderson 1981), and revealing the geographical distribution and source characteristics of earthquakes.

Deep dipping zones of seismicity had already been found beneath magmatic arcs, suggesting that these regions were sites of convergent motion and recycling of mass into the mantle (Benioff 1954). The coincidence of shallow earthquakes with oceanic ridges and fracture zones, and their first-motion solutions, when combined with the ‘transform fault hypothesis’ (Wilson 1965), further supported the idea that the seafloor was diverging from ridge axes, parallel to fracture zones.

Recognition that these motions could be described by rotations of rigid shells on a sphere was widespread (McKenzie & Parker 1967; Le Pichon 1968; Morgan 1968). Broad aseismic expanses were ‘plates’ bounded by narrow seismically active ‘margins.’ This established the ‘mobile lithosphere hypothesis’ (Isacks *et al.* 1968), spurring wide interest in seafloor studies and fuelling speculation about the dynamics of global plate motions.

Application of these ideas to deformation in mountain belts soon followed (Dewey & Bird 1970). The early classification of convergent margins into ocean–ocean (OO), ocean–continent (OC) and continent–continent (CC) types continues to be useful. However, the data supporting these distinctions are more diverse than at oceanic ridges and are decidedly ‘geological’: that is, not readily reduced to simple physics-based concepts. A legacy of the last 60 years is the tremendous volume of detailed work by specialists characterizing the pressure–temperature–time and structural histories of these prominent seismically active regions.

Given the fact that Earth is rather round, spherical harmonics provide a convenient set of orthogonal basis functions for representing global gravity and topography observations. Spherical harmonic models comprise tables of coefficients, organized by degree *l* and order *m*, obtained through Fourier analysis of gridded observations. The dimensionless wavenumbers *l* and *m* are associated with latitude and longitude in a spherical-polar coordinate system. Spherical harmonics are therefore convenient for studying the spectral characteristics of the fields (Wieczorek 2015). Through the application of accurate recurrence relationships – as implemented in SHTOOLS, for example – these models are readily converted to Cartesian grids for visualization as maps.

A rule of thumb for interpreting these models is Jeans’ relation: *λ* is dimensional wavelength and *R* is the planetary radius. Long wavelengths are represented by low-degree harmonics, while short wavelengths are given by high-degree ones.

Gravity observations had only a peripheral impact on plate tectonics, not for the lack of scientific relevance but rather because of spatial resolution (McKenzie 1967; Kaula 1969). The formal spatial resolution of a spherical harmonic model is indicated by the largest degree, *l*_{max}, at which model coefficients of all degrees and orders are present and can be estimated using Jeans’ relation. However, because the construction of such models involves least-squares estimation, the actual resolution varies with geographical location as dictated by the point density of *in situ* observations and/or the orbital characteristics of satellite-based sensing systems (Floberghagen 2002).

A decade after the October 1957 launch of *Sputnik-1*, the global gravity field was resolved to *l*_{max} = 16 (Gaposchkin & Lambeck 1971): that is, wavelengths greater than about 2500 km. Regional gravity fields were known to a higher resolution, but their interpretation suffered from many of the same issues involved in topography interpretation. By 1996, global gravity models of *l*_{max} = 360 were possible, combining long-wavelength satellite-tracking observations with gridded block means from terrestrial gravimetry (Nerem *et al.* 1995). The intensive numerical work needed to realize these models was possible because of advances in digital technology. Today, combined gravity models of *l*_{max} = 2159 and higher (EGM2008: Pavlis *et al.* 2012, 2013) and complimentary topographical models (RET2014: Hirt & Rexer 2015) are available for download from the Internet.

Seismicity catalogues are available from numerous agencies, also via anonymous download. Here, we use the global centroid moment tensor (CMT) database (Ekström *et al.* 2012). CMT inversion incorporates a wider range of seismic frequencies than hypocentre location. As a result, fewer CMT solutions are available with generally larger moment magnitudes. Still, the variation of moment release v. depth from CMTs appears more ‘natural’ than from hypocentres. Differences between these two methods can be attributed to their relative sensitivity to lateral heterogeneity in seismic source regions, systematic data-selection biases and their operational goals.

Modern high-resolution geophysical datasets provide crucial constraints for models of planetary structure and dynamics. Combined with predictions of rock-like modes of deformation from the mechanics of DG-2 materials (see the subsection ‘Diharmonic structural theory’ in Appendix A), we examine the global spatial and spectral distribution of rock folding and faulting.

## Characteristic rock structures

A structural renaissance dawned with the publication of *Folding and Fracturing of Rocks* (Ramsay 1967), which focused attention on the importance of detailed mapping for understanding Earth's deformation fabrics. Incorporating insights from rock mechanics, *An Outline of Structural Geology* (Hobbs *et al.* 1976) expanded this new approach to the collection, modelling and interpretation of structural data.

Over the years, we have found the concept of characteristic structures to be useful for interpreting deformation. Examples include linear-planar tectonite shape fabrics, single and multilayer folds, contact strain zones, and shear zones. Each of these has an analogue at the scale of plate tectonics. This is indeed a case where the small scale reflects the large.

Linear-planar tectonite shape fabrics principally reflect the amount of strain of the rock. This is often expressed as a fold axial-plane cleavage or schistosity, with or without a linear fabric. In some rocks, the effects of fracturing, recrystalization and/or partial melting are recognized. However, in general, it is a strain fabric, sometimes modified by subsequent thermal events, that is observed. This is true across different ages, structural settings, lithologies and metamorphic grades. Often, a series of superposed fabrics reflecting changes in kinematics and metamorphic conditions are preserved.

Folds, both single layer and multilayer, are characteristic ductile structures with wavelength-to-thickness ratios typically in the range 2 < *ξ* < 7. This applies to all sizes of folds, and hints that a scaling relationship might be in play. Multilayer folds may be harmonic, disharmonic or polyharmonic as a function of layer thickness and competence contrast between layers. Folds exhibit contact strain zones, typically as wide as an arc length from the most competent layer, or layer package.

Interpretation of gravity data from orogenic belts has revealed the presence of long-wavelength multilayer folds, often with polyharmonic folding (Burov *et al.* 1993; Jin *et al.* 1994; Caporali 2000). This style of folding was recognized and described at outcrop scale by Ramsay and his graduate students (e.g. Hudleston 1973; Ramsay & Huber 1987; Watkinson & Thiessen 1988). Using the characteristic attributes of these folding systems, from small to large scale, we find wavelength-to-thickness ratios of about 4.

Shear zones and shear bands form in different rock types, mountain belts and metamorphic grades. They show similar geometries, strain gradients, grain-size diminution and high strain localization. These features have been observed in exposures of upper-mantle to low-grade metamorphic rocks and also in exhumed xenoliths from deeper mantle sources (Boyd 1973).

Field structural analysis of rocks exposed in the deeper levels of mountain belts has demonstrated a surprising degree of structural order. There are high degrees of homogeneous domain symmetry over significant areas. Early formed structures are frequently observed and preserved, and they control the attitudes and styles of superimposed later episodes, phases and progressively formed structures. This is solid-state ductile deformation.

Shear localization fabrics from electron microscope images, thin sections and regional shear zones of tens-of-kilometres width all show grain-size reduction and more intense strain within the shear zone, compared to without. Small parasitic folds have similar geometries and orientation to the very large-scale nappe structures in polyharmonic fold systems (e.g. Helvetic nappe piles of the Swiss Alps 10–15 km thick).

Arguably, the most significant patterns in deformed rocks are localized zones of high strain and/or major displacement: that is, ductile or ductile–brittle to brittle shear zones. Ramsay & Graham (1970) spurred a global recognition of the importance of these features. Their strain patterns range from transtension to simple shear to transpression, with volume changes, and their study led to an understanding of strain compatibility and grain-size reduction. These are now recognized in all kinds of rocks, preserved in very old orogenic belts, and occur in extensional, strike-slip and contractional kinematic regimes.

Broadband seismic instrumentation has now been deployed in many continental regions. Teleseismic shear-wave splitting and receiver function analysis of these data resolves the presence of numerous shear-zone-like reflectors in the deep continental lithosphere and upper mantle. Cooper *et al.* (2017) observed that at depths of between about 80 and 130 km lithospheric rocks are ‘riddled’ with seismic reflectors.

In summary, characteristic structures show that rocks deform irreversibly in the solid state. While deformation changes the shape of individual mineral grains and domains, it does not typically obliterate pre-existing fabrics. Structure persists at all scales because of rock competence.

## DG-2 materials

A DG-2 material is an ideal ‘continuum mechanical’ model of non-linear viscoelastic behaviour exhibiting certain geologically interesting responses to loading. Specifically, it manifests both folding and shear-banding modes of deformation in pure shear, parameterized by the ratio of thermal *κ* to mechanical *χ* diffusivities (Patton & Watkinson 2010). We refer to this parameter as ‘thermomechanical competence’.

Thermodynamics shows that a DG-2 material will deform adiabatically and homogeneously for low values of thermomechanical competence. Furthermore, it will behave rigidly for *κ*/*χ* ≥ 1 (Patton & Watkinson 2013). Because thermomechanical competence is inversely proportional to normalized wavelength *ξ* = *λ*/*z*, these transitions from homogeneous deformation to folding to shear banding to rigidity occur at

The transition between folding and shear banding, at *κ*/*χ* = 1/2, is scale-invariant. This powerful scaling property arises from the fact that the two-dimensional stress balance of DG-2 materials is governed, at leading order in a small parameter *ε*, by a diharmonic equation that admits these instabilities as general solutions. Mathematically, a self-consistent rescaling transformation exists, which resolves what would otherwise be a bald singularity (Patton & Watkinson 2010). This feature, by itself, recommends DG-2 materials over many other rheological models we might consider, including the popular Maxwell viscoelastic model.

As with all linearized solutions, it is important to bear in mind that they are strictly correct only if the small parameter satisfies the relationship 0 < *ε* ≪ 1. For folds, *ε* is given by the ratio of fold amplitude to fold wavelength; while for shear bands, it is given by the ratio of band width to band length (Patton & Watkinson 2010).

Recently, we extended this work to four-dimensional spacetime by examining poloidal–toroidal coupling in a self-gravitating DG-2 material body (see the subsection ‘Poloidal–toroidal coupling’ in Appendix A). The toroidal field, at leading order, is governed by an inhomogeneous diharmonic equation. This result is independent from the assumption of incompressibility made in our earlier work (Patton & Watkinson 2005; Patton & Watkinson 2010). In the present context, a suitable small parameter is given by the deviation of solid topography from the geoid.

Toroidal motions are rigid-body rotations about Euler vectors on a sphere (local strike-slip). They are non-linearly coupled to kinematic expansion, and they occur by simple shearing along an equipotential surface. The ‘weight’ of the material above the shear zone provides the normal stress needed to suppress the otherwise dominant dynamic rescaling mechanism (Patton & Watkinson 2010).

On the other hand, poloidal motions are upwards and downwards mass movements (local dip-slip) explicitly coupled to kinematic expansion. Hence, because we need not assume incompressibility, we find that poloidal and toroidal motions can coexist in a self-gravitating DG-2 material body. Moreover, these motions are coupled via terms proportional to the first normal stress coefficient (Rivlin & Ericksen 1955), which is subsumed by the mechanical diffusivity *χ*. Material gradients, like the strong lateral variations of effective viscosity assumed in conventional geodynamics, are not required.

Because Earth's topography, *h*, is only 10^{−}^{3} of the mean planetary radius, *R*, it conforms to the linearization conditions imposed above: that is, 0 < *ε* = *h*/*R* ≪ 1. Hence, it is reasonable to expect that the scaling relationships identified in our earlier study of fold-like and fault-like structures in DG-2 materials should be expressed in terrestrial gravity and topography fields as a function of wavelength.

The diharmonic equation can be written as a wave equation (see the subsection ‘Seismogenesis’ in Appendix A). Similar equations model the kinematics of fault rupture (Aki & Richards 2002). In-plane motions generate P-SV waves, while anti-plane motions generate SH waves. Equating dimensionless coefficients between these models, we find that crack-rupture speed, as a function of thermomechanical competence, is greater than the shear-wave speed for *κ*/*χ* > 1/2 (Fig. 3). This is consistent with spontaneous instabilities occuring in this range of thermomechanical competence.

In elastodynamic rupture models, stress and velocity singularities at the crack tip are controlled by a hypothetical anelastic cohesive zone of small but finite size (Aki & Richards 2002). The size of this zone plays a role similar to the rescaling modulus *|α|* in DG-2 materials theory (Patton & Watkinson 2010). However, because ruptures predicted by the stress balance of DG-2 materials are localized within a larger thermomechanical continuum, it is reasonable to expect that the same scaling relationships as applied in the analysis of gravity and topography would be expressed in the three-dimensional distribution of earthquakes. Despite their similar rupture dynamics, the large-scale behaviour of a DG-2 material is quite different from a Hooke solid. In the large, gravity prevails.

## Data characterization

### Statistical measures

The fact that gravity and topography are measured only at discrete points on the Earth's surface dictates that they be treated as random variables. Consequently, statistical methods are required to interpolate between measurements, extrapolate into poorly observed areas, and account for uncertainties and errors. Compounded with diverse gravimetric and geodetic methods, we begin to appreciate the difficulty of the inversion problems that lead to useful high-resolution model products like EGM2008 (Pavlis *et al.* 2012) and RET2014 (Hirt & Rexer 2015).

Keeping these issues in mind, we shall treat gravity and topography as known fields. Of these two stochastic processes, topography is probably better known and therefore serves as our reference signal. We then interpret the empirical relationship between gravity and topography based on geologically reasonable assumptions.

Such has been the modus operandi since the first efforts to determine Earth's shape. Returning from Condamine's expedition, Bouguer (*c.* 1749) expressed puzzlement at the fact that the attraction of the mountains near Quito was measurably less than expected. This phenomenon, called isostatic compensation, plays a significant role in Earth's gravity field (Watts 2001).

Two independent measures of the relationship between gravity and topography are the correlation *γ*(*l*) and the admittance *Z*(*l*) (Wieczorek 2015). The correlation is usually treated as a normalized cross-variance that measures both the strength and direction of this relationship. It is biased downwards by uncorrelated noise in the gravity signal. The correlation is often combined with confidence limits to indicate the significance of the relationship (Bills & Lemoine 1995).

The admittance is a linear transfer function that measures how well the topography predicts gravity (Dorman & Lewis 1970). It is usually reported in units of mgal/m. Assuming transverse isotropy, this function is computed by averaging over spherical harmonic order at each degree and thus remains sensitive to anisotropy with depth.

Consistent with our interest in rock competence, we normalize the empirical admittance *Z*_{u}(*l*) (equation 5.8 of Watts 2001). Model parameters include a mean surface density of 2450 kg m^{−3} and a mean ocean depth of 4800 *±* 200 m. The result, a response function

Directional isotropy of the gravity and topography fields can also be measured statistically. Bills & Lemoine (1995) showed how directional derivatives in latitude and longitude can be used to formulate an isotropy measure *R*(*l*), defined by a weighted ratio of north to east mean slopes. Because lines of latitude define spherical harmonic zones, while lines of longitude define spherical harmonic sectors, this isotropy measure is actually a weighted ratio of zone to sector mean slopes. No averaging over spherical harmonic orders is done. If the isotropy hypothesis was universally true, then the isotropy measure would equal unity at all wavelengths.

In the absence of sampling artefacts, a scalar function on a sphere is independent of the coordinate system. Hence, we are free to rotate the coordinate system to any arbitrary orientation and measure the zone to sector isotropy about that axis. By sampling the data over a range of orientations, we find an overall spread of anisotropy, which can be interpreted tectonically.

Spherical harmonic rotations, which scale as *l*^{3}, are relatively memory-intensive compared to map processing, which scales as *l*^{2}. Thus, while the memory of off-the-shelf computers is adequate for manipulating two-dimensional arrays of double-precision floats to the full resolution of the data, it limits rotations to about *l*_{max} = 800. Consequently, our analysis is limited to that range.

Earthquakes are highly localized spontaneous shearing instabilities within the Earth that generate elastic waves. The radiated waves are detected and recorded at seismometers, which are typically part of global and regional networks. The inferred source characteristics of earthquakes (e.g. location, origin time, first motions, rise time, pressure/null/tension axes, moment magnitude) reflect the *a priori* unknown properties of the source region, as well as the medium through which the waves propagate. All of this information must be treated as random variables.

We use normalized depth-moment release (NDMR) curves to probe Earth's layer competence. They are an elaboration upon similar curves presented by Frohlich (1989). Specifically, we compute the sum of log seismic moment, log *M*_{0}, normalized by the number of events, *N*, in the geographical sample as a function of depth. This statistic measures the mean magnitude of events in the sample region. By comparing the NDMR of various tectonic subsets with the global average, we can study variations of layer competence in three dimensions.

NDMR curves for the NEIC Preliminary Determination of Epicentres (https://earthquake.usgs.gov) exhibit large spikes at 5, 10, 15 and 35 km depths. These are default depths for events with poor depth control. This catalogue is not ideal for probing competent layering at typical crustal depths. In contrast, NDMR curves for the CMT catalogue (http://www.globalcmt.org/) exhibit a more ‘natural’ appearance. We used the CMT catalogue in our work.

### Competent layering

We know from seismology that the Earth is radially stratified, with a largely solid-state crust and mantle overlying a liquid outer core (Dziewonski & Anderson 1981). We know this because shear waves propagate through solids but not through liquids, consistent with global patterns in converted seismic waveforms (Udías 1999). We also know from isostasy that the stronger lithosphere (Watts 2001) overlies the weaker asthenosphere (Barrell 1914). Is there evidence of other layers in Earth's gravity–topography spectrum?

Subsequently, we make the following two assumptions in our global study. First, we assume that the rocks comprising the crust and mantle can be modelled as DG-2 materials (Rivlin & Ericksen 1955). Therefore, to the extent that Earth's observed gravity disturbance (Fig. 1b) is related to the folding and faulting modes of deformation in competent layers, the structural characteristics of this ideal material provide a means to uniquely interpret the observed field.

Second, we adopt the radial structure and labelling convention established for the reference model ThERM (Thermomechanical Earth Reference Model: Patton & Watkinson 2010) to refer to layers of thickness *z* in the crust and mantle (Table 1). Specifically, we focus on the brittle crust (H1−H4), the lithosphere (L1−L4) and the upper mantle (M1−M4). We augment this list with the mantle itself, layer CMB.

Combining these assumptions, we can predict the degree at which a change in deformation mode should occur for a layer of given thickness. We focus on the layered radial structure of ThERM because its least-squares misfit to PREM (Preliminary Reference Earth Model: Dziewonski & Anderson 1981) is minimized for a 100 km thick ‘lithosphere’ (L4). Specifically, by eliminating dimensional wavelength between the normalized wavelength and Jeans’ relation we find:

Rounded to the nearest degree, these ThERM-scaled transitional harmonics are reported in Table 1.

### Origin of long-wavelength gravity disturbances

It is well known that gravity and topography are poorly correlated at long wavelengths (e.g. Kaula 1969). In contrast, the correlation is strongly positive and consistently exceeds the 99% confidence limit for wavelengths less than about 3000 km (*l* ≥ 13) (Fig. 4) (Bills & Lemoine 1995).

This dichotomy has led to the conclusion that the broad variations in gravity (Fig. 1c) have little connection with other surface phenomena (e.g. Kaula 1969; Panasyuk & Hager 2000; Flament *et al.* 2013). Kaula (1969, p 4807) lamented this fact, wishing that it were based on ‘considerably more geologic evidence’. Consequently, gravity studies have tended to emphasize either global ‘dynamic’ or regional ‘isostatic’ problems based on wavelength (Nerem *et al.* 1995).

The long-wavelength low-amplitude gravity disturbance (Fig. 1c) strongly influences the visual impression of Earth's gravity residual to the normal field (Fig. 1b). While this reflects the distinctly non-Gaussian distribution of the measured values (Kaula 1993), it does not explain the origin of the bias.

Given the phenomenon of plate tectonics, it is natural to seek an explanation that involves the upwards and downwards (poloidal) movements of mass in the mantle (Hess 1962). Subducted slabs, inferred from the distribution of earthquakes, penetrate to depths of at least 700 km and are well correlated with this disturbance (Hager 1984). Hence, it is puzzling that dynamic models that account for slab loading predict some 10^{3} m of topography (Panasyuk & Hager 2000; Flament *et al.* 2013). Scaled using the mean vertical gravity gradient of 0.3082 mgal/m, this implies gravity disturbances more than 4–6 times larger than those observed.

This in turn has led to the conclusion that the dynamic mantle signal is masked by the much greater noise of the near-surface layers (Panasyuk & Hager 2000). Presumably, this is due to uncertainties in models of isostatic compensation used to reduce surface observations, or to uncorrelated gravity signals that are discarded by the isotropic averaging procedure.

On the other hand, the fact that the long-wavelength bias can be removed simply by filtering degrees *l*<13 from the disturbance field suggests another explanation. Observe that Earth's correlation spectrum exhibits distinct lows at *l* = 2, 5 and 12 (Fig. 4). These lows occur on the long-wavelength side of the transitional harmonics predicted for layer CMB (Table 1). The negative correlation at *l* = 2 is explained by Earth's rotationally induced ellipsoidal shape (Fig. 1a). Similarly, the other lows suggest that the poor correlation between gravity and topography at long wavelengths is the result of folding and shear-banding deformation of the entire mantle, which is out of phase with surface topography. Subducted slabs in this case act as internal loads, emplaced via poloidal motion, which perturb the ellipsoidal reference shape of the planet. This is not a direct mass heterogeneity effect, but one consistent with the regional flexure of a competent layer called the mantle.

Scaling the long-wavelength gravity disturbance (Fig. 1c) by the mean vertical gravity gradient and referring it to the geoid, we obtain the shape perturbation shown in Figure 5. Added to the reference ellipsoid, this map gives a higher-order approximation to Earth's shape.

The dynamic range of the shape perturbation is *−*120 to 170 m. This is less than 10^{−}^{2} of the rotation-induced polar flattening. If indeed *Ut tensio, sic vis* (Hooke *c.* 1678), this suggests that the forces driving global plate motions are less than 1/100th of the centrifugal ones.

### Regional gravity disturbances

Joint analysis of EGM2008 (Pavlis *et al.* 2012) and RET2014 (Hirt & Rexer 2015) using a normalized admittance approach (Watts 2001) reveals a strong wavelength dependence in Earth's response to loading. The response function is generally consistent with isostatic compensation of topography at longer wavelengths, the absence of compensation at medium wavelengths and dynamic enhancement at shorter wavelengths (Fig. 6). Given the inverse distance dependence of the gravity potential, this pattern suggests a variation in material competence with depth. DG-2 materials theory provides us with the means to make a unique global interpretation of these patterns.

The lower part of Figure 6 shows a graphical representation of the ThERM-scaled transitions (Table 1). Coloured bars represent degree ranges where the predicted response of a given layer is dominated by homogeneous deformation, folding, shear banding or rigidity. Vertical placement of these bars mirrors layer thickness, where the thinnest layer is on top with progressively thicker layers beneath.

As shown earlier, the mantle responds to surface loading only at the longest wavelengths. Additionally, we see that ThERM predicts Earth to be less competent on the outside, particularly at long wavelengths. ThERM also predicts distinct ‘folding-only’ spectral ranges where shear banding and earthquakes should not occur. Similarly, it predicts distinct ‘mixed-mode’ spectral ranges where folding, shear banding and seismicity are to be expected. Given the characteristic importance of shear localization to rock deformation at all scales, we tabulate the degree ranges and implied depths of shearing (Table 2). Each of these depth ranges has interpretive value for geology.

ThERM does not explain every aspect of Earth's response to loading. Since admittance analysis combines information from two stochastic processes, it is unreasonable to expect that it would. Still, it is remarkable that all of the predicted transitions in this spectral range appear to be expressed in the response function. We submit that this establishes ThERM as a credible interpretive tool for understanding global plate dynamics and related phenomena. In the following subsections, we catalogue the geological insights and dynamic implications of these findings.

#### Compensated disturbances

Gravity is less than predicted by the uncompensated topography model for wavelengths greater than about 96 km (*l* < 416) (Fig. 6). ThERM suggests that compensation results from folding in the lower crust (L1−L2), and from folding and shear banding of the lithosphere (L3−L4) and upper mantle (M1−M4). Shallower layers of this model are relatively incompetent, transparent to topographical loads; while the lower mantle itself is rigid, not responding except at wavelengths greater than about 3000 km.

A map of the spectral band 13 < *l* < 417 reveals numerous tectonic features (Fig. 1d). The largest gravity disturbances are associated with convergent margins. The Andes and Himalaya have central highs flanked by broad lows. The Hawaiian–Emperor seamount chains also exhibit symmetrical patterns. In contrast, island arcs have a high on the overriding side flanked by a single low on the underthrust side. The smallest disturbances tend to occur in the ocean basins. Although there is some hint of oceanic ridges and fracture zones, their signatures are remarkably muted. Overall, the visual impression of the ocean basins differs markedly from the continents, suggesting different modes of compensation.

The muted gravity signature of the ocean basins can be explained by thermal isostasy. There are two types of models for this phenomenon. The first invokes the progressive cooling and thickening of a boundary layer (Turcotte & Oxburgh 1967). The second invokes the cooling of a plate of given thickness possessing an isothermal lower boundary (McKenzie 1967; Parsons & Sclater 1977; Stein & Stein 1992). For relatively young seafloor, these models predict the same thing: viz. the mean depth of the ocean should increase as the square root of seafloor age. However, for older seafloor, their predictions diverge, the former calling for continued thermal subsidence, while the latter implies flattening of the seafloor.

Observed ocean depth generally follows a square root of age trend, but for ages greater than about 60–80 myr it tends to flatten (Parsons & Sclater 1977). Consequently, the plate model can be used to estimate the thermal thickness of the oceanic lithosphere: Parsons & Sclater (1977) found *z* = 125 km, while Stein & Stein (1992) found *z* = 95 km. The latter value is practically equivalent to the L3 layer, while the former is one of the thermal modes from ThERM (Patton & Watkinson 2010). Using equation (1), a 125 km-thick layer should manifest diharmonic transitions at *l*∈{46, 160, 320}.

Note that the peak of the empirical admittance occurs at about *l* = 320. Furthermore, observe that the response function trends downwards over the range 320 < *l* < 417. This is spectral evidence in favour of the hypothetical isothermal boundary condition applied at the base of the plate model by Parsons & Sclater (1977). Thus, it appears that mature oceanic lithosphere influences Earth's gravity–topography spectrum in this range. Moreover, folding of oceanic lithosphere should affect the spectrum in the range 57 < *l* < 208 (Table 1). Measurements of trench–outer-rise flexural wavelengths are consistent with these values (Fig. 6, purple arrows).

Large-amplitude symmetrical disturbances in active mountain belts suggest regional folding of the lithosphere, while their coincidence with convergent plate margins suggests that localized poloidal shearing in the deeper mantle drives shortening (Benioff 1954). In the mixed-mode degree range 30 < *l* < 60, we see that the continental lithosphere (M1−M2) folds; while at deeper depths, we see the upper-mantle shear bands. Spectrally, this coincides with a smoother segment of the response function. The Andes and Himalaya appear to be short-scale orogenic loads nestled in flexures of the continental lithosphere above downgoing slabs.

Gravity disturbances in the continents appear to reflect tectonic history. Rifts and sutures are evident as elongated highs. Prominent examples include the mid-continental rift of North America, the Baikal Rift, the East Africa Rift System, the Urals, the Caledonide suture and the ancestral Rockies.

#### Uncompensated disturbances

Gravity is close to the uncompensated topography model for wavelengths in the range 29–96 km (416 < *l* < 1390) (Fig. 6). A map of this spectral band reveals major mountain belts, rift zones and seamount chains (Fig. 1e). Reasonably, this implies finite strength for many of Earth's major topographical features (Gilbert 1889; Barrell 1914; Watts 2001).

#### Enhanced disturbances

Gravity is enhanced with respect to the uncompensated topography model for wavelengths of less than about 29 km (*l* > 1389) (Fig. 6). ThERM suggests that enhancement is due to shear banding in the middle (L1) and upper (H3−H4) crust, and due to folding in the uppermost layers (H1−H2). Discernible features on a map of this spectral band are Earth's major orogenic belts (Fig. 1f): the Himalaya, the Andes, the Cordillera, the New Guinea Highlands and the Southern Alps.

A map of the Himalaya–Tibet region shows that most of the enhanced gravity signal occurs in the mountains (Fig. 7). It is also coincident with earthquakes in the 14–37 km depth range. Deeper seismicity is associated with structural syntaxes of the Burma arc and Hindu Kush. The distribution of these deeper events is consistent with ThERM. This short-wavelength-enhanced disturbance is nestled in a trough of the long-wavelength shape perturbation, and lies above a seismically inferred remnant of subducted oceanic lithosphere.

### Directional isotropy

Having found that Earth's overall ellipsoidal shape is dynamically perturbed at long wavelengths (Fig. 5), we wish to know if this is also expressed in competent structures at shorter wavelengths. The fact that global plate motions diverge from oceanic ridges and converge on subduction zones suggests that this is a good hypothesis, but how can we test it?

Bills & Lemoine (1995) demonstrated that Earth's isotropy spectrum, in the standard geographical projection, is close to unity for *l* < 20, but at shorter wavelengths of up to *l*_{max} = 70 it tends to fall below the line. This indicates sectoral anisotropy, which they attributed to the generally north–south trend of oceanic ridges and coastlines, and to the fact that Euler poles describing plate motions are clustered at high latitudes (Argus *et al.* 2011).

Well-sampled scalar fields on a sphere are independent of coordinates. Therefore, we are free to rotate the coordinate system to any arbitrary orientation. Although the individual coefficients in each orientation differ from one another, their synthesis reproduces the original model field. This is a well-known fact of differential geometry (Guggenheimer 1977). Furthermore, because the isotropy measure actually quantifies zone-to-sector isotropy, it can be used in any orientation, without modification.

We computed isotropy spectra for the geoid and topography at 42 orientations defined by the 12 vertices of a regular icosahedron and the midpoints of its 30 edges. Due to rotational symmetry, only 21 of these orientations are unique. As mentioned earlier, the formal resolution of these rotated fields was limited to *l*_{max} = 800 by available hardware memory. As a hedge against sampling artefacts, we retained all of the results.

Resolution of EGM2008 over Antarctica is limited to *l*_{max} = 180 due to the GRACE (Gravity Recovery and Climate Experiment)-only data available there (Pavlis *et al.* 2012). Furthermore, over large parts of Asia, Africa and South America, residual terrain modelling (RTM) was used to fill in areas where access to *in situ* gravity observations was proprietary. However, RTM was applied only beyond *l*_{max} = 900 (Pavlis *et al.* 2013) and should not affect our results.

Isotropy measures for topography and the geoid range from about 0.8 to 1.4 and higher, particularly at long wavelengths (Fig. 8). Recall that if the isotropy hypothesis was universally true, then all of these curves would plot at *R*(*l*) = 1. The fact that they do not indicates significant broadband anisotropy in Earth's structure.

The isotropy measure in the standard geographical orientation generally falls at the bottom of this range. In contrast, isotropy about the shape-perturbation minimum, at 2.0° N, 119.5° E, generally falls in the top of this range. Is this an artefact of Earth's dominantly ellipsoidal shape? If so, all five of the equatorial spectra should appear as a single curve. However, as indicated by the yellow curves in the geoid isotropy spectrum, this is not the case.

Thus, the zonal energy in Earth's plate system is maximized about the shape-perturbation minimum, particularly for wavelengths shorter than about 700 km (*l* > 58). This is a remarkably broadband indication of global structural symmetry about an axis extending from eastern Borneo to Brazil, which itself is defined by plate dynamics.

The isotropy measure for topography at short wavelengths *l* > 400 is maximized about 0° N, 342° E, in the Romanche Fracture Zone. This is practically coincident with Earth's least principal inertia axis (Chen & Shen 2010). The affected spectral range coincides with uncompensated topography, suggesting that these features are carried passively by the rigid lithosphere.

The equatorial positions of the eastern Borneo and Romanche Fracture Zone poles are strong indications for coupling between global plate motions and true polar wander. These findings complement work by Chao & Gross (1995) that shows feedback between earthquakes and rotational energy, and they complement work by Chao *et al.* (1995) that demonstrates feedback between earthquakes and the gravity field.

## Geological interpretation

### Crustal detachments

Metamorphic core complexes (MCCs) comprise mid- to lower-crustal footwall rocks that have been brought to the surface along crustal detachments. Detailed structural mapping and geochronological analysis have led to the recognition that crustal anatexis is commonly coeval with the extension and exhumation of these characteristic rock packages (Kruckenberg *et al.* 2008). MCCs are found in all mountain belts, but they are particularly well exposed in the Cordillera (Coney & Harms 1984; DeCelles 2004). The question is why?

The global-enhanced disturbance (Fig. 1f) shows five prominent orogenic belts: the Himalaya, the New Guinea Highlands, the Southern Alps, the Andes and the Cordillera. The first two are coincident with CC-type convergent margins, the third with an OO-type margins, the fourth with an OC-type margin, while the last parallels a major transform margin. Is this significant? Global plate dynamics says that it is.

As shown earlier, the Himalaya–Tibet region features a short-scale gravity-enhanced orogen nestled in a regional flexure of the continental lithosphere, overlying a former subduction zone. We infer this from the complementary nature of the compensated and enhanced gravity disturbances and their spatial coincidence with seismicity. The mountain load is isostatically compensated by the flexed lithosphere, which is about 255 km thick. As erosion proceeds, this prominent highland rises, and it probably will continue to do so for some time. This is a particularly spectacular example of a ‘collisional orogen’ (Cloos 1993).

We can interpret the complementary gravity disturbances in the Andes in the same way. However, the continental lithosphere there appears to be thinner, only about 173 km thick, based on the shorter wavelength of the regional compensated disturbance and the depth-distribution of earthquakes.

During the Mesozoic, the entire western margin of North America was a subduction zone consuming oceanic lithosphere of the Farallon Plate. In the late Paleocene–early Eocene, the Pacific Plate came into contact with the North American Plate, producing a transform margin. By the Miocene, subduction had ceased entirely, replaced by the San Andreas megashear (Atwater 1970).

Over this period, the dynamic boundary condition for the North American Cordillera changed from convergence and shortening to the present-day transtensional regime of the Great Basin. Thus, a gravity-enhanced orogenic welt that had been built up over some 40 myr began to collapse (Coney & Harms 1984; Constenius 1996; DeCelles 2004). Subsequent crustal extension and isostatic rebound have brought partially melted mid- to lower-crustal rocks to the surface.

Peak metamorphic pressures for MCCs in the Cordillera are consistent with depths predicted by ThERM for the mixed-mode spectral range *l* = 541*−*2900 (Table 2) (Patton & Watkinson 2013). Furthermore, equilibrium melting is likely given the reversion to thermodynamic symmetry at sub-H3 depths (see the subsection ‘Heat transfer’ in Appendix A). In light of this it seems unlikely that the association of these characteristic rock packages with crustal-scale shear zones – detachments – is mere coincidence. Rather, it seems likely to be a natural consequence of globally correlated shear-localization processes inferred from the mechanics of DG-2 materials.

### Subduction interface

The association of topography, seismicity and volcanism has long aroused geological curiosity (e.g. Dutton 1889). Today, it is recognized as evidence for plate convergence and subduction (Benioff 1954; Isacks *et al.* 1968). This process and its role in crustal differentiation have been studied extensively. But how do you melt rocks that are cooler than their surroundings?

The gravity disturbance in the Palau–Yap–Mariana Islands region of the SW Pacific exhibits three arcuate lows, which increase in size from west to east (Fig. 9). Palau lies at the western cusp, while the deepest point in the oceans, the Challenger Deep, lies about 3° east of the eastern cusp. Gravity disturbances in this region are negatively skewed, in contrast to the global gravity disturbance which is positively skewed by a factor of 2 (Fig. 1b).

Seismicity in Palau and the Yap Islands is relatively shallow, less than about 37 km. In contrast, seismicity beneath the Mariana Islands extends to depths in excess of 255 km, particularly in the north. The shallow seismicity of the Yap Islands bifurcates at the Challenger Deep, defining plate margins in the Mariana forearc and back-arc (Bird 2003). Volcanoes occur in the Mariana arc.

Thus, while subduction is well established in the Mariana arc, this does not appear to be the case in the Yap Trench. However, the similar arcuate map pattern there can be attributed to shortening in the crustal detachment zone, as discussed earlier. Both the depth of earthquakes and their mechanisms are consistent with this idea.

Hatherton & Dickinson (1969) found a relationship between the potassium content of andesites and the depth to earthquakes beneath island arcs. They reported that the range of depths above which andesitic rocks are found is 80–290 km. This is practically equivalent to the sub-L3 depths predicted for the mixed-mode spectral range *l* = 79*−*416 (Table 2), and is coincident with the low-velocity zone of PREM (Dziewonski & Anderson 1981) (Fig. 10).

The integrated pattern suggests some degree of self-similarity in the deformation mechanisms accommodating poloidal plate motions. One can anticipate shortening and thickening of the overriding layers. Another self-similar pattern, albeit expressed in ancient crustal thickness, was noted in the Yilgarn and Pilbara (Patton 2007). What is the connection to crustal differentiation?

Chao & Gross (1995) demonstrated a strong tendency for earthquakes to decrease global gravitational energy. Their estimates of energy feedback to Earth's rotation and seismic-wave energy release were three orders of magnitude smaller than the observed decrease. Assuming that elastic energy cannot build up continuously over time, they argued that earthquakes convert most of this energy to heat.

The mechanics of DG-2 materials shows that earthquakes are dynamically necessary (see the subsection ‘Seismogenesis’ in Appendix A), limiting the local build-up of strain energy. Furthermore, the thermodynamics of DG-2 materials can revert to classical energy–entropy relationships, manifested as Joule heating, for depths below L3 (Table 1) (Patton & Watkinson 2013) (see the subsection ‘Heat transfer’ in Appendix A). Thus, excess gravitational energy can be liberated as heat, which in turn causes melting. In light of this, the curious association of topography, seismicity and volcanism at convergent margins arises from the local superposition of time-like energetic feedbacks tied to the subduction interface.

### Shear-mediated phase changes

Seismic discontinuities in the mantle transition zone are widely attributed to phase changes in the olivine–wadsleyite–ringwoodite–perovskite system (Ringwood 1991). However, it is possible that these phase changes are shear-mediated.

For example, mineral physics experiments have shown that differential stresses affect the equilibrium pressure–temperature conditions for phase changes in quartz and enstatite (Coe & Paterson 1969; Coe & Kirby 1975). This would appear to result from the mechanical coupling of expansion to shear that is expressed in kinematic relationships (see equation A3 in Appendix A) and modelled by finite deformation theory (Truesdell & Noll 2004). A geologically relevant example of this coupling is found in DG-2 materials (Rivlin & Ericksen 1955).

Mixed-mode response in the spectral range *l* = 30*−*60 (Fig. 6) predicts that at depths of between 663 and 690 km deformation should be dominated by shear banding. Seismic discontinuities in the mantle transition zone have average depths of 415, 520 and 660 km. Of these, the 660 km-discontinuity exhibits the largest impedance contrast (Shearer 2000). Perhaps its global observability reflects the influence of localized shearing at that depth.

### Sunda-centric symmetry

Present-day subduction zones compose two distinct interconnected groups: one centred on the Sunda Plate and the other on the South America Plate. These groups lie at opposite ends of the axis defined by the minimum in Earth's dynamically induced shape perturbation (Fig. 5). Subduction zones tend to lie on either sectoral (radial) or zonal (concentric) arcs about this axis, as shown in Figure 11.

For example, the Middle America and Chile trenches lie on sectoral arcs about 120° from one another, while the Peru trench lies on a zonal arc about 15° from the Brazil antipole. These geometric observations are consistent with 3–3 sectoral symmetry, and 24–0 zonal symmetry, about the shape minimum pole.

Similarly, the tectonically complex margins extending from Borneo to Fiji and from Borneo to the Gulf of Alaska approximate sectoral arcs about 90° from one another. Furthermore, the Phillipine, Sunda and Tonga trenches approximate zonal arcs at distances of about 12°, 15° and 70° from the eastern Borneo pole. These geometrical observations are consistent with 2–2 sectoral symmetry, and with 30–0, 24–0 and 5–0 zonal symmetry, respectively, about the shape-minimum pole. The Alpine–Himalayan belt generally lies parallel to one of the 2–2 sectoral arcs.

Degree-equivalent radii of curvature for subduction zones/island arcs lie in the range 5–100 (Table 3). The Sandwich arc has the shortest radius at about 400 km (*l* = 100), while the Chile and Tonga arcs have the longest radii at about 8000–8500 km (*l* = 5). Most of the arcs have radii in the range 800–3000 km (12 < *l* < 51). These observations are consistent with subduction zones forming as a result of folding and shear banding in the upper mantle alone (i.e. above M4). Significantly, the latter range is where ThERM predicts layer CMB to respond rigidly (Fig. 4). The lower mantle need not be involved in global plate motions.

The zonal gravitational energy of Earth's present plate system, at wavelengths of less than about 700 km, is maximized around the shape-minimum pole (Fig. 8, *R*_{N}(*l*)). Note that the present motion of the Pacific Plate, in the Sunda Plate frame, is parallel to one of the 2–2 sectoral symmetry planes. It is also parallel to the trend of the Hawaiian Island chain. Given the known age progression of volcanic rocks comprising these islands, the dynamically induced symmetry of Earth's tectonic structure about the shape-minimum pole is likely to have influenced global plate motions for the last 45 myr.

## Global plate dynamics

Coupling of poloidal to toroidal motions in the self-gravitating Earth is essential for global plate dynamics. Pointing out that ‘there is no generally accepted definition of convection’, Elsasser (1967, p 227) defined convection as ‘motion which results from density differences in a gravitational field’. His ‘stress guide hypothesis’ invoked plastic deformation in the asthenosphere to explain plate mobility in large aspect ratio convection cells. Dense subducting slabs pulled the rest of the plate along, uncoupled from the mantle beneath.

McKenzie (1969, p 29) explored this hypothesis but could not reconcile it with known plate geometry. He found that ‘the only useful simplification of the convection problem appears to be the idea of plates’. Reluctantly, it seems, he came to the conclusion that there is ‘no obvious alternative to numerical solution of the general non-linear convection equations’. Reinforced by the effective viscosity of rock mechanics, geodynamics thus became synonymous with convection in viscous fluids.

Kaula (1980) and Bercovici *et al.* (2000) have presented evolution equations for poloidal–toroidal coupling in viscous fluids. They found that viscosity gradients were needed to achieve coupling. Incompressibilty was assumed so that the Helmholtz decomposition could be used to express the results in terms of separate poloidal and toroidal potentials. Coupling to gravity as a body force was through the usual Boussinesq approximation, where temperature-dependent volume changes perturb the pressure alone.

Numerical simulations based on these concepts demonstrate that radial viscosity gradients alone are not sufficient for coupling. This leads to the conclusion that lateral viscosity gradients are essential (Forte & Peltier 1994). This point of view is reinforced by the conviction that the dynamic mantle signal is hidden beneath geological noise (Kaula 1969; Panasyuk & Hager 2000; Flament *et al.* 2013).

Evolution equations for poloidal and toroidal motions in viscous fluids can be recorded completely in terms of the shearing and the expansion. In contrast, evolution equations for DG-2 materials require both the rotation and time-like motions, in addition to the other two kinematic elements. This allows DG-2 materials to manifest the geologically interesting behaviour discussed in this paper, which is simply not possible in viscous fluids.

Algebraically, the toroidal evolution (equation A19 in Appendix A) admits a potential ψ. In turn, at leading order in a small parameter, *ε*, this potential is governed by an inhomogeneous diharmonic equation (equation A26 in Appendix A), a balance equation for the rotation. This result is independent from the incompressibility assumption. The balance of angular momentum implied by the toroidal potential provides a boundary condition for the poloidal field.

A deeper understanding of these results and their application to global plate dynamics can be attained by considering three limiting kinematic cases: viz. incompressibility, symmetry and anti-symmetry. Global incompressibility requires that all vertical derivatives of the vertical velocity vanish. Practically, this means a stably stratified structure in which only toroidal motion may be present. Because these hypothetical motions are not coupled to expansion in any way, they will eventually decay to nothing, save for rotation of the body as a whole.

Global symmetry requires that all components of the rotation vanish. Practically, this means that no dislocation or localized shearing motions can be present. This simplifies the problem drastically, reducing it to the viscous fluid case (Kaula 1980; Bercovici *et al.* 2000). In turn, one would have to appeal to viscosity gradients in order to couple poloidal to toroidal motions.

Global anti-symmetry requires that all shearing motions vanish. Consequently, only dislocation motions are possible. A body rotating on its axis, exhibiting no other motions of any kind, is also possible. Each of these is a decidedly unnatural model for a dynamic planet, like Earth.

In summary, it appears that as a rocky self-gravitating body cools, its symmetrical modes of deformation begin to ‘freeze out’. A manifestation of this freezing is the emergence of finite strength and the ability to transmit shear waves. Frozen parts retain locally heterogeneous elastic fields, consistent with diharmonic scaling. The ultimate end state is not hydrostatic, but rather lithostatic. The final shape of a terrestrial body reflects its integrated history of deformation, providing ample grist for the geologist's mill.

## Conclusion

We began this contribution by asking four questions. We shall conclude by answering them in light of the insights offered by the mechanics and themodynamics of DG-2 materials. Our answers are informed not only by statistically measurable patterns in global geophysical data, but also by the characteristic structures found in naturally deformed tectonites:

Are extreme variations in effective viscosity necessary to drive global plate motions? Based on our study of poloidal–toroidal coupling in self-gravitating heat conducting DG-2 materials (see the subsection ‘Poloidal–toroidal coupling’ in Appendix A), the answer is categorically no. If normal-stress differences are included in geodynamic models, they render all material gradient terms, including viscosity gradients, superfluous at leading order. The non-linear coupling of shearing to all four kinematic elements (rotation, shearing, expansion and rupture) in these materials offers a rich physical theory for interpreting diverse geological phenomena. One can imagine numerous avenues for future research based on this work and its logical extensions to higher-order terms in suitably chosen small parameters. These same features are wholly absent in the mechanics of viscous fluids.

The failure of conventional geodynamics to explain global patterns in geophysical as well as geological observations, to the extent that DG-2 materials do, can be attributed to a wretched approximation, similar to the one that prompted Prandtl (

*c.*1904) to include viscosity in models of air flow over wings. At the time, appropriately scaled physical theory showed that the effects of viscosity occurred only in terms multiplied by a small parameter. Yet, if these terms were neglected, computations predicted there would be no drag on the wing, contrary to experience. Subsequently, Prandtl found that the viscosity terms dominate flows near the boundaries of solid objects – in boundary layers – and thus cannot be neglected.The strength of rocks has traditionally been discounted in geodynamics, and justifiably so. Rock mechanics experiments show that the influx of heat and fluids rapidly leads to irreversible deformation under load. Furthermore, the small deviation of Earth's topography from the reference ellipsoid testifies to the feeble strength of rock at the planetary scale.

On the other hand, it is not clear that dynamic similarity exists between the laboratory and planetary scales. After all, what is the planetary scale equivalent to the laboratory loading frame? More importantly, everyday experience shows that rock strength cannot be neglected. The structures observed and described in natural tectonites, by structural geologists the world over, are just that – rock-like. The lesson? A small cause can lead to a large effect, as shown in this contribution.

To what extent do Earth's topography and near-surface density variations mask the dynamics of the deep interior? Based on the fact that the only known deep dynamic signatures are the subducted slab flux, which is reasonably explained as perturbing Earth's overall ellipsoidal shape, and hotspots, which could be sourced from the upper mantle alone (see the subsection ‘Heat transfer’ in Appendix A), the answer appears to be not much. The mechanics of DG-2 materials suggest that the lower mantle is effectively rigid at all but the longest wavelengths. Furthermore, given the fact that the Prandtl number of geodynamics is about 10

^{22}– enormous – it is difficult to imagine how viscous fluid dynamics is relevant to the solid-state vertical transport of fertile rocks in the mantle. In contrast, the thermodynamics of self-gravitating heat conducting DG-2 materials offers a quantitative theory for these phenomena that can be tested, both observationally and experimentally.How is the azimuthal anisotropy of Earth's gravity and topography, first identified by Bills & Lemoine (1995), related to the dynamics of global plate motions? Based on our study of directional isotropy and its apparent symmetry about the eastern Borneo–Brazil shape-perturbation axis, we can say that tectonic structures at wavelengths of less than about 700 km directly reflect plate dynamics. This, by itself, is not surprising. What is surprising is the broadband nature of this mechanically induced pattern of folded and faulted rocks.

Is there evidence in global high-resolution gravity–topography spectra for competent layering in Earth's radial structure? Based on our normalized admittance study of modern high-resolution gravity and topography data, the answer is categorically yes. Earth's radial and lateral structure, apparent in global geophysical and geological data, looks remarkably like ThERM. This pertains to the interpretation of newly acquired data from broadband instruments deployed across the globe, as well as to observations from older studies of continental structure; all of which hint that geological structures exist throughout the upper mantle. One wonders what could be learned by applying the mechanics of DG-2 materials to other terrestrial bodies.

## Acknowledgements

This work is dedicated with great respect to John Ramsay, for his friendship, mentoring and insights. He set the course of a modern analytical approach to structural geology based on meticulous fieldwork and observation. Our paper is an attempt to integrate that structural geology into a data-driven geophysical, global model application. The authors gratefully acknowledge comments and suggestions from Bruce E. Hobbs and an anonymous reviewer which helped to improve and clarify the presentation. All work was completed on Linux using open-source code, particularly SHTOOLS, The Generic Mapping Tools, LaTeX, and JabRef. To mom and dad – if these were silent, the rocks would cry out.

## Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

## Appendix A: Mechanics of DG-2 materials

Analysis of the mechanics of DG-2 materials (Rivlin & Ericksen 1955) in four-dimensional spacetime (Misner *et al.* 1973) leads naturally to a theory for the thermomechanical evolution of matter energy in a self-gravitating spheroidal body. Structural predictions from this work are used to interpret the spectral and spatial characteristics of Earth's gravity, topography and seismicity, as presented in the main text. For convenience, the mathematical symbols used, their meanings and their physical dimensions are presented in Table A1.

#### Spacetime

We model spacetime as a four-dimensional manifold endowed with a metric *ds*^{2} of signature (*−*, +, +, +). Consequently, time-like vectors have negative length, while space-like vectors have positive length. For convenience in applications, the metric can be decomposed into spatial and temporal parts using the projection tensor *P _{αβ}* =

*g*+

_{αβ}*u*, where

_{α}u_{β}*u*are components of the velocity, viz.:

_{α}Greek indices run from 0 to 3, corresponding to one time and three space dimensions, respectively. Repeated indices imply summation.

#### Kinematics in spacetime

An understanding of complex mechanical systems begins with kinematics. For the study of relative motion between material particles within deforming solid bodies, a description based on velocity gradients is superior to other simpler, more intuitive approaches.

In spacetime, there are 16 components of the velocity gradients tensor *u _{α}*

_{;β}.

*A priori*, this tensor exhibits no particular symmetry. Any such tensor can be decomposed uniquely into anti-symmetrical

*Y*and symmetrical

*S*parts, so that

*u*

_{α}_{;β}=

*Y*+

_{αβ}*S*. Furthermore, it can be shown that the anti-symmetrical part is trace-free (i.e. free from changes in volume). The symmetrical part, in turn, can be uniquely decomposed into two parts: the first of which is trace-free,

_{αβ}*F*, and the other accounting for changes in volume (1/3)

*θP*. Thus,

_{αβ}*u*

_{α}_{;β}=

*Y*+

_{αβ}*F*+ (1/3)

_{αβ}*θP*.

_{αβ}Rotational motion is anti-symmetrical, given by the rotation two-form:

Rotation is essential for dislocation motion, and has interpretive value for dynamic rupture processes and seismogenesis. The form of the shearing tensor shows that shearing and expansion are coupled additively. The significance of these relationships for the interpretation of solid-state deformation cannot be overstated.

Separating spatial from temporal parts, the velocity gradients can be written as:

Here, the rotation, shearing and expansion are space-like elements of motion, while the last term is time-like. The overdot denotes a proper derivative.

The distinction between active and preserved deformation depends on whether the observed motions are time-like or space-like. Seismic ruptures are a particularly violent example of the former, while static fault offsets and displacement gradients (strain) are examples of the latter. The spacetime point of view is natural for describing the kinematics of a complex macroscopic system.

It is important to recognize that kinematic restrictions imposed either mathematically (e.g. incompressibility) or experimentally (e.g. simple shear) limit the potential response of a material to loading. Of course, such restrictions can render tractable an otherwise difficult theoretical problem and therefore yield insight, but in the laboratory they tend to obscure the phenomena of interest. In the field of rock mechanics, this problem was pointed out by Hobbs (1972).

#### Dynamics of DG-2 materials

We have attempted to circumvent the restrictions of the laboratory by studying finite deformation theoretically. We have invested significant effort in understanding DG-2 materials (Patton & Watkinson 2005, 2010, 2013) specifically because they manifest rock-like characteristic modes of deformation, which are absent in other simpler rheologies. Here, we extend our earlier work to include self-gravitation and heat conduction.

Einstein (1916) accounted for gravitation using the field equation:
*G* takes its source in a finite stress-energy density *T*. By construction, the mean curvature satisfies the contracted Bianchi identity. Consequently, its divergence must vanish (∇· *G* ≡ 0). In other words, the geometry of spacetime is conserved. This, in turn, constrains the evolution of stress energy to the equations of motion contained in ∇· *T* = 0.

Here, we assume that stress-energy density is given by:

Equation (A8) is the constitutive equation for a compressible DG-2 material (Rivlin & Ericksen 1955), where *A*_{(1)} is the velocity-strain tensor, defined by:
*A*_{(2)} is the acceleration-strain tensor, defined by:

Normal stress effects proportional to *β*_{1} and *β*_{2} are expected to be of the same order of magnitude (Dunn & Rajagopal 1995). The former modulus is associated with finite material strength (Patton & Watkinson 2005, 2010, 2013), while the latter is associated with non-linear viscosity. In another context, these moduli are associated with the absorption and phase shift of plane shear waves (Truesdell 1964).

The gross spheroidal symmetry of self-gravitating bodies, like Earth, is neatly explained by exact solutions of equation (A5) (Boyer & Lindquist 1967), which build on the seminal solution by Schwarzschild (*c.* 1916). Thus, the following developments govern local interactions of matter energy within a self-gravitating body of decidedly terrestrial nature.

#### Equations of motion

Upon taking the divergence of the stress-energy density (equation A6) and projecting it into the rest frame of an observer with velocity *u ^{α}*, we obtain the energy relationships:

The former relationship is the continuity equation for mass energy, while the latter, combined with the definition of the heat flux vector, is a convection–diffusion equation for heat.

Similarly, contracting the divergence of the stress energy with the projection tensor, we obtain:

The first relationship is a balance equation for interactions between heat flux and deformation, while the latter is the usual balance equation for linear momentum.

#### Heat transfer

Although it is well known that temperatures tend to increase with depth in the Earth, we require a specific model for such variation.

In his study of the thermodynamics of inhomogeneous systems of a single variable, Lavenda (1995) showed that a self-gravitating body, like a planet or star, is subject to an unconventional breaking of the symmetry expected between the entropy and energy descriptions. Thermodynamic symmetry is associated with systems at equilibrium. Consequently, he found that the energy density for such a body is given by a monotonically decreasing convex function of the entropy density. Significantly, the energy density is also inversely proportional to the body's radius, being high in the centre and decreasing outwards.

Earlier, Lavenda (1995) had shown that similar symmetry-breaking applied in the case of a stretched elastic band. Patton & Watkinson (2013) applied this argument to the thermodynamics of strained elastic solids, and thereby showed that when the energy density is a monotonically decreasing convex function of the entropy density, it is also proportional to length. This is consistent with first-order expectations for strain energy in rocks – *Ut tensio, sic vis* (Hooke *c.* 1678).

In each of these cases, the temperature of the system is proportional to the slope of the energy density function. For a self-gravitating body, both pressure and temperature are predicted to increase with depth (Fig. 10, *δU*_{D}). Consequently, heat will tend to flow radially, against the temperature gradient. For the strained-solid case, heat will be absorbed locally by the material and only be released upon encountering a phase change of some kind. Two kinds of phase changes are possible.

The first is an equilibrium phase change leading to conventional Joule heating (Chandrasekhar 1967). Because equilibrium implies a reversion to thermodynamic symmetry, we expect that heating of this kind will be associated with the dashed branch of the localization curve *δU*_{L} shown in Figure 10. Consequently, we predict that partial melting will occur at or below the H3, L3 and M3 modes of ThERM (Patton & Watkinson 2010), at depths of 13.8, 96 and 663 km, respectively. Timescales for conductive heating from these depths are about 6, 290 and 14 000 myr, respectively. Clearly, heat from the first two of these hypothetical sources could be expressed in the 4500 myr of Earth's history, while associated buoyant melts could be erupted on even shorter timescales.

The second is a non-equilibrium phase change which, based upon the dynamic rescaling of stress in DG-2 materials, leads to spontaneous shear rupture. The pulverization and comminution of particles in the process zone of dynamic ruptures can therefore be expected to evolve heat. The fact that pseudotachylite is often found associated with fault surfaces or injected into fault walls is consistent with these ideas.

Thus, we anticipate that heat will be evolved locally within the Earth system, in response to fault ruptures, as well as regionally, in response to tectonic forces (Chao & Gross 1995; Chao *et al.* 1995). These tectonic contributions to heat flow supplement those expected from radioactive decay and secular cooling.

Neglecting material gradient terms in equations (A7) and (A12), the constants can be collected into the thermal diffusivity *κ* = *k*/*ρC*. Straightforward non-dimensionalization yields a convenient reference thermal timescale for solid-state deformation problems, viz. *t*_{thermal} = *d*^{2}/*κ*.

#### Quasistatic momentum balance

Now we focus on equations (A8) and (A14). By introducing scales for length, time and stress as *d*, *d*^{2}/*κ*, and *ηκ*/*d*^{2}, respectively, and multiplying through by *d*^{3}/*ηκ*, we obtain:

This equation features three dimensionless numbers, set off in square brackets from the (now) dimensionless variables. To motivate further exploration of this equation in the context of the solid-state deformation of terrestrial planets, we present a simple physics-of-length-scales argument.

The first dimensionless number, *κ*/*ν*, is the reciprocal of the Prandtl number. Its square root gives the ratio of thicknesses of thermal to viscous boundary layers in the model (Tozer 1973). Given that the Prandtl number for Earth's mantle is universally believed to be enormous, any thermal boundary layer must be much thinner than a viscous one. This justifies neglect of the first (inertial) term of equation (A15) in the following developments.

The second dimensionless number, *κ*/*χ*, is ‘thermomechanical competence’ (Patton & Watkinson 2010). Its square root gives the ratio of thicknesses of thermal to mechanical boundary layers in the model. *A priori*, the size of this number is unknown. The third dimensionless number is the ratio of the normal stress moduli *ζ* = *β*_{2}/*β*_{1}. The square root of the product *ζκ*/*χ* therefore gives the ratio of thicknesses of thermal to shear attenuation boundary layers in the model. *A priori*, the size of this number is unknown. We can therefore conclude that convection–diffusion models based on equation (A15) will exhibit some behaviours altogether absent in the classical Rayleigh–Bénard problem.

#### Poloidal–toroidal coupling

An evolution equation is a differential law governing the temporal development of a system. This subsection outlines the derivation of such equations for the poloidal and toroidal velocity fields of a DG-2 material exhibiting a very large Prandtl number. These relationships can be extracted from the quasistatic balance of linear momentum, discussed above, by projecting successive curls of the divergence of stress onto a radial unit vector

Residual terms in these expressions are potential drivers for the targeted field. Results for viscous fluids (Kaula 1980; Bercovici *et al.* 2000), defined by equation (A8) with *β*_{1} = *β*_{2} = 0, are recovered as special cases. The results reported here for DG-2 materials exhibit explicit dependence on the natural time *β*_{1}/*η* (Truesdell 1964), even in the quasi-static limit. This is important for the nucleation and propagation of shear cracks, as is shown in the subsection on ‘Seismogenesis’ later in this Appendix.

#### Poloidal evolution

The evolution equation for the poloidal field can be found by substituting equation (A8) into the relationship:

The explicit result is rather lengthy and therefore not reproduced here. Implicitly, this relationship is a fifth-order partial differential equation of the vertical velocity component *w* alone, of the form:
_{h} is the two-dimensional horizontal Laplacian. Residual forcing terms are non-linear, quadratic in derivatives of the velocity field, and involve first-, second- and third-order derivatives of the viscosity and normal-stress moduli. For reasons presented below in connection with the toroidal field, these material gradient terms are suppressed here. The remaining terms of the function Π, proportional to *β*_{1}, can be organized by the velocity gradient components and the next higher spatial derivatives of them. Note the presence of the expansion *θ* and the absence of rotation *ω*. Very little algebraic collapse occurs in this expression, and what does can be summarized as follows.

Each of the velocity gradient components appears as a cofactor for a sum of fourth-order mixed partial derivatives of the velocity components. An obvious pattern is factors of the form *−*ΔΔ_{h}*u*, *−*ΔΔ_{h}*v*, and *−*ΔΔ_{h}*w* multiplying spatial derivatives of the vertical velocity *∂w*/*∂x*, *∂w*/*∂y* and *∂w*/*∂z*, respectively. These terms are associated with the convergence and divergence of material at horizontal boundaries.

Similarly, each of the second-order velocity derivatives appears as a cofactor for sums of mixed third-order partial derivatives of the velocity components. Of these, only ∂^{2}*u*/∂*z*^{2} and ∂^{2}*v*/∂*z*^{2} multiply factors of the form ∂^{2}*θ*/∂*x*∂*z* and ∂^{2}*θ*/∂*y*∂*z*. These terms, involving the expansion, provide coupling of poloidal motions to thermal and chemical density changes, and thereby give rise to non-isotropic pressure gradients in the mantle, which, in turn, drive such motion.

#### Toroidal evolution

The evolution equation for the toroidal field can be obtained by substituting equation (A8) into the relationship:

Again, the complete result is rather lengthy and therefore omitted here. First- and second-order derivatives of the viscosity and normal stress coefficients appear explicitly, but only the standard gradients are exhibited here in the general relationship:

While this is a fifth-order partial differential equation of the toroidal potential ψ, it is only a fourth-order equation in the velocity components. The residual forcing terms Λ are non-linear, quadratic in derivatives of the velocity components. Note the presence of all kinematic elements.

In the absence of normal stress differences, equation (A19) reduces to:

While all of the shearing components, except for *σ*_{33}, appear in this expression, the rotation components implicit in equation (A19) do not. Additionally, if the viscosity gradients were taken to vanish, then the right-hand side of equation (A20) would also vanish. Thus, in conventional mantle convection models, lateral viscosity gradients are needed to drive toroidal motions (Kaula 1980; Bercovici *et al.* 2000).

Further analysis, however, shows this to be neither necessary nor sufficient in the presence of finite normal-stress differences. Upon neglecting all material gradient terms, equation (A19) reduces to:

Note that rotation terms are retained, consistent with the documented dynamic connection between normal-stress differences and rotation (Patton & Watkinson 2010). Therefore, residual terms proportional to *β*_{1} can drive toroidal motions. Viscosity gradients are not required. Moreover, it appears that material gradients of any kind are not required. This suggests that a useful simplification of this highly non-linear system can be achieved by neglecting these terms. We adopt this geodynamic approximation in the following developments and applications.

In any case, material gradients lack intrinsic physical meaning, except in relation to a particular configuration (Truesdell & Noll 2004). This raises the issue of what defines a meaningful reference configuration. Material gradients clearly play a role in engineered systems and in the detailed structure of natural tectonites. Mathematically, these terms compose the majority of the terms found in the preceding analysis. What, then, is their purpose? Speculatively, they serve to ‘record’ the structure of complex mechanical systems.

#### Leading-order behaviour

Given that equation (A21) governs the toroidal evolution of the solid physical surface of a thermomechanical matter-energy configuration – a terrestrial planet – we can use perturbations about a pure-shearing plane-strain base velocity field to linearize this expression. This yields the somewhat simpler partial differential equation, in dimensionless variables, viz.:

Upon dividing through by 1 + 2*κ*/*χ* and defining the function:

Factoring the last three terms on the left-hand side, taking advantage of the rescaling operator defined in Patton & Watkinson (2010) as:
*ω*, equation (A24) can be written as:

This is a third-order partial differential equation where the first term is a proper derivative of the diffusion of rotation and the second term is diharmonic.

Thus, at steady-state, in a system possessing a large Prandtl number, we find that the leading-order behaviour of equation (A21) reduces to a homogeneous diharmonic equation (Patton & Watkinson 2005, 2010, 2013), viz.:

Under the arguably isothermal conditions pertaining to the Earth's surface, this result suggests that, at any moment in time, toroidal motions should manifest localized shear-band structures, as predicted by the dynamic rescaling theorem of DG-2 materials and its corollary (Patton & Watkinson 2010). The implications of this result for global plate dynamics and for the analysis of gravity–topography cross-spectra require an appreciation of its connections to the dynamics of crack propagation.

#### Seismogenesis

Equation (A27) can be interpreted as a vector wave equation, in wavefront-normalized coordinates, for shearing motions occurring in the near field of a propagating in-plane crack, provided that the following relationship holds, viz.:

Here, *v*_{c} is the rupture speed of the crack, *μ* and *ρ* are the shear rigidity and density, respectively. Note that compressional waves are also generated by an in-plane dislocation (Aki & Richards 2002).

Together, equations (A23) and (A28) predict that crack-rupture speed depends on the thermomechanical competence *κ*/*χ* according to:

Thus, as thermomechanical competence increases, rupture speed asymptotically approaches the Eshelby speed, *V*_{S} at *κ*/*χ* = 1/2 (Fig. 3). Consequently, DG-2 dynamics predicts supershear rupture speeds for *κ*/*χ* > 1/2. This supershear domain coincides with that for which shear banding is possible, as summarized in the following subsection on ‘Diharmonic structural theory’. Furthermore, the onset of shear strain localization is spontaneous, as predicted by the toroidal evolution equation (equation A26). In this light, at plate margins or other highly stressed regions, earthquakes appear to be a dynamic necessity.

In order to eliminate the stress and velocity singularities present in homogeneous dynamic rupture models, Barenblatt (1959) introduced an anelastic cohesive process zone of width *d* at the tip of a propagating rupture. Comparing that work with Patton & Watkinson (2010), we find the relationship:
*τ*_{c} is cohesive stress. Later, Ida (1972) introduced slip-rate dependent cohesion, with critical slip distance *D*. Comparing that work with Patton & Watkinson (2010), we find the relationship:

These two relationships are the same when the critical slip distance is given by *D* = 2Γ/*τ*_{c}. Equations (A30) and (A31) suggest that laboratory-scale rock mechanics experiments could be used to estimate the rescaling modulus *|α|* found in the theory of DG-2 materials (Patton & Watkinson 2005). Note that dynamic rescaling is exact in the rotation *ω* (Patton & Watkinson 2010).

#### Diharmonic structural theory

DG-2 materials manifest fold and shear-band structures as a function of thermomechanical competence (Patton & Watkinson 2005, 2010, 2013). Folds are possible for all values of *κ*/*χ*, but only for *κ*/*χ* < 1/2 are the four distinct roots of equation (A27) purely real. This range corresponds to material folds, with normalized wavelengths *ξ* = *λ*/*z* > 2. Folds with even longer normalized wavelengths correspond to lower values of thermomechanical competence, with a transition at *κ*/*χ* ≥ 1/2. These shear bands form symmetrically about the shortening axis, at angles ranging from 0° up to 45°, as *κ*/*χ* increases.

For *κ*/*χ* ≥ 1, DG-2 materials are isothermal and effectively rigid, and in this range admit complementary shear bands under simple-shearing plane strain. With respect to the shearing direction, angles for primary shear bands range from 45° down to 0°, as *κ*/*χ* increases. Secondary shear bands form at right angles to the primaries.

Thus, the dynamic rescaling of stress and thermodynamics of DG-2 materials predicts a range of mechanical responses as a function of normalized wavelength: viz. homogeneous deformation, folding, banding and rigidity (Fig. 2).Throughout this contribution, we use these convenient word labels to refer to relative ranges of *κ*/*χ* – thermomechanical competence. If semantic questions arise, the mathematical relationships that underpin these concepts and the context of their usage should be the final arbiter. The main article applies these relationship to the interpretation of Earth's gravity–topography cross-spectrum and its correlation with seismicity.

- © 2019 The Author(s). Published by The Geological Society of London. All rights reserved