Frederik J. Simons

Probabilistic assessment of sea level during the last interglacial stage

With polar temperatures ∼3–5 °C warmer than today, the last interglacial stage (∼125 kyr ago) serves as a partial analogue for 1–2 °C global warming scenarios. Geological records from several sites indicate that local sea levels during the last interglacial were higher than today, but because local sea levels differ from global sea level, accurately reconstructing past global sea level requires an integrated analysis of globally distributed data sets. Here we present an extensive compilation of local sea level indicators and a statistical approach for estimating global sea level, local sea levels, ice sheet volumes and their associated uncertainties. We find a 95% probability that global sea level peaked at least 6.6 m higher than today during the last interglacial; it is likely (67% probability) to have exceeded 8.0 m but is unlikely (33% probability) to have exceeded 9.4 m. When global sea level was close to its current level (≥-10 m), the millennial average rate of global sea level rise is very likely to have exceeded 5.6 m kyr-1 but is unlikely to have exceeded 9.2 m kyr-1. Our analysis extends previous last interglacial sea level studies by integrating literature observations within a probabilistic framework that accounts for the physics of sea level change. The results highlight the long-term vulnerability of ice sheets to even relatively low levels of sustained global warming.The Last Interglacial (LIG) stage (ca. 130 115 ka), with polar temperatures likely 3 5◦C warmer than today, serves as a partial analogue for low-end future warming scenarios. Multiple indicators suggest that LIG global sea level (GSL) was higher than at present; based upon a small set of local sea level indicators, the Intergovernmental Panel on Climate Change (IPCC)s Fourth Assessment Report inferred an elevation of approximately 4 6 m. While this estimate may be correct, it is based upon overly simplistic assumptions about the relationship between local sea level and global sea level. Sea level is often viewed as a simple function of changing global ice volume. This perspective neglects local variability, which arises from several factors, including the distortion of the geoid and the elastic and isostatic deformation of the solid Earth by shifting ice masses. Accurate reconstruction of past global and local sea levels, as well as ice sheet volumes, therefore requires integrating globally distributed data sets of local sea level indicators. To assess the robustness of the IPCCs global estimate and search for patterns in local sea level that are diagnostic of meltwater sources, we have compiled a comprehensive database that includes a variety of local sea level indicators from 47 localities, as well as a global sea level record derived from oxygen isotopes. We generate a global synthesis from these data using a novel statistical approach that couples Gaussian process regression to Markov Chain Monte Carlo simulation of geochronological errors. Our analysis strongly supports the hypothesis that global sea level during the Last Interglacial was higher than today, probably peaking between 6 9 m above the present level. This level is close to that expected from the complete melting of the Greenland Ice Sheet, or from major melting of both the Greenland and West Antarctic Ice Sheets. In the period when sea level was within 10 m of the modern value, the fastest rate of sea level rise sustained for a 1 ky period was likely about 80 110 cm per century. Combined with the evidence for mildly higher temperatures during the LIG, our results highlight the vulnerability of ice sheets to even relatively low levels of sustained global warming.

Spatiospectral Concentration on a Sphere

We pose and solve the analogue of Slepians time-frequency concentration problem on the surface of the unit sphere to determine an orthogonal family of strictly bandlimited functions that are optimally concentrated within a closed region of the sphere or, alternatively, of strictly spacelimited functions that are optimally concentrated in the spherical harmonic domain. Such a basis of simultaneously spatially and spectrally concentrated functions should be a useful data analysis and representation tool in a variety of geophysical and planetary applications, as well as in medical imaging, computer science, cosmology, and numerical analysis. The spherical Slepian functions can be found by solving either an algebraic eigenvalue problem in the spectral domain or a Fredholm integral equation in the spatial domain. The associated eigenvalues are a measure of the spatiospectral concentration. When the concentration region is an axisymmetric polar cap, the spatiospectral projection operator commutes with a Sturm--Liouville operator; this enables the eigenfunctions to be computed extremely accurately and efficiently, even when their area-bandwidth product, or Shannon number, is large. In the asymptotic limit of a small spatial region and a large spherical harmonic bandwidth, the spherical concentration problem reduces to its planar equivalent, which exhibits self-similarity when the Shannon number is kept invariant.

The deep structure of the Australian continent from surface wave tomography

Abstract We present a new model of 3-D variations of shear wave speed in the Australian upper mantle, obtained from the dispersion of fundamental and higher-mode surface waves. We used nearly 1600 Rayleigh wave data from the portable arrays of the Skippy project and from permanent stations (from Agso , Iris and Geoscope ). Agso data have not been used before and provide better data coverage of the Archean cratons in western Australia. Compared to previous studies we improved the vertical parameterization, the weighting scheme that accounts for variations in data quality and reduced the influence of epicenter mislocation on velocity structure. The dense sampling by seismic waves provides for unprecedented resolution of continental structure, but the wave speed beneath westernmost Australia is not well constrained. Global compilations of geological and seismological data (using regionalizations based on tectonic behavior or crustal age) suggest a correlation between crustal age and the thickness and composition of the continental lithosphere. However, the age and the tectonic history of crustal elements vary on wavelengths much smaller than have been resolved with global seismological studies. Using our regional upper mantle model we investigate how the seismic signature of tectonic units changes with increasing depth. At large wavelengths, and to a depth of about 200 km, the inferred velocity anomalies corroborate the global pattern and display a progression of wave speed with crustal age: slow wave propagation prevails beneath the Paleozoic fold belts in eastern Australia and wave speeds increase westward across the Proterozoic and reach a maximum in the Archean cratons. The high wave speeds associated with Precambrian shields extend beyond the Tasman Line, which marks the eastern limit of Proterozoic outcrop. This suggests that parts of the Paleozoic fold belts are underlain by Proterozoic lithosphere. We also infer that the North Australia craton extends off-shore into Papua New Guinea and beneath the Indian Ocean. For depths in excess of 200 km a regionalization with smaller units reveals that some tectonic subregions of Proterozoic age are marked by pronounced velocity highs to depths exceeding 300 km, but others do not and, surprisingly, the Archean units do not seem to be marked by such a thick high wave speed structure either. The Precambrian cratons that lack a thick high wave speed “keel” are located near passive margins, suggesting that convective processes associated with continental break-up may have destroyed a once present tectosphere. Our study suggests that deep lithospheric structure varies as much within domains of similar crustal age as between units of different ages, which hampers attempts to find a unifying relationship between seismic signature and lithospheric age.

Isostatic response of the Australian lithosphere: Estimation of effective elastic thickness and anisotropy using multitaper spectral analysis

Gravity and topography provide important insights regarding the degree and mechanisms of isostatic compensation. The azimuthally isotropic coherence function be- tween the Bouguer gravity anomaly and topography evolves from high to low for increasing wavenumber, a diagnostic that can be predicted for a variety of lithospheric loading models and used in inversions for flexural rigidity thereof. In this study we investigate the isostatic response of continental Australia. We consider the effects of directionally anisotropic plate strength on the coherelce. The anisotropic coherence function is calculated for regions of Australia that have distinctive geological and geophysical properties. The coherence estimation is performed by the Thomson multiple-Slepian-taper spectral analysis method extended to two-dimensional fields. Our analysis reveals the existence of flexural anisotropy in central Australia, indicative of a weaker N-S direction of lower Te. This observation is consistent with the suggestion that the parallel faults in that area act to make the lithosphere weaker in the direction perpendicular to them. It can. also be related to the N-S direction of maximum stress and possibly the presence of E-W running zones weakened due to differential sediment burial rates. We also demonstrate that the multitaper method has distinct advantages for computing the isotropic coherence function. The ability to make many independent estimates of the isostatic response that are minimally affected by spectral leakage results in a coherence that is more robust than with modified periodogram methods, particularly at low wavenumbers. Our analysis elucidates the reasons for discrepancies in previous estimates of effective elastic thickness Te of the Australian lithosphere. In isotropic inversions for Te, we obtain values that are as much as a factor of 2 less than those obtained in standard inversions of the periodogram coherence using Bouguer gravity and topography but greater than those obtained by inversions that utilize free-air rather than Bouguer gravity and ignore the presence of subsurface loads. However, owing to the low spectral power of the Australian topography, the uncertainty on any estimate of Te is substantial.

Spherical Slepian functions and the polar gap in geodesy

The estimation of potential fields such as the gravitational or magnetic potential at the surface of a spherical planet from noisy observations taken at an altitude over an incomplete portion of the globe is a classic example of an ill-posed inverse problem. We show that this potential-field estimation problem has deep-seated connections to Slepians spatiospectral localization problem which seeks bandlimited spherical functions whose energy is optimally concentrated in some closed portion of the unit sphere. This allows us to formulate an alternative solution to the traditional damped least-squares spherical harmonic approach in geodesy, whereby the source field is now expanded in a truncated Slepian function basis set. We discuss the relative performance of both methods with regard to standard statistical measures such as bias, variance and mean squared error, and pay special attention to the algorithmic efficiency of computing the Slepian functions on the region complementary to the axisymmetric polar gap characteristic of satellite surveys. The ease, speed, and accuracy of our method make the use of spherical Slepian functions in earth and planetary geodesy practical.

Spatiospectral localization of isostatic coherence anisotropy in Australia and its relation to seismic anisotropy: Implications for lithospheric deformation

We investigate the two-dimensional (2-D) nature of the coherence between Bouguer gravity anomalies and topography on the Australian continent. The coherence function or isostatic response is commonly assumed to be isotropic. However, the fossilized strain field recorded by gravity anomalies and their relation to topography is manifest in a degree of isostatic compensation or coherence which does depend on direction. We have developed a method that enables a robust and unbiased estimation of spatially, directionally, and wavelength-dependent coherence functions between two 2-D fields in a computationally efficient way. Our new multispectrogram method uses orthonormalized Hermite functions as data tapers, which are optimal for spectral localization of nonstationary, spatially dependent processes, and do not require solving an eigenvalue problem. We discuss the properties and advantages of this method with respect to other techniques. We identify regions on the continent marked by preferential directions of isostatic compensation in two wavelength regimes. With few exceptions, the short-wavelength coherence anisotropy is nearly perpendicular to the major trends of the suture zones between stable continental domains, supporting the geological observation that such zones are mechanically weak. Mechanical anisotropy reflects lithospheric strain accumulation, and its presence must be related to the deformational processes affecting the lithosphere integrated over time. Three-dimensional models of seismic anisotropy obtained from surface wave inversions provide an independent estimate of the lithospheric fossil strain field, and simple models have been proposed to relate seismic anisotropy to continental deformation. We compare our measurements of mechanical anisotropy with our own model of the azimuthally anisotropic seismic wave speed structure of the Australian lithosphere. The correlation of isostatic anisotropy with directions of fast wave propagation gleaned from the azimuthal anisotropy of surface waves decays with depth. This may support claims that above ?200 km, internally coherent deformation of the entire lithosphere is responsible for the anisotropy present in surface wave speeds or split shear waves.

Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere

Abstract We interpret the three-dimensional seismic wave-speed structure of the Australian upper mantle by comparing its azimuthal anisotropy to estimates of past and present lithospheric deformation. We infer the fossil strain field from the orientation of gravity anomalies relative to topography, bypassing the need to extrapolate crustal measures, and derive the current direction of mantle deformation from present-day plate motion. Our observations provide the depth resolution necessary to distinguish fossil from contemporaneous deformation. The distribution of azimuthal seismic anisotropy is determined from multi-mode surface-wave propagation. Mechanical anisotropy, or the directional variation of isostatic compensation, is a proxy for the fossil strain field and is derived from a spectral coherence analysis of digital gravity and topography data in two wavelength bands. The joint interpretation of seismic and tectonic data resolves a rheological transition in the Australian upper mantle. At depths shallower than ∼150–200 km strong seismic anisotropy forms complex patterns. In this regime the seismic fast axes are at large angles to the directions of principal shortening, defining a mechanically coupled crust–mantle lid deformed by orogenic processes dominated by transpression. Here, seismic anisotropy may be considered ‘frozen’, which suggests that past deformation has left a coherent imprint on much of the lithospheric depth profile. The azimuthal seismic anisotropy below ∼200 km is weaker and preferentially aligned with the direction of the rapid motion of the Indo-Australian plate. The alignment of the fast axes with the direction of present-day absolute plate motion is indicative of deformation by simple shear of a dry olivine mantle. Motion expressed in the hot-spot reference frame matches the seismic observations better than the no-net-rotation reference frame. Thus, seismic anisotropy supports the notion that the hot-spot reference frame is the most physically reasonable. Independently from plate motion models, seismic anisotropy can be used to derive a best-fitting direction of overall mantle shear.

Minimum-Variance Multitaper Spectral Estimation on the Sphere

We develop a method to estimate the power spectrum of a stochastic process on the sphere from data of limited geographical coverage. Our approach can be interpreted either as estimating the global power spectrum of a stationary process when only a portion of the data are available for analysis, or estimating the power spectrum from local data under the assumption that the data are locally stationary in a specified region. Restricting a global function to a spatial subdomain—whether by necessity or by design—is a windowing operation, and an equation like a convolution in the spectral domain relates the expected value of the windowed power spectrum to the underlying global power spectrum and the known power spectrum of the localization window. The best windows for the purpose of localized spectral analysis have their energy concentrated in the region of interest while possessing the smallest effective bandwidth as possible. Solving an optimization problem in the sense of Slepian (1960) yields a family of orthogonal windows of diminishing spatiospectral localization, the best concentrated of which we propose to use to form a weighted multitaper spectrum estimate in the sense of Thomson (1982). Such an estimate is both more representative of the target region and reduces the estimation variance when compared to estimates formed by any single bandlimited window. We describe how the weights applied to the individual spectral estimates in forming the multitaper estimate can be chosen such that the variance of the estimate is minimized.

Spectral estimation on a sphere in geophysics and cosmology

We address the problem of estimating the spherical-harmonic power spectrum Sl of a statistically isotropic scalar signal s(r) from noise-contaminated data d(r) = s(r) + n(r) on a region R of the unit sphere. Three different methods of spectral estimation are considered: (i) the spherical analogue of the 1-D periodogram, (ii) the maximum likelihood method, and (iii) a spherical analogue of the 1-D multitaper method. The periodogram exhibits strong spectral leakage, especially for small regions of area A << 4π, and is generally unsuitable for spherical spectral analysis applications, just as it is in 1-D. The maximum likelihood method is particularly useful in the case of nearly whole-sphere coverage, A=4π, and has been widely used in cosmology to estimate the spectrum of the cosmic microwave background radiation from spacecraft observations. The spherical multitaper method affords easy control over the fundamental tradeoff between spectral resolution and variance, and is easily implemented, requiring neither non-linear iteration nor large-scale matrix inversion. As a result, the method is ideally suited for routine applications in geophysics, geodesy or planetary science, where the objective is to obtain a spatially localized estimate of the spectrum of a signal s(r) from data d(r) = s(r)+n(r) within a pre-selected and typically small region R.

Mapping Greenland’s mass loss in space and time

The melting of polar ice sheets is a major contributor to global sea-level rise. Early estimates of the mass lost from the Greenland ice cap, based on satellite gravity data collected by the Gravity Recovery and Climate Experiment, have widely varied. Although the continentally and decadally averaged estimated trends have now more or less converged, to this date, there has been little clarity on the detailed spatial distribution of Greenland’s mass loss and how the geographical pattern has varied on relatively shorter time scales. Here, we present a spatially and temporally resolved estimation of the ice mass change over Greenland between April of 2002 and August of 2011. Although the total mass loss trend has remained linear, actively changing areas of mass loss were concentrated on the southeastern and northwestern coasts, with ice mass in the center of Greenland steadily increasing over the decade.