Author | Title | Year | Journal/Proceedings | Reftype | DOI/URL |
---|---|---|---|---|---|

Amirbekyan, A., Michel, V. and Frederik J. Simons |
Parameterizing surface-wave tomographic models with harmonic spherical splines | 2008 | Geophys. J. Int. Vol. 174(2), pp. 617–628 |
article | DOI URL |

Abstract: We present a mathematical framework and a new methodology for the parametrization of surface wave phase-speed models, based on traveltime data. Our method is neither purely local, like block-based approaches, nor is it purely global, like those based on spherical harmonic basis functions. Rather, it combines the well-known theory and practical utility of the spherical harmonics with the spatial localization properties of spline basis functions. We derive the theoretical foundations for the application of harmonic spherical splines to surface wave tomography and summarize the results of numerous numerical tests illustrating the performance of a practical inversion scheme based upon them. Our presentation is based on the notion of reproducing-kernel Hilbert spaces, which lends itself to the parametrization of fully 3-D tomographic earth models that include body waves as well. |
|||||

BibTeX:
@article{Amirbekyan+2008b, author = {Abel Amirbekyan and Volker Michel and |
|||||

Beggan, C.D., Saarimäki, J., Whaler, K.A. and Frederik J. Simons |
Spectral and spatial decomposition of lithospheric magnetic field models using spherical Slepian functions | 2013 | Geophys. J. Int. Vol. 193(1), pp. 136–148 |
article | DOI URL |

Abstract: Global magnetic field models are typically expressed as spherical-harmonic expansion coefficients. Slepian functions are linear combinations of spherical harmonics that produce new basis functions, which vanish approximately outside chosen geographical boundaries but also remain orthogonal within the spatial region of interest. Hence, they are suitable for decomposing spherical-harmonic models into portions that have significant magnetic field strength only in selected areas. Slepian functions are spatio-spectrally concentrated, balancing spatial bias and spectral leakage. Here, we employ them as a basis to decompose the global lithospheric magnetic field model MF7 up to degree and order 72, into two distinct regions. One of the resultant fields is concentrated within the ensemble of continental domains, and the other is localized over its complement, the oceans. Our procedure neatly divides the spectral power at each harmonic degree into two parts. The field over the continents dominates the overall crustal magnetic field, and each region has a distinct power-spectral signature. The oceanic power spectrum is approximately flat, while that of the continental region shows increasing power as the spherical-harmonic degree increases. We provide a further breakdown of the field into smaller, non-overlapping continental and oceanic regions, and speculate on the source of the variability in their spectral signatures. |
|||||

BibTeX:
@article{Beggan+2013, author = {Ciarán D. Beggan and Jarno Saarimäki and Kathryn A. Whaler and |
|||||

Beveridge, A.K., Harig, C. and Frederik J. Simons |
The changing mass of glaciers on the Tibetan Plateau, 2002–2016, using time-variable gravity from the GRACE satellite mission |
2018 | J. Geod. Sci. Vol. 8, pp. 83–97 |
article | DOI URL |

Abstract: The Tibetan Plateau is the largest region of high elevation in the world. The source of water for a number of important rivers, the Himalayan region is vital to the billions of inhabitants of the Asian continent. Over the last fifty years, the climate in the region has warmed more rapidly than anywhere else at the same latitude. Causes and effects, and the geographical details of these alarming warming trends are as yet not fully known. One way of assessing the effects of climate change in this area is to measure the change in glacier volume in the region, but estimates made on the basis of different techniques have not been conclusive to date, and remain difficult to reconcile. We examine the temporal behavior of the mass flux integrated over four distinct groupings of Tibetan glaciers using satellite gravimetry from the Gravity Recovery and Cli- mate Experiment (GRACE). We use a technique known as spatio-spectral localization using spherical Slepian functions to convert global spherical harmonic expansions of the time-dependent geopotential into monthly estimates of mass changes over the Tibetan Plateau. Subsequent reductions are aimed at interpreting this mass change as due to gains or losses in ice mass. We find that (ice) mass has been decreasing on the Tibetan Plateau between 2002 and 2016 but with significant spatial variability throughout the region. Specifically, in the regions of Himalaya, Pamir, Qilian, and Tien Shan, glaciers have been losing ice mass at a rate of −11±3, −1±2, +8±2, and −6±1 Gt/yr, respectively, over the last decade. |
|||||

BibTeX:
@article{Beveridge+2018, author = {Alyson Kathryn Beveridge and Christopher Harig and |
|||||

Bevis, M., Harig, C., Khan, S.A., Brown, A., Frederik J. Simons, Willis, M., Fettweis, X., van den BroekeBroeke, M.R., Madsen, F.B., Kendrick, E.C., II, D.J.C., van Dam, T., Knudsen, P. and Nylen, T. |
Accelerating changes in ice mass within Greenland, and the ice sheet's sensitivity to atmospheric forcing | 2019 | Proc. Natl. Acad. Sc. Vol. 116(6), pp. 1934–1939 |
article | DOI URL |

Abstract: From early 2003 to mid-2013, the total mass of ice in Greenland declined at a progressively increasing rate. In mid-2013, an abrupt reversal occurred, and very little net ice loss occurred in the next 12–18 months. Gravity Recovery and Climate Experiment (GRACE) and global positioning system (GPS) observations reveal that the spatial patterns of the sustained acceleration and the abrupt deceleration in mass loss are similar. The strongest accelerations tracked the phase of the North Atlantic Oscillation (NAO). The negative phase of the NAO enhances summertime warming and insolation while reducing snowfall, especially in west Greenland, driving surface mass balance (SMB) more negative, as illustrated using the regional climate model MAR. The spatial pattern of accelerating mass changes reflects the geography of NAO-driven shifts in atmospheric forcing and the ice sheet’s sensitivity to that forcing. We infer that southwest Greenland will become a major future contributor to sea level rise. |
|||||

BibTeX:
@article{Bevis+2019, author = {Michael Bevis and Christopher Harig and Shfaqat A. Khan and Abel Brown and |
|||||

Borisov, D., Modrak, R., Rusmanugroho, H., Yuan, Y.O., Gao, F., Frederik J. Simons and Tromp, J. |
Spectral-element based 3D elastic full waveform inversion of surface waves in the presence of complex topography using an envelope-based misfit function | 2016 | SEG Tech. Prog. Expanded Abstracts, pp. 1211–1215 | inproceedings | DOI URL |

Abstract: Full-waveform inversion (FWI) is a data fitting technique used to estimate properties of the Earth from seismic data by minimizing the misfit between observed and simulated seismograms. Because of very high computational cost, this technique has so far been used either in a 2D fully elastic formulation or in a 3D acoustic formulation, when applied to active-source surveys in order to image the shallow subsurface (i.e., down to the first few kilometers). However, the Earth is three-dimensional, (visco)elastic and highly heterogeneous. Therefore, obtaining more accurate models requires solving the full 3D elastic wave equation. In this study, we use an envelope-based misfit function to construct shallow 3D models of shear wavespeed while inverting surface waves. The envelope-based misfit function has proven to be effective for inverting surface waves, which are particularly exposed to the cycle-skipping problem. To accurately model the wavefield in the presence of complex topography, we use a spectral-element wave propagation code. A synthetic example on the SEAM Phase II foothills model illustrates that inversion of surface waves at the initial stages in such a challenging environment allows us to obtain an improved shear wavespeed starting model for traditional FWI. |
|||||

BibTeX:
@inproceedings{Borisov+2016, author = {Dmitry Borisov and Ryan Modrak and Herurisa Rusmanugroho and Yanhua O. Yuan and Fuchun Gao and |
|||||

Borisov, D., Gao, F., Williamson, P., Frederik J. Simons and Tromp, J. |
Robust surface-wave full-waveform inversion | 2019 | SEG Tech. Prog. Expanded Abstracts, pp. 5005-5009 | inproceedings | DOI URL |

Abstract: Estimation of subsurface seismic properties is important in civil engineering, oil & gas exploration, and global seismology. We present a method and an application of robust surface-wave inversion in the context of 2D elastic waveform inversion of an active-source onshore dataset acquired on irregular topography. The lowest available frequency is 5 Hz. The recorded seismograms at relatively near offsets are dominated by dispersive surface waves. In exploration seismology, surface waves are generally treated as noise (“ground roll”) and removed as part of the data processing. In contrast, here, we invert surface waves to constrain the shallow parts of the shear wavespeed model. To diminish the dependence of surface waves on the initial model, we use a layer-stripping approach combined with an envelope-based misfit function. Surface waves are initially inverted using short offsets (up to 0.6 km) and over a high-frequency range (12.5–15 Hz) to constrain the shallow parts of the model. The lower-frequency components and longer offsets, which can sample deeper parts of the model, are gradually added to the process as the inversion proceeds. At the final stage, surface waves are inverted using offsets of up to 1.5 km over a band between 5–15 Hz. The final Vs model includes high-resolution features in the near surface, and shows good agreement with results from dispersion-curve analysis. The data fit is also greatly improved. |
|||||

BibTeX:
@inproceedings{Borisov+2019, author = {Dmitry Borisov and Fuchun Gao and Paul Williamson and |
|||||

Burky, A., Irving, J.C.E. and Frederik J. Simons |
Mantle transition zone receiver functions for Bermuda: Automation, quality control, and interpretation |
2021 | J. Geophys. Res. Vol. 126(3), pp. e2020JB020177 |
article | DOI URL |

Abstract: The origin of the Bermuda rise remains ambiguous, despite, or perhaps because of, the existence of sometimes incongruous seismic wave‐speed and discontinuity models in the sub‐Bermudian mantle. Hence, whether Bermuda is the surface manifestation of a mantle plume remains in question. Using the largest data set of seismic records from Bermuda to date, we estimate radial receiver functions at the Global Seismographic Network (GSN) station BBSR in multiple frequency bands, using iterative time‐domain deconvolution. Motivated by synthetic experiments using axisymmetric spectral‐element forward waveform modeling, we devise a quality metric for our receiver functions to aid in the automation and reproduction of mantle transition zone discontinuity studies. We interpret the complex signals we observe by considering the mineralogical controls on mantle transition zone discontinuity structure, and conclude that our results are likely to be indicative of a thicker than average mantle transition zone. Our result is incompatible with the canonical model of a whole mantle plume in an olivine dominated mantle; however, considerations of phase transitions in the garnet system would allow us to reconcile our observations with the possible presence of a through‐going hot thermal anomaly beneath Bermuda. |
|||||

BibTeX:
@article{Burky+2021, author = {Alexander Burky and Jessica C. E. Irving and |
|||||

Charléty, J., Voronin, S., Nolet, G., Loris, I., Frederik J. Simons, Sigloch, K. and Daubechies, I.C. |
Global seismic tomography with sparsity constraints: Comparison with smoothing and damping regularization | 2013 | J. Geophys. Res. Vol. 118(9), pp. 4887–4899 |
article | DOI URL |

Abstract: We present a realistic application of an inversion scheme for global seismic tomography that uses as prior information the sparsity of a solution, defined as having few nonzero coefficients under the action of a linear transformation. In this paper, the sparsifying transform is a wavelet transform. We use an accelerated iterative soft-thresholding algorithm for a regularization strategy, which produces sparse models in the wavelet domain. The approach and scheme we present may be of use for preserving sharp edges in a tomographic reconstruction and minimizing the number of features in the solution warranted by the data. The method is tested on a data set of time delays for finite-frequency tomography using the USArray network, the first application in global seismic tomography to real data. The approach presented should also be suitable for other imaging problems. From a comparison with a more traditional inversion using damping and smoothing constraints, we show that (1) we generally retrieve similar features, (2) fewer nonzero coefficients under a properly chosen representation (such as wavelets) are needed to explain the data at the same level of root-mean-square misfit, (3) the model is sparse or compressible in the wavelet domain, and (4) we do not need to construct a heterogeneous mesh to capture the available resolution. |
|||||

BibTeX:
@article{Charlety+2013, author = {Jean Charléty and Sergey Voronin and Guust Nolet and Ignace Loris and |
|||||

Dahlen, F.A. and Frederik J. Simons |
Spectral estimation on a sphere in geophysics and cosmology | 2008 | Geophys. J. Int. Vol. 174(3), pp. 774–807 |
article | DOI URL |

Abstract: We address the problem of estimating the spherical-harmonic power spectrum of a statistically isotropic scalar signal from noise-contaminated data on a region of the unit sphere. Three different methods of spectral estimation are considered: (i) the spherical analogue of the one-dimensional (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 trade-off between spectral resolution and variance, and is easily implemented regardless of the region size, requiring neither non-linear iteration nor large-scale matrix inversion. As a result, the method is ideally suited for most applications in geophysics, geodesy or planetary science, where the objective is to obtain a spatially localized estimate of the spectrum of a signal from noisy data within a pre-selected and typically small region. |
|||||

BibTeX:
@article{Dahlen+2008, author = {F. Anthony Dahlen and |
|||||

Freeden, W., Michel, V. and Frederik J. Simons |
Spherical-Harmonics Based Special Function Systems and Constructive Approximation Methods |
2018 | Handbook of Mathematical Geodesy, pp. 753–819 | incollection | DOI URL |

Abstract: Special function systems are reviewed that reflect particular properties of the Legendre polynomials, such as spherical harmonics, zonal kernels, and Slepian functions. The uncertainty principle is the key to their classification with respect to their localization in space and frequency/momentum. Methods of constructive approximation are outlined such as spherical harmonic and Slepian expansions, spherical spline and wavelet concepts. Regularized Functional Matching Pursuit is described as an approximation technique of combining heterogeneous systems of trial functions to a kind of a ‘best basis’. |
|||||

BibTeX:
@incollection{Freeden+2018, author = {Willi Freeden and Volker Michel and |
|||||

Galanti, E., Kaspi, Y., Frederik J. Simons, Durante, D., Parisi, M., Scott and Bolton, J. |
Determining the depth of Jupiter's Great Red Spot with Juno: A Slepian approach | 2019 | Astroph. J. Lett. Vol. 874, pp. L24 |
article | DOI URL |

Abstract: One of Jupiter’s most prominent atmospheric features, the Great Red Spot (GRS), has been observed for more than two centuries, yet little is known about its structure and dynamics below its observed cloud level. While its anticyclonic vortex appearance suggests it might be a shallow weather-layer feature, the very long time span for which it was observed implies it is likely deeply rooted, otherwise it would have been sheared apart by Jupiter’s turbulent atmosphere. Determining the GRS depth will shed light not only on the processes governing the GRS, but on the dynamics of Jupiter’s atmosphere as a whole. The Juno mission single flyby over the GRS (PJ7) discovered using microwave radiometer measurements that the GRS is at least a couple hundred kilometers deep. The next flybys over the GRS (PJ18 and PJ21), will allow high-precision gravity measurements that can be used to estimate how deep the GRS winds penetrate below the cloud level. Here we propose a novel method to determine the depth of the GRS based on the new gravity measurements and a Slepian function approach that enables an effective representation of the wind-induced spatially confined gravity signal, and an efficient determination of the GRS depth given the limited measurements. We show that with this method the gravity signal of the GRS should be detectable for wind depths deeper than 300 km, with reasonable uncertainties that depend on depth (e.g., ±100 km for a GRS depth of 1000 km). |
|||||

BibTeX:
@article{Galanti+2019, author = {Eli Galanti and Yohai Kaspi and |
|||||

Goes, S., Frederik J. Simons and Yoshizawa, K. |
Seismic constraints on temperature of the Australian uppermost mantle | 2005 | Earth Planet. Sci. Lett. Vol. 236(1–2), pp. 227–237 |
article | DOI URL |

Abstract: We derive estimates of temperature of the Australian continental mantle between 80 and 350 km depth from two published S-velocity models. Lithospheric temperatures range over about 1000 ◦C, with a large-scale correlation between temperature and tectonic age. In detail however, variations ranging from 200 to 700 ◦C occur within each tectonic province. At the current seismic resolution, strictly Proterozoic and Archean blocks do not have substantially different temperatures, nor does the Phanerozoic lithosphere east and west of the Tasman line. Temperatures close to an average (moist) MORB source mantle solidus characterize the eastern seaboard and its offshore. Differences between the temperatures derived from the two velocity models illustrate the importance of well-constrained absolute velocities and gradients for physical interpretation. The large range of lithospheric temperatures cannot be explained solely with documented variability in crustal heat production, but requires significant variations in mantle heat flow as well. |
|||||

BibTeX:
@article{Goes+2005, author = {Saskia Goes and |
|||||

Gualtieri, L., Bachmann, E., Frederik J. Simons and Tromp, J. |
The origin of secondary microseism Love waves | 2020 | Proc. Natl. Acad. Sc. Vol. 117(47), pp. 29504–-29511 |
article | DOI URL |

Abstract: The interaction of ocean surface waves produces pressure fluctuations at the seafloor capable of generating seismic waves in the solid Earth. The accepted mechanism satisfactorily explains secondary microseisms of the Rayleigh type, but it does not justify the presence of transversely polarized Love waves, nevertheless widely observed. An explanation for two-thirds of the worldwide ambient wave field has been wanting for over a century. Using numerical simulations of global-scale seismic wave propagation at unprecedented high frequency, here we explain the origin of secondary microseism Love waves. A small fraction of those is generated by boundary force-splitting at bathymetric inclines, but the majority is generated by the interaction of the seismic wave field with three-dimensional heterogeneity within the Earth. We present evidence for an ergodic model that explains observed seismic wave partitioning, a requirement for full-wave field ambient-noise tomography to account for realistic source distributions. |
|||||

BibTeX:
@article{Gualtieri+2020, author = {Lucia Gualtieri and Etienne Bachmann and |
|||||

Gualtieri, L., Bachmann, E., Frederik J. Simons and Tromp, J. |
Generation of secondary microseism Love waves: effects of bathymetry, 3D structure, and source seasonality | 2021 | Geophys. J. Int. Vol. 226(1), pp. 192-219 |
article | DOI URL |

Abstract: Secondary microseisms are ubiquitous ambient noise vibrations due to ocean activity, domi- nating worldwide seismographic records at seismic periods between 3 and 10 s. Their origin is a heterogeneous distribution of pressure fluctuations along the ocean surface. In spherically symmetric earth models, no Love surface waves are generated by such a distributed surface source. We present global-scale modelling of three-component secondary microseisms using a spectral-element method, which naturally accounts for a realistic distribution of surface sources, topography and bathymetry, and 3-D heterogeneity in Earth’s crust and mantle. Seis-mic Love waves emerge naturally once the system reaches steady state. The ergodic origin of Love waves allows us to model the horizontal components of secondary microseisms for the first time. Love waves mostly originate from the interaction of the seismic wavefield with heterogeneous Earth structure in which the mantle plays an important role despite the short periods involved. Bathymetry beneath the source region produces weak horizontal forces that are responsible for a weak and diffuse Love wavefield. The effect of bathymetric force splitting into radial and horizontal components is overall negligible when compared to the effect of 3-D heterogeneity. However, we observe small and well-focused Love-wave arrivals at seis- mographic stations in Europe due to force splitting at the steepest portion of the North Atlantic Ridge and the ocean–continent boundary. The location of the sources of Love waves is seasonal at periods shorter than about 7 s, while seasonality is lost at the longer periods. Sources of Rayleigh and Love waves from the same storm may be located very far away, indicating that energy equipartitioning might not hold in the secondary microseism period band. |
|||||

BibTeX:
@article{Gualtieri+2021, author = {Lucia Gualtieri and Etienne Bachmann and |
|||||

Han, S.-C. and Frederik J. Simons |
Spatiospectral localization of global geopotential fields from the Gravity Recovery and Climate Experiment (GRACE) reveals the coseismic gravity change owing to the 2004 Sumatra-Andaman earthquake | 2008 | J. Geophys. Res. Vol. 113, pp. B01405 |
article | DOI URL |

Abstract: Regional mass fluxes owing to transport and adjustment within the Earth system that are implicitly contained in the monthly Gravity Recovery and Climate Experiment (GRACE) global geopotential coefficients are revealed by localizing global spectra using spatiospectrally concentrated window functions. We have analyzed 45 monthly global GRACE harmonic coefficient series in order to find the coseismic signature associated with the 2004 great Sumatra-Andaman earthquake. A significant gravity change after the earthquake is found in the time series of the GRACE coefficients after localization with a single band-limited window centered near the north of the island of Sumatra. This change is undetectable from the original global coefficients or from coefficients localized elsewhere on the globe. A step function with its discontinuity at 26 December 2004 usefully models the coseismic gravity change. The localized GRACE coefficients contain the jumps (associated with the earthquake) up to degree and order 55, although not all of them within this band produce changes that are statistically significant. The gravity change calculated from the localized GRACE coefficients displays 30 mGal peak-to-peak variations that are very well correlated with an independently derived seismic model based on elastic dislocation theory. |
|||||

BibTeX:
@article{Han+2008a, author = {Shin-Chan Han and |
|||||

Harig, C., Zhong, S. and Frederik J. Simons |
Constraints on upper-mantle viscosity inferred from the flow-induced pressure gradient across a continental keel | 2010 | Geochem. Geophys. Geosys. Vol. 11(6), pp. Q06004 |
article | DOI URL |

Abstract: The thickness of continental lithosphere varies considerably from tectonically active to cratonic regions, where it can be as thick as 250–300 km. Embedded in the upper mantle like a ship, when driven to move by a velocity imposed at the surface, a continental keel is expected to induce a pressure gradient in the mantle. We hypothesize that the viscosity of the asthenosphere or the shear coupling between lower lithosphere and asthenosphere should control this pressure effect and thus the resulting dynamic topography. We perform three‐dimensional finite element calculations to examine the effects of forcing a continental keel by an imposed surface velocity, with the Australian region as a case study. When the upper mantle is strong but still weaker than the lower mantle, positive dynamic topography is created around the leading edge, and negative dynamic topography is created around the trailing edge of the keel, which is measurable by positive and negative geoid anomalies, respectively. For a weak upper mantle the effect is much reduced. We analyze geoidal and gravity anomalies in the Australian region by spatiospectral localization using Slepian functions. The method allows us to remove a best fit estimate of the geographically localized low spherical harmonic degree contributions. Regional geoid anomalies thus filtered are on the order of ±10 m across the Australian continent, with a spatial pattern similar to that predicted by the models. The comparison of modeled and observed geoid anomalies places constraints on mantle viscosity structure. Models with a two‐layer mantle cannot sufficiently constrain the ratio of viscosity between the upper and lower mantle. The addition of a third, weak, upper mantle layer, an asthenosphere, amplifies the effects of keels. Our three‐layer models, with lower mantle viscosity of 3 × 10^22 Pa s, suggest that the upper mantle (asthenosphere) is 300 times weaker than the lower mantle, while the transition zone (400–670 km depths) has a viscosity varying between 10^21 and 10^22 Pa s. |
|||||

BibTeX:
@article{Harig+2010, author = {Christopher Harig and Shijie Zhong and |
|||||

Harig, C. and Frederik J. Simons |
Mapping Greenland's mass loss in space and time | 2012 | Proc. Natl. Acad. Sc. Vol. 109(49), pp. 19934–19937 |
article | DOI URL |

Abstract: 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. |
|||||

BibTeX:
@article{Harig+2012, author = {Christopher Harig and |
|||||

Harig, C. and Frederik J. Simons |
Accelerated West Antarctic ice mass loss continues to outpace East Antarctic gains | 2015 | Earth Planet. Sci. Lett. Vol. 415, pp. 134–141 |
article | DOI URL |

Abstract: While multiple data sources have confirmed that Antarctica is losing ice at an accelerating rate, different measurement techniques estimate the details of its geographically highly variable mass balance with different levels of accuracy, spatio-temporal resolution, and coverage. Some scope remains for methodological improvements using a single data type. In this study we report our progress in increasing the accuracy and spatial resolution of time-variable gravimetry from the Gravity Recovery and Climate Experiment (GRACE). We determine the geographic pattern of ice mass change in Antarctica between January 2003 and June 2014, accounting for glacio-isostatic adjustment (GIA) using the IJ05_R2 model. Expressing the unknown signal in a sparse Slepian basis constructed by optimization to prevent leakage out of the regions of interest, we use robust signal processing and statistical estimation methods. Applying those to the latest time series of monthly GRACE solutions we map Antarctica’s mass loss in space and time as well as can be recovered from satellite gravity alone. Ignoring GIA model uncertainty, over the period 2003–2014, West Antarctica has been losing ice mass at a rate of −121±8 Gt/yr and has experienced large acceleration of ice mass losses along the Amundsen Sea coast of −18±5 Gt/yr2, doubling the mass loss rate in the past six years. The Antarctic Peninsula shows slightly accelerating ice mass loss, with larger accelerated losses in the southern half of the Peninsula. Ice mass gains due to snowfall in Dronning Maud Land have continued to add about half the amount of West Antarctica’s loss back onto the continent over the last decade. We estimate the overall mass losses from Antarctica since January 2003 at −92±10 Gt/yr. |
|||||

BibTeX:
@article{Harig+2015a, author = {Christopher Harig and |
|||||

Harig, C., Lewis, K.W., Plattner, A. and Frederik J. Simons |
A suite of software analyzes data on the sphere | 2015 | Eos Trans. AGU Vol. 96(6), pp. 18–22 |
article | DOI URL |

Abstract: The software improves data analysis over small portions of a spherical planetary surface. Among other applications, it has helped track Greenland’s ice loss over time. |
|||||

BibTeX:
@article{Harig+2015b, author = {Christopher Harig and Kevin W. Lewis and Alain Plattner and |
|||||

Harig, C. and Frederik J. Simons |
Ice mass loss in Greenland, the Gulf of Alaska, and the Canadian Archipelago: Seasonal cycles and decadal trends | 2016 | Geophys. Res. Lett. Vol. 43, pp. 3150–3159 |
article | DOI URL |

Abstract: Over the past several decades mountain glaciers and ice caps have been significant contributors to sea level rise. Here we estimate the ice mass changes in the Canadian Archipelago, the Gulf of Alaska, and Greenland since 2003 by analyzing time-varying gravimetry data from the Gravity Recovery and Climate Experiment. Prior to 2013, interannual ice mass variability in the Gulf of Alaska and in regions around Greenland remains within the average estimated over the whole data span. Beginning in summer 2013, ice mass in regions around Greenland departs positively from its long-term trend. Over Greenland this anomaly reached almost 500 Gt through the end of 2014. Overall, long-term ice mass loss from Greenland and the Canadian Archipelago continues to accelerate, while losses around the Gulf of Alaska region continue but remain steady with no significant acceleration. |
|||||

BibTeX:
@article{Harig+2016, author = {Christopher Harig and |
|||||

Kalnins, L.M., Frederik J. Simons, Kirby, J.F., Wang, D.V. and Olhede, S.C. |
On the robustness of estimates of mechanical anisotropy in the continental lithosphere: A North American case study and global reanalysis | 2015 | Earth Planet. Sci. Lett. Vol. 419, pp. 43–51 |
article | DOI URL |

Abstract: Lithospheric strength variations both influence and are influenced by many tectonic processes, including orogenesis and rifting cycles. The long, complex, and highly anisotropic histories of the continental lithosphere might lead to a natural expectation of widespread mechanical anisotropy. Anisotropy in the coherence between topography and gravity anomalies is indeed often observed, but whether it corresponds to an elastic thickness that is anisotropic remains in question. If coherence is used to estimate flexural strength of the lithosphere, the null-hypothesis of elastic isotropy can only be rejected when there is significant anisotropy in both the coherence and the elastic strengths derived from it, and if interference from anisotropy in the data themselves can be plausibly excluded. We consider coherence estimates made using multitaper and wavelet methods, from which estimates of effective elastic thickness are derived. We develop a series of statistical and geophysical tests for anisotropy, and specifically evaluate the potential for spurious results with synthetically generated data. Our primary case study, the North American continent, does not exhibit meaningful anisotropy in its mechanical strength. Similarly, a global reanalysis of continental gravity and topography using multitaper methods produces only scant evidence for lithospheric flexural anisotropy. |
|||||

BibTeX:
@article{Kalnins+2015, author = {Lara M. Kalnins and |
|||||

Kazei, V., Ovcharenko, O., Alkhalifah, T. and Frederik J. Simons |
Realistically textured random velocity models for deep learning applications | 2019 | EAGE Conf. Proc., pp. 1-5 | inproceedings | DOI URL |

Abstract: Deep learning can be used to help reconstruct low frequencies in seismic data, and to directly infer velocity models in simple cases. In order to succeed with deep learning, a good training set of velocity models is critical. We present a new way to design random models that are statistically similar to a given guiding model. Our approach is based on shuffling the coefficients of a wavelet packet decomposition (WPD) of the guiding model, allowing us to replicate realistic textures from a synthetic model. We generate realistically random models from the BP 2004 and Marmousi II models for neural network training, and utilize the trained network to extrapolate low frequencies for the SEAM model. We apply full-waveform inversion to the extrapolated data to understand the limitations of our approach. |
|||||

BibTeX:
@inproceedings{Kazei+2019, author = {Vladimir Kazei and Oleg Ovcharenko and Tariq Alkhalifah and |
|||||

Kopp, R.E., Frederik J. Simons, Mitrovica, J.X., Maloof, A.C. and Oppenheimer, M. |
Probabilistic assessment of sea level during the last interglacial stage | 2009 | Nature Vol. 462, pp. 863–867 |
article | DOI URL |

Abstract: 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 (greater than or equal to 10 m), the millennial average rate of global sea level rise is very likely to have exceeded 5.6 m/kyr but is unlikely to have exceeded 9.2 m/kyr. 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. |
|||||

BibTeX:
@article{Kopp+2009, author = {R. E. Kopp and |
|||||

Kopp, R.E., Frederik J. Simons, Mitrovica, J.X., Maloof, A.C. and Oppenheimer, M. |
A probabilistic assessment of sea level variations within the last interglacial stage | 2013 | Geophys. J. Int. Vol. 193(2), pp. 711–716 |
article | DOI URL |

Abstract: The last interglacial stage (LIG; ca. 130–115 ka) provides a relatively recent example of a world with both poles characterized by greater-than-Holocene temperatures similar to those expected later in this century under a range of greenhouse gas emission scenarios. Previous analyses inferred that LIG mean global sea level (GSL) peaked 6–9 m higher than today. Here, we extend our earlier work to perform a probabilistic assessment of sea level variability within the LIG highstand. Using the terminology for probability employed in the Intergovernmental Panel on Climate Change assessment reports, we find it extremely likely (95 per cent probability) that the palaeo-sea level record allows resolution of at least two intra-LIG sea level peaks and likely (67 per cent probability) that the magnitude of low-to-high swings exceeded 4 m. Moreover, it is likely that there was a period during the LIG in which GSL rose at a 1000-yr average rate exceeding 3 m/kyr, but unlikely (33 per cent probability) that the rate exceeded 7 m/kyr and extremely unlikely (5 per cent probability) that it exceeded 11 m/kyr. These rate estimates can provide insight into rates of Greenland and/or Antarctic melt under climate conditions partially analogous to those expected in the 21st century. |
|||||

BibTeX:
@article{Kopp+2013, author = {R. E. Kopp and |
|||||

Lewis, K.W. and Frederik J. Simons |
Spatial variability of the Martian crustal magnetic field | 2011 | 42nd Lunar Planetary Science Conference, pp. 2621 | inproceedings | URL |

Abstract: The crust of Mars retains heterogenous remanent magnetism. Magnetic power spectra can provide constraints on the depths and strengths of magnetic sources. We use a spatiospectral windowing approach to map local variability across the planet. |
|||||

BibTeX:
@inproceedings{Lewis+2011, author = {Kevin W. Lewis and |
|||||

Lewis, K.W. and Frederik J. Simons |
Local spectral variability and the origin of the Martian crustal magnetic field | 2012 | Geophys. Res. Lett. Vol. 39, pp. L18201 |
article | DOI URL |

Abstract: The crustal remanent magnetic field of Mars remains enigmatic in many respects. Its heterogeneous surface distribution points to a complex history of formation a modification, and has been resistant to attempts at identifying magnetic paleopoles and constraining the geologic origin of crustal sources. We use a multitaper technique to quantify the spatial diversity of the field via the localized magnetic power spectrum, allowing us to isolate more weakly magnetized regions and characterize them spectrally for the first time. We find clear geographical differences in spectral properties and parameterize them in terms of source strengths and equivalent-layer decorrelation depths. Depths to the magnetic layer in our model correlate with independent crustal-thickness estimates. The correspondence indicates that a significant fraction of the martian crustal column may contribute to the observed field, as would be consistent with an intrusive magmatic origin. We identify several anomalous regions, and propose geophysical mechanisms for their spectral signatures. |
|||||

BibTeX:
@article{Lewis+2012, author = {Kevin W. Lewis and |
|||||

Lewis, K.W., Frederik J. Simons and Eggers, G.L. |
Maximum-likelihood estimation of lithospheric thickness on Venus | 2013 | 44th Lunar Planetary Science Conference, pp. 2612 | inproceedings | URL |

Abstract: We utilize a new, maximum likelihood-based technique to estimate elastic thickness and loading characteristics of the venusian lithosphere. |
|||||

BibTeX:
@inproceedings{Lewis+2013, author = {Kevin W.~Lewis and |
|||||

Lewis, K.W., Frederik J. Simons, Olhede, S.C. and Eggers, G.L. |
Maximum-likelihood analysis of planetary roughness | 2017 | 48th Lunar Planetary Science Conference, pp. 2608 | inproceedings | URL |

Abstract: We have developed a new computational method for the characterization of the spatial statistics of planetary data fields, including topographic roughness. |
|||||

BibTeX:
@inproceedings{Lewis+2017, author = {Kevin W. Lewis and |
|||||

Liu, Z., Hoffmann, J., Frederik J. Simons and Tromp, J. |
Elastic full waveform inversion of VSP data from a complex anticline in northern Iraq | 2021 | SEG Tech. Prog. Expanded Abstracts, pp. XXX | inproceedings | |

Abstract: We demonstrate an application of isotropic elastic Full-Waveform Inversion (FWI) to a field data set of Vertical Seismic Profiles (VSP) from a structurally complex narrow anticline in Northern Iraq. A practical elastic FWI workflow is developed to invert offset VSP seismic data. Synthetic tests indicate that this approach can give accurate wave speed updates for the target anticlinal structure. We use inverted wave speed models to perform elastic Reverse-Time Migration (RTM) of the data. The RTM results show that the shear wave speed (VS) image has a higher resolution than the compressional speed (VP) image for the target structure, owing to the presence of interpretable P-to-S converted waves. The results of the field test show that elastic FWI produces elastic models from which we obtain clear VS images of the thrust fault of interest. |
|||||

BibTeX:
@inproceedings{Liu+2021, author = {Zhaolun Liu and Jürgen Hoffmann and |
|||||

Maloof, A.C., Rose, C.V., Calmet, C.C., Beach, R., Samuels, B.M., Erwin, D.H., Poirier, G.R., Yao, N. and Frederik J. Simons |
Possible animal body-fossils from pre-Marinoan limestones, South Australia | 2010 | Nature Geosci. Vol. 3(9), pp. 653–659 |
article | DOI URL |

Abstract: The Neoproterozoic era was punctuated by the Sturtian (about 710 million years ago) and Marinoan (about 635 million years ago) intervals of glaciation. In South Australia, the rocks left behind by the glaciations are separated by a succession of limestones and shales, which were deposited at tropical latitudes. Here we describe millimetre-to centimetre-scale fossils from the Trezona Formation, which pre-dates the Marinoan glaciation. These weakly calcified fossils occur as anvil, wishbone, ring and perforated slab shapes and are contained within stromatolitic limestones. The Trezona Formation fossils pre-date the oldest known calcified fossils of this size by 90 million years, and cannot be separated from the surrounding calcite matrix or imaged by traditional X-ray-based tomographic scanning methods. Instead, we have traced cross-sections of individual fossils by serially grinding and scanning each sample at a resolution of 50.8 μm. From these images we constructed three-dimensional digital models of the fossils. Our reconstructions show a population of ellipsoidal organisms without symmetry and with a network of interior canals that lead to circular apertures on the fossil surface. We suggest that several characteristics of these reef-dwelling fossils are best explained if the fossils are identified as sponge-grade metazoans. |
|||||

BibTeX:
@article{Maloof+2010, author = {Adam C. Maloof and Catherine V. Rose and Claire C. Calmet and Robert Beach and Bradley M. Samuels and Douglas H. Erwin and Gerald R. Poirier and Nan Yao and |
|||||

McGuire, J.J., Frederik J. Simons and Collins, J.A. |
Analysis of seafloor seismograms of the 2003 Tokachi-Oki earthquake sequence for earthquake early warning | 2008 | Geophys. Res. Lett. Vol. 35, pp. L14310 |
article | DOI URL |

Abstract: Earthquake Early Warning (EEW) algorithms estimate the magnitude of an underway rupture from the first few seconds of the P-wave to allow hazard assessment and mitigation before the S-wave arrival. Many large subduction-zone earthquakes initiate 50–150 km offshore, potentially allowing seafloor instruments sufficient time to identify large ruptures before the S-waves reach land. We tested an EEW algorithm using accelerograms recorded offshore Hokkaido in the region of the 2003 Mw 8.1 Tokachi-Oki earthquake and its aftershocks. A wavelet transform of the first 4 s of the P-wave concentrates information about earthquake magnitude from both waveform amplitude and frequency content. We find that wavelets with support of a few seconds provide discriminants for EEW that are both accurate enough to be useful and superior to peak acceleration or peak velocity. Additionally, we observe a scaling of wavelet coefficient magnitude above Mw 6.0 indicating that, at least for the mainshock (Mw 8.1) and largest aftershock (Mw 7.1), the final size of a rupture could have been estimated from the initial portion of the seismogram. |
|||||

BibTeX:
@article{McGuire+2008, author = {Jeffrey J. McGuire and |
|||||

Michel, V. and Frederik J. Simons |
A general approach to regularizing inverse problems with regional data using Slepian wavelets |
2017 | Inv. Probl. Vol. 33, pp. 125016 |
article | DOI URL |

Abstract: Slepian functions are orthogonal function systems that live on subdomains (for example, geographical regions on the Earth’s surface, or bandlimited portions of the entire spectrum). They have been firmly established as a useful tool for the synthesis and analysis of localized (concentrated or confined) signals, and for the modeling and inversion of noise-contaminated data that are only regionally available or only of regional interest. In this paper, we consider a general abstract setup for inverse problems represented by a linear and compact operator between Hilbert spaces with a known singular-value decomposition (svd). In practice, such an svd is often only given for the case of a global expansion of the data (e.g. on the whole sphere) but not for regional data distributions. We show that, in either case, Slepian functions (associated to an arbitrarily prescribed region and the given compact operator) can be determined and applied to construct a regularization for the ill-posed regional inverse problem. Moreover, we describe an algorithm for constructing the Slepian basis via an algebraic eigenvalue problem. The obtained Slepian functions can be used to derive an svd for the combination of the regionalizing projection and the compact operator. As a result, standard regularization techniques relying on a known svd become applicable also to those inverse problems where the data are regionally given only. In particular, wavelet-based multiscale techniques can be used. An example for the latter case is elaborated theoretically and tested on two synthetic numerical examples. |
|||||

BibTeX:
@article{Michel+2017, author = {Volker Michel and |
|||||

Montési, L.G.J., di Toro, G., Frederik J. Simons, Akber-Knutson, S., Becker, T.W., Billen, M., Deschamps, A. and Kellogg, J.B. |
Young scientists focus on the dynamics of the lithosphere | 2006 | Eos Trans. AGU Vol. 87(44), pp. 482 |
article | DOI URL |

Abstract: Young researchers face a wide choice of scientific approaches and directions that may shape their careers. The Earth sciences, in particular, offer a broad range of topics to study and techniques to use that exceeds what any current scientist was exposed to at the schools they attended. At early stages of their careers, researchers need to gain confidence in their expertise, and they often can benefit by expanding their scientific horizons and forging new collaborations. |
|||||

BibTeX:
@article{Montesi+2006, author = {Laurent G. J. Montési and Giulio di Toro and |
|||||

Nolet, G., Hello, Y., van der Lee, S., Bonnieux, S., Ruiz, M.C., Pazmino, N.A., Deschamps, A., Regnier, M.M., Font, Y., Chen, Y.J. and Frederik J. Simons |
Imaging the Galápagos mantle plume with an unconventional application of floating seismometers | 2019 | Sci. Rep. Vol. 9, pp. 1326 |
article | DOI URL |

Abstract: We launched an array of nine freely floating submarine seismometers near the Galápagos islands, which remained operational for about two years. P and PKP waves from regional and teleseismic earthquakes were observed for a range of magnitudes. The signal-to-noise ratio is strongly influenced by the weather conditions and this determines the lowest magnitudes that can be observed. Waves from deep earthquakes are easier to pick, but the S/N ratio can be enhanced through filtering and the data cover earthquakes from all depths. We measured 580 arrival times for different raypaths. We show that even such a limited number of data gives a significant increase in resolution for the oceanic upper mantle. This is the first time an array of floating seismometers is used in seismic tomography to improve the resolution significantly where otherwise no seismic information is available. We show that the Galápagos Archipelago is underlain by a deep (about 1900 km) 200–300 km wide plume of high temperature, with a heat flux very much larger than predicted from its swell bathymetry. The decrease of the plume temperature anomaly towards the surface indicates that the Earth’s mantle has a subadiabatic temperature gradient. |
|||||

BibTeX:
@article{Nolet+2019, author = {Guust Nolet and Yann Hello and Suzan van der Lee and Sébastien Bonnieux and Mario C. Ruiz and Nelson |
|||||

Plattner, A., Frederik J. Simons and Wei, L. |
Analysis of real vector fields on the sphere using Slepian functions | 2012 | 2012 IEEE Statistical Signal Processing Workshop (SSP'12) | inproceedings | DOI URL |

Abstract: We pose and solve the analogue of Slepian's time-frequency concentration problem for vector fields on the surface of the unit sphere, to determine an orthogonal family of strictly bandlimited vector fields that are optimally concentrated within a closed region of the sphere or, alternatively, of strictly spacelimited functions that are optimally concentrated in the vector 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. |
|||||

BibTeX:
@inproceedings{Plattner+2012b, author = {Alain Plattner and |
|||||

Plattner, A. and Frederik J. Simons |
A spatiospectral localization approach for analyzing and representing vector-valued functions on spherical surfaces | 2013 | Wavelets and Sparsity XV, pp. 88580N | inproceedings | DOI URL |

Abstract: We review the construction of three different Slepian bases on the sphere, and illustrate their theoretical behavior and practical use for solving ill-posed satellite inverse problems. The first basis is scalar, the second vectorial, and the third suitable for the vector representation of the harmonic potential fields on which we focus our analysis. When data are noisy and incompletely observed over contiguous domains covering parts of the sphere at satellite altitude, expanding the unknown solution in terms of a Slepian basis and seeking truncated expansions to achieve least-squares data fit has advantages over conventional approaches that include the ease with which the solutions can be computed, and a clear statistical understanding of the competing effects of solution bias and variance in modulating the mean squared error, as we illustrate with several new examples. |
|||||

BibTeX:
@inproceedings{Plattner+2013, author = {Alain Plattner and |
|||||

Plattner, A. and Frederik J. Simons |
Spatiospectral concentration of vector fields on a sphere | 2014 | Appl. Comput. Harmon. Anal. Vol. 36, pp. 1–22 |
article | DOI URL |

Abstract: We construct spherical vector bases that are bandlimited and spatially concentrated, or alternatively, spacelimited and spectrally concentrated, suitable for the analysis and representation of real-valued vector fields on the surface of the unit sphere, as arises in the natural and biomedical sciences, and engineering. Building on the original approach of Slepian, Landau, and Pollak we concentrate the energy of our function bases into arbitrarily shaped regions of interest on the sphere and within a certain bandlimit in the vector spherical-harmonic domain. As with the concentration problem for scalar functions on the sphere, which has been treated in detail elsewhere, the vector basis can be constructed by solving a finite-dimensional algebraic eigenvalue problem. The eigenvalue problem decouples into separate problems for the radial, and tangential components. For regions with advanced symmetry such as latitudinal polar caps, the spectral concentration kernel matrix is very easily calculated and block-diagonal, which lends itself to efficient diagonalization. The number of spatiospectrally well-concentrated vector fields is well estimated by a Shannon number that only depends on the area of the target region and the maximal spherical harmonic degree or bandwidth. The Slepian spherical vector basis is doubly orthogonal, both over the entire sphere and over the geographic target regions. Like its scalar counterparts it should be a powerful tool in the inversion, approximation and extension of bandlimited fields on the sphere: vector fields such as gravity and magnetism in the earth and planetary sciences, or electromagnetic fields in optics, antenna theory and medical imaging. |
|||||

BibTeX:
@article{Plattner+2014a, author = {Alain Plattner and |
|||||

Plattner, A. and Frederik J. Simons |
High-resolution local magnetic field models for the Martian South Pole from Mars Global Surveyor data | 2015 | J. Geophys. Res. Vol. 120, pp. 1543–1566 |
article | DOI URL |

Abstract: We present two high-resolution local models for the crustal magnetic field of the Martian south polar region. Models SP130 and SP130M were derived from three-component measurements made by Mars Global Surveyor at nighttime and at low altitude (<200 km). The availability area for these data covers the annulus between latitudes −76◦ and −87◦ and contains a strongly magnetized region (southern parts of Terra Sirenum) adjacent to weakly magnetized terrains (such as Prometheus Planum). Our localized field inversions take into account the region of data availability, a finite spectral bandlimit (spherical harmonic degree L = 130), and the varying satellite altitude at each observation point. We downward continue the local field solutions to a sphere of Martian polar radius 3376 km. While weakly magnetized areas in model SP130 contain inversion artifacts caused by strongly magnetized crust nearby, these artifacts are largely avoided in model SP130M, a mosaic of inversion results obtained by independently solving for the fields over individual subregions. Robust features of both models are magnetic stripes of alternating polarity in southern Terra Sirenum that end abruptly at the rim of Prometheus Planum, an impact crater with a weak or undetectable magnetic field. From a prominent and isolated dipole-like magnetic feature close to Australe Montes, we estimate a paleopole with a best fit location at longitude 207◦ and latitude 48◦. From the abruptly ending magnetic field stripes, we estimate average magnetization values of up to 15 A/m. |
|||||

BibTeX:
@article{Plattner+2015a, author = {Alain Plattner and |
|||||

Plattner, A. and Frederik J. Simons |
Mars' heterogeneous South Polar magnetic field revealed using altitude vector Slepian functions | 2015 | 46th Lunar Planetary Science Conference, pp. 1794 | inproceedings | URL |

Abstract: We present a high-resolution regional model of the martian south polar magnetic field by taking advantage of locally available high-quality data. |
|||||

BibTeX:
@inproceedings{Plattner+2015b, author = {Alain Plattner and |
|||||

Plattner, A. and Frederik J. Simons |
Potential-field estimation using scalar and vector Slepian functions at satellite altitude | 2015 | Handbook of Geomathematics, pp. 2003–2055 | incollection | DOI URL |

Abstract: In the last few decades a series of increasingly sophisticated satellite missions has brought us gravity and magnetometry data of ever improving quality. To make optimal use of this rich source of information on the structure of Earth and other celestial bodies, our computational algorithms should be well matched to the specific properties of the data. In particular, inversion methods require specialized adaptation if the data are only locally available, their quality varies spatially, or if we are interested in model recovery only for a specific spatial region. Here, we present two approaches to estimate potential fields on a spherical Earth, from gradient data collected at satellite altitude. Our context is that of the estimation of the gravitational or magnetic potential from vector-valued measurements. Both of our approaches utilize spherical Slepian functions to produce an approximation of local data at satellite altitude, which is subsequently transformed to the Earth's spherical reference surface. The first approach is designed for radial-component data only, and uses scalar Slepian functions. The second approach uses all three components of the gradient data and incorporates a new type of vectorial spherical Slepian functions which we introduce in this chapter. |
|||||

BibTeX:
@incollection{Plattner+2015c, author = {Alain Plattner and |
|||||

Plattner, A. and Frederik J. Simons |
Internal and external potential field estimation from regional vector data at varying satellite altitude | 2017 | Geophys. J. Int. Vol. 211, pp. 207–238 |
article | DOI URL |

Abstract: When modeling satellite data to recover a global planetary magnetic or gravitational potential field, the method of choice remains their analysis in terms of spherical harmonics. When only regional data are available, or when data quality varies strongly with geographic location, the inversion problem becomes severely ill-posed. In those cases, adopting explicitly local methods is to be preferred over adapting global ones (e.g., by regularization). Here, we develop the theory behind a procedure to invert for planetary potential fields from vector observations collected within a spatially bounded region at varying satellite altitude. Our method relies on the construction of spatiospectrally localized bases of functions that mitigate the noise amplification caused by downward continuation (from the satellite altitude to the source) while balancing the conflicting demands for spatial concentration and spectral limitation. The `altitude-cognizant’ gradient vector Slepian functions (AC-GVSF) enjoy a noise tolerance under downward continuation that is much improved relative to the `classical’ gradient vector Slepian functions (CL-GVSF), which do not factor satellite altitude into their construction. Furthermore, venturing beyond the realm of their first application, published in a preceding paper, in the present article we extend the theory to being able to handle both internal and external potential-field estimation. Solving simultaneously for internal and external fields under the limitation of regional data availability reduces internal-field artifacts introduced by downward-continuing unmodeled external fields, as we show with numerical examples. We explain our solution strategies on the basis of analytic expressions for the behavior of the estimation bias and variance of models for which signal and noise are uncorrelated, (essentially) space- and bandlimited, and spectrally (almost) white. The AC-GVSF are optimal linear combinations of vector spherical harmonics. Their construction is not altogether very computationally demanding when the concentration domains (the regions of spatial concentration) have circular symmetry, e.g., on spherical caps or rings — even when the spherical-harmonic bandwidth is large. Data inversion proceeds by solving for the expansion coefficients of truncated function sequences, by least-squares analysis in a reduced-dimensional space. Hence, our method brings high-resolution regional potential-field modeling from incomplete and noisy vector-valued satellite data within reach of contemporary desktop machines. |
|||||

BibTeX:
@article{Plattner+2017a, author = {Alain Plattner and |
|||||

Plattner, A., Golabek, G.J. and Frederik J. Simons |
A spectral view of the Terra Sirenum / Cimmeria crustal magnetic field |
2017 | 48th Lunar Planetary Science Conference, pp. 1627 | inproceedings | URL |

Abstract: Spatial distribution of regional power spectra provides a new perspective of the strong crustal magnetic field of the Terra Sirenum/Cimmeria region on Mars. |
|||||

BibTeX:
@inproceedings{Plattner+2017b, author = {Alain Plattner and Gregor J. Golabek and |
|||||

Reuber, G.S. and Frederik J. Simons |
Multi-physics adjoint modeling of Earth structure: combining gravimetric, seismic, and geodynamic inversions | 2018 | Intern. J. Geomath. Vol. 11(30), pp. 1-38 |
article | DOI URL |

Abstract: We discuss the resolving power of three geophysical imaging and inversion techniques, and their combination, for the reconstruction of material parameters in the Earth’s subsurface. The governing equations are those of Newton and Poisson for gravitational problems, the acoustic wave equation under Hookean elasticity for seismology, and the geodynamics equations of Stokes for incompressible steady-state flow in the mantle. The observables are the gravitational potential, the seismic displacement, and the surface velocity, all measured at the surface. The inversion parameters of interest are the mass density, the acoustic wave speed, and the viscosity. These systems of partial differential equations and their adjoints were implemented in a single Python code using the finite-element library FeNICS. To investigate the shape of the cost functions, we present a grid search in the parameter space for three end-member geological settings: a falling block, a subduction zone, and a mantle plume. The performance of a gradient-based inversion for each single observable separately, and in combination, is presented. We furthermore investigate the performance of a shape-optimizing inverse method, when the material is known, and an inversion that inverts for the material parameters of an anomaly with known shape. |
|||||

BibTeX:
@article{Reuber+2020, author = {Georg S. Reuber and |
|||||

Sambridge, M., Beghein, C., Frederik J. Simons and Snieder, R. |
How do we understand and visualize uncertainty? | 2006 | The Leading Edge Vol. 25(5), pp. 542–546 |
article | DOI URL |

Abstract: Geophysicists are often concerned with reconstructing subsurface properties using observations collected at or near the surface. |
|||||

BibTeX:
@article{Sambridge+2006, author = {Sambridge, M. and Beghein, C. and |
|||||

Simon, J.D., Frederik J. Simons and Nolet, G. |
Multiscale estimation of event arrival times and their uncertainties in hydroacoustic records from autonomous oceanic floats | 2020 | Bull. Seism. Soc. Am. Vol. 100(3), pp. 970-997 |
article | DOI URL |

Abstract: We describe an algorithm to pick event onsets in noisy records, characterize their error distributions, and derive confidence intervals on their timing. Our method is based on an Akaike information criterion that identifies the partition of a time series into a noise and a signal segment that maximizes the signal-to-noise ratio. The distinctive feature of our approach lies in the timing uncertainty analysis, and in its application in the time domain and in the wavelet timescale domain. Our novel data are records collected by freely floating Mobile Earthquake Recording in Marine Areas by Independent Divers (MERMAID) instruments, midcolumn hydrophones that report triggered segments of ocean-acoustic time series. |
|||||

BibTeX:
@article{Simon+2020, author = {Joel D. Simon and |
|||||

Frederik J. Simons, Zuber, M.T. and Korenaga, J. |
Isostatic response of the Australian lithosphere: Estimation of effective elastic thickness and anisotropy using multitaper spectral analysis | 2000 | J. Geophys. Res. Vol. 105(B8), pp. 19163–19184 |
article | DOI URL |

Abstract: Gravity and topography provide important insights regarding the degree and mechanisms of isostatic compensation. The azimuthally isotropic coherence function between 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 coherence. 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. |
|||||

BibTeX:
@article{Simons+2000, author = { |
|||||

Frederik J. Simons and van der Hilst, R.D. |
Age-dependent seismic thickness and mechanical strength of the Australian lithosphere | 2002 | Geophys. Res. Lett. Vol. 29(11), pp. 1529 |
article | DOI URL |

Abstract: We present constraints on the regional variations of the seismic and mechanical thickness of the Australian lithosphere. We infer the seismic thickness from a waveform tomographic model of S-wave speed, and as a proxy for the elastic thickness we use the wavelength at which the coherence of surface topography and Bouguer gravity drops below half of its long-wavelength maximum. Our results show that on scales <1000 km the relationship between the age of the crust and the thickness of the lithosphere is more complicated than longer-wavelength or global averages suggest. Recent geochemical and geodynamical evidence for small-scale secular variations of the composition and stability of continental cratons further illustrates the complexity of the age dependence of seismo-mechanical lithospheric properties on regional scales. |
|||||

BibTeX:
@article{Simons+2002a, author = { |
|||||

Frederik J. Simons, van der Hilst, R.D., Montagner, J.-P. and Zielhuis, A. |
Multimode Rayleigh wave inversion for heterogeneity and azimuthal anisotropy of the Australian upper mantle | 2002 | Geophys. J. Int. Vol. 151(3), pp. 738–754 |
article | DOI URL |

Abstract: We present an azimuthally anisotropic 3-D shear-wave speed model of the Australian upper mantle obtained from the dispersion of fundamental and higher modes of Rayleigh waves. We compare two tomographic techniques to map path-average earth models into a 3-D model for heterogeneity and azimuthal anisotropy. Method I uses a rectangular surface cell parametrization and depth basis functions that represent independently constrained estimates of radial earth structure. It performs an iterative inversion with norm damping and gradient regularization. Method II uses a direct inversion of individual depth layers constrained by Bayesian assumptions about the model covariance. We recall that Bayesian inversions and discrete regularization approaches are theoretically equivalent, and with a synthetic example we show that they can give similar results. The model we present here uses the discrete regularized inversion of independent path constraints of Method I, on an equal-area grid. With the exception of westernmost Australia, we can retrieve structure on length scales of about 250 km laterally and 50 km in the radial direction, to within 0.8 per cent for the velocity, 20 per cent for the anisotropic magnitude and 20◦ for its direction. On length scales of 1000 km and longer, down to about 200 km, there is a good correlation between velocity heterogeneity and geologic age. At shorter length scales and at depths below 200 km, however, this relationship breaks down. The observed magnitude and direction of maximum anisotropy do not, in general, appear to be correlated to surface geology. The pattern of anisotropy appears to be rather complex in the upper 150 km, whereas a smoother pattern of fast axes is obtained at larger depth. If some of the deeper directions of anisotropy are aligned with the approximately N–S direction of absolute plate motion, this correspondence is not everywhere obvious, despite the fast (7 cm per yr) northward motion of the Australian plate. More research is needed to interpret our observations in terms of continental deformation. Predictions of SKS splitting times and directions, an integrated measure of anisotropy, are poorly matched by observations of shear-wave birefringence. |
|||||

BibTeX:
@article{Simons+2002b, author = { |
|||||

Frederik J. Simons, van der Hilst, R.D. and Zuber, M.T. |
Spatiospectral localization of isostatic coherence anisotropy in Australia and its relation to seismic anisotropy: Implications for lithospheric deformation | 2003 | J. Geophys. Res. Vol. 108(B5), pp. 2250 |
article | DOI URL |

Abstract: 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 |
|||||

BibTeX:
@article{Simons+2003a, author = { |
|||||

Frederik J. Simons and van der Hilst, R.D. |
Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere | 2003 | Earth Planet. Sci. Lett. Vol. 211(3–4), pp. 271–286 |
article | DOI URL |

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. |
|||||

BibTeX:
@article{Simons+2003b, author = { |
|||||

Frederik J. Simons, Becker, T.W., Kellogg, J.B., Billen, M., Hardebeck, J., Lee, C.-T.A., Montési, L.G.J., Panero, W. and Zhong, S. |
Young solid Earth researchers of the world unite! | 2004 | Eos Trans. AGU Vol. 85(16), pp. 160–161 |
article | DOI URL |

Abstract: In early January 2004, one of us attended a workshop on ``science priorities and educational opportunities that can be addressed using ocean observatories.'' The attendees constituted a broad group—men and women, scientists, engineers, educators, representatives from the private and public sector—but lacked diversity in at least one important aspect: age. |
|||||

BibTeX:
@article{Simons+2004, author = { |
|||||

Frederik J. Simons, Becker, T.W., Kellogg, J.B., Billen, M., Lee, C.-T.A., Montési, L.G.J., Panero, W. and Zhong, S. |
MYRES: A program to unite young solid earth researchers | 2005 | Eos Trans. AGU Vol. 86(5), pp. 48–49 |
article | DOI URL |

Abstract: The first Meeting of Young Researchers in the Earth Sciences (MYRES-I),held in August of 2004, focused on “Heat, helium, hotspots, and whole mantle convection.”Biennial meetings, with MYRES-I as the first, are one of the ways the MYRES initiative is building an “international, interdisciplinary, open and unbiased community of colleagues who interact regularly to informally exchange ideas, data, and tools, and formulate new collaborative research projects” (see Young Solid Earth Researchers of the World Unite! published in Eos, 85(16), 160, 2004). This article reports on our first workshop, discusses what is happening in the community, and calls for proposals to keep MYRES funded. |
|||||

BibTeX:
@article{Simons+2005, author = { |
|||||

Frederik J. Simons, Dahlen, F.A. and Wieczorek, M.A. |
Spatiospectral concentration on a sphere | 2006 | SIAM Rev. Vol. 48(3), pp. 504–536 |
article | DOI URL |

Abstract: We pose and solve the analogue of Slepian’s 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. |
|||||

BibTeX:
@article{Simons+2006a, author = { |
|||||

Frederik J. Simons and Dahlen, F.A. |
Spherical Slepian functions and the polar gap in geodesy | 2006 | Geophys. J. Int. Vol. 166(3), pp. 1039–1061 |
article | DOI URL |

Abstract: 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 Slepian's 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. |
|||||

BibTeX:
@article{Simons+2006b, author = { |
|||||

Frederik J. Simons, Nolet, G., Babcock, J.M., Davis, R.E. and Orcutt, J.A. |
A future for drifting seismic networks | 2006 | Eos Trans. AGU Vol. 87(31), pp. 305 \& 307 |
article | DOI URL |

Abstract: Earth models in which seismic wave speeds vary only with depth are sufficiently well constrained to accurately locate earthquakes and calculate the paths followed by seismic rays. The differences between observations and theoretical predictions of seismograms in such one-dimensional Earth models can be used to reconstruct the three-dimensional (3-D) wave speed distribution in the regions sampled by the seismic waves by a procedure known as seismic tomography, a technique akin to medical CAT scanning. Unequal geographical data coverage continues to fundamentally limit the quality of tomographic reconstructions of seismic wave speeds in the interior of the Earth. As a possible solution to gaining equal geographic data coverage, a prototype of a mobile receiver that will serve as a floating seismometer has recently been developed. This type of instrument could provide an easy, cost-effective way to collect seismic data in the ocean. |
|||||

BibTeX:
@article{Simons+2006c, author = { |
|||||

Frederik J. Simons, Dando, B.D.E. and Allen, R.M. |
Automatic detection and rapid determination of earthquake magnitude by wavelet multiscale analysis of the primary arrival | 2006 | Earth Planet. Sci. Lett. Vol. 250(1–2), pp. 214–223 |
article | DOI URL |

Abstract: Earthquake early warning systems must save lives. It is of great importance that networked systems of seismometers be equipped with reliable tools to make rapid determinations of earthquake magnitude in the few to tens of seconds before the damaging ground motion occurs. A new fully automated algorithm based on the discrete wavelet transform detects as well as analyzes the incoming first arrival with great accuracy and precision, estimating the final magnitude to within a single unit from the first few seconds of the P wave. |
|||||

BibTeX:
@article{Simons+2006d, author = { |
|||||

Frederik J. Simons and Dahlen, F.A. |
A spatiospectral localization approach to estimating potential fields on the surface of a sphere from noisy, incomplete data taken at satellite altitudes | 2007 | Wavelets XII, pp. 670117 | inproceedings | DOI URL |

Abstract: Satellites mapping the spatial variations of the gravitational or magnetic fields of the Earth or other planets ideally fly on polar orbits, uniformly covering the entire globe. Thus, potential fields on the sphere are usually expressed in spherical harmonics, basis functions with global support. For various reasons, however, inclined orbits are favorable. These leave a ``polar gap'': an antipodal pair of axisymmetric polar caps without any data coverage, typically smaller than 10◦ in diameter for terrestrial gravitational problems, but 20◦ or more in some planetary magnetic configurations. The estimation of spherical harmonic field coefficients from an incompletely sampled sphere is prone to error, since the spherical harmonics are not orthogonal over the partial domain of the cut sphere. Although approaches based on wavelets have gained in popularity in the last decade, we present a method for localized spherical analysis that is firmly rooted in spherical harmonics. We construct a basis of bandlimited spherical functions that have the majority of their energy concentrated in a subdomain of the unit sphere by solving Slepian's (1960) concentration problem in spherical geometry, and use them for the geodetic problem at hand. Most of this work has been published by us elsewhere. Here, we highlight the connection of the ``spherical Slepian basis'' to wavelets by showing their asymptotic self-similarity, and focus on the computational considerations of calculating concentrated basis functions on irregularly shaped domains. |
|||||

BibTeX:
@inproceedings{Simons+2007, author = { |
|||||

Frederik J. Simons, Nolet, G., Georgief, P., Babcock, J.M., Regier, L.A. and Davis, R.E. |
On the potential of recording earthquakes for global seismic tomography by low-cost autonomous instruments in the oceans | 2009 | J. Geophys. Res. Vol. 114, pp. B05307 |
article | DOI URL |

Abstract: We describe the development and testing of an autonomous device designed to revolutionize Earth structure determination via global seismic tomography by detecting earthquakes at teleseismic distances in the oceans. One prototype MERMAID, short for Mobile Earthquake Recording in Marine Areas by Independent Divers, was constructed and tested at sea. The instrument combines two readily available, relatively low-cost but state-of-the-art components: a Sounding Oceanographic Lagrangian Observer, or SOLO float, and an off-the-shelf hydrophone, with custom-built data logging hardware. We report on the development of efficient wavelet-based algorithms for the detection and discrimination of seismic events and analyze three time series of acoustic pressure collected at a depth of 700 m in pilot experiments conducted offshore San Diego, CA. In these tests, over 120 hours of data were gathered, and five earthquakes, of which one was teleseismic, were recorded and identified. Quantitative estimates based on these results suggest that instruments of the MERMAID type may collect up to a hundred tomographically useful teleseismic events per year. The final design will also incorporate a Global Positioning System receiver, onboard signal processing software optimized for low-power chips, and high-throughput satellite communication equipment for telemetered data transfer. With these improvements, we hope to realize our vision of a global array of autonomous floating sensors for whole-earth seismic tomography. |
|||||

BibTeX:
@article{Simons+2009a, author = { |
|||||

Frederik J. Simons, Hawthorne, J.C. and Beggan, C.D. |
Efficient analysis and representation of geophysical processes using localized spherical basis functions | 2009 | Wavelets XIII, pp. 74460G | inproceedings | DOI URL |

Abstract: While many geological and geophysical processes such as the melting of icecaps, the magnetic expression of bodies emplaced in the Earth’s crust, or the surface displacement remaining after large earthquakes are spatially localized, many of these naturally admit spectral representations, or they may need to be extracted from data collected globally, e.g. by satellites that circumnavigate the Earth. Wavelets are often used to study such nonstationary processes. On the sphere, however, many of the known constructions are somewhat limited. And in particular, the notion of 'dilation' is hard to reconcile with the concept of a geological region with fixed boundaries being responsible for generating the signals to be analyzed. Here, we build on our previous work on localized spherical analysis using an approach that is firmly rooted in spherical harmonics. We construct, by quadratic optimization, a set of bandlimited functions that have the majority of their energy concentrated in an arbitrary subdomain of the unit sphere. The ‘spherical Slepian basis’ that results provides a convenient way for the analysis and representation of geophysical signals, as we show by example. We highlight the connections to sparsity by showing that many geophysical processes are sparse in the Slepian basis. |
|||||

BibTeX:
@inproceedings{Simons+2009b, author = { |
|||||

Frederik J. Simons |
Afloat on a sea of noise | 2009 | Planet Earth Magazine Vol. Winter, pp. 28–29 |
article | URL |

Abstract: Global geophysicists map the Earth’s interior using seismic waves. But until recently the oceans have been largely off-limits. Frederik Simons describes creating a floating instrument that he hopes will change that. |
|||||

BibTeX:
@article{Simons+2009c, author = { |
|||||

Frederik J. Simons |
Turning freshmen into scientists with field research and quantitative analysis of geoscientific data | 2010 | MATLAB Digest Academic Edition Vol. October, pp. 1–3 |
article | URL |

Abstract: For students who arrive at Princeton expecting their courses to be heavily dependent on textbook learning, the Freshman Seminar Earth’s Changing Surface and Climate is a welcome surprise. Professor Adam Maloof and I give students practical experience in geological and geophysical research. Our approach is rooted in the conviction that anyone with a basic curiosity can ask scientifically interesting questions. Moreover, with the right tools and advice, students can be scientists from the minute they walk into the lab. |
|||||

BibTeX:
@article{Simons+2010, author = { |
|||||

Frederik J. Simons and Wang, D.V. |
Spatiospectral concentration in the Cartesian plane | 2011 | Intern. J. Geomath. Vol. 2(1), pp. 1–36 |
article | DOI URL |

Abstract: We pose and solve the analogue of Slepian’s time–frequency concentration problem in the two-dimensional plane, for applications in the natural sciences. We determine an orthogonal family of strictly bandlimited functions that are optimally concentrated within a closed region of the plane, or, alternatively, of strictly spacelimited functions that are optimally concentrated in the Fourier domain. The Cartesian Slepian functions can be found by solving a Fredholm integral equation whose associated eigenvalues are a measure of the spatiospectral concentration. Both the spatial and spectral regions of concentration can, in principle, have arbitrary geometry. However, for practical applications of signal representation or spectral analysis such as exist in geophysics or astronomy, in physical space irregular shapes, and in spectral space symmetric domains will usually be preferred. When the concentration domains are circularly symmetric in both spaces, the Slepian functions are also eigenfunctions of a Sturm–Liouville operator, leading to special algorithms for this case, as is well known. Much like their one-dimensional and spherical counterparts with which we discuss them in a common framework, a basis of functions that are simultaneously spatially and spectrally localized on arbitrary Cartesian domains will be of great utility in many scientific disciplines, but especially in the geosciences. |
|||||

BibTeX:
@article{Simons+2011a, author = { |
|||||

Frederik J. Simons, Loris, I., Nolet, G., Daubechies, I.C., Voronin, S., Vetter, P.A., Charléty, J. and Vonesch, C. |
Solving or resolving global tomographic models with spherical wavelets, and the scale and sparsity of seismic heterogeneity | 2011 | Geophys. J. Int. Vol. 187(2), pp. 969–988 |
article | DOI URL |

Abstract: We propose a class of spherical wavelet bases for the analysis of geophysical models and for the tomographic inversion of global seismic data. Its multiresolution character allows for modelling with an effective spatial resolution that varies with position within the Earth. Our procedure is numerically efficient and can be implemented with parallel computing. We discuss two possible types of discrete wavelet transforms in the angular dimension of the cubed sphere. We describe benefits and drawbacks of these constructions and apply them to analyse the information in two published seismic wave speed models of the mantle, using the statistics of wavelet coefficients across scales. The localization and sparsity properties of wavelet bases allow finding a sparse solution to inverse problems by iterative minimization of a combination of the l2 norm of the data residuals and the l1 norm of the model wavelet coefficients. By validation with realistic synthetic experiments we illustrate the likely gains from our new approach in future inversions of finite-frequency seismic data. |
|||||

BibTeX:
@article{Simons+2011b, author = { |
|||||

Frederik J. Simons, Loris, I., Brevdo, E. and Daubechies, I.C. |
Wavelets and wavelet-like transforms on the sphere and their application to geophysical data inversion | 2011 | Wavelets and Sparsity XIV, pp. 81380X | inproceedings | DOI URL |

Abstract: Many flexible parameterizations exist to represent data on the sphere. In addition to the venerable spherical harmonics, we have the Slepian basis, harmonic splines, wavelets and wavelet-like Slepian frames. In this paper we focus on the latter two: spherical wavelets developed for geophysical applications on the cubed sphere, and the Slepian "tree", a new construction that combines a quadratic concentration measure with wavelet-like multiresolution. We discuss the basic features of these mathematical tools, and illustrate their applicability in parameterizing large-scale global geophysical (inverse) problems. |
|||||

BibTeX:
@inproceedings{Simons+2011c, author = { |
|||||

Frederik J. Simons and Olhede, S.C. |
Maximum-likelihood estimation of lithospheric flexural rigidity, initial-loading fraction, and load correlation, under isotropy | 2013 | Geophys. J. Int. Vol. 193(3), pp. 1300–1342 |
article | DOI URL |

Abstract: Topography and gravity are geophysical fields whose joint statistical structure derives from interface-loading processes modulated by the underlying mechanics of isostatic and flexural compensation in the shallow lithosphere. Under this dual statistical-mechanistic viewpoint an estimation problem can be formulated where the knowns are topography and gravity and the principal unknown the elastic flexural rigidity of the lithosphere. In the guise of an equivalent ``effective elastic thickness'', this important, geographically varying, structural parameter has been the subject of many interpretative studies, but precisely how well it is known or how best it can be found from the data, abundant nonetheless, has remained contentious and unresolved throughout the last few decades of dedicated study. The popular methods whereby admittance or coherence, both spectral measures of the relation between gravity and topography, are inverted for the flexural rigidity, have revealed themselves to have insufficient power to independently constrain both it and the additional unknown initial-loading fraction and load-correlation factors, respectively. Solving this extremely ill-posed inversion problem leads to non-uniqueness and is further complicated by practical considerations such as the choice of regularizing data tapers to render the analysis sufficiently selective both in the spatial and spectral domains. Here, we rewrite the problem in a form amenable to maximum-likelihood estimation theory, which we show yields unbiased, minimum-variance estimates of flexural rigidity, initial-loading fraction and load correlation, each of those separably resolved with little a posteriori correlation between their estimates. We are also able to separately characterize the isotropic spectral shape of the initial loading processes. Our procedure is well-posed and computationally tractable for the two-interface case. The resulting algorithm is validated by extensive simulations whose behavior is well matched by an analytical theory with numerous tests for its applicability to real-world data examples. |
|||||

BibTeX:
@article{Simons+2013, author = { |
|||||

Frederik J. Simons and Plattner, A. |
Scalar and vector Slepian functions, spherical signal estimation and spectral analysis | 2015 | Handbook of Geomathematics, pp. 2563–2608 | incollection | DOI URL |

Abstract: It is a well-known fact that mathematical functions that are timelimited (or spacelimited) cannot be simultaneously bandlimited (in frequency). Yet the finite precision of measurement and computation unavoidably bandlimits our observation and modeling of scientific data, and we often only have access to, or are only interested in, a study area that is temporally or spatially bounded. In the geosciences we may be interested in spectrally modeling a time series defined only on a certain interval, or we may want to characterize a specific geographical area observed using an effectively bandlimited measurement device. It is clear that analyzing and representing scientific data of this kind will be facilitated if a basis of functions can be found that are ``spatiospectrally'' concentrated, i.e. ``localized'' in both domains at the same time. Here, we give a theoretical overview of one particular approach to this ``concentration'' problem, as originally proposed for time series by Slepian and coworkers, in the 1960s. We show how this framework leads to practical algorithms and statistically performant methods for the analysis of signals and their power spectra in one and two dimensions, and, particularly for applications in the geosciences, for scalar and vectorial signals defined on the surface of a unit sphere. |
|||||

BibTeX:
@incollection{Simons+2015, author = { |
|||||

Frederik J. Simons, Verhelst, F. and Swennen, R. |
Quantitative characterization of coal by means of microfocal X-ray computed microtomography (CMT) and color image analysis (CIA) | 1997 | Intern. J. Coal Geol. Vol. 34(1–2), pp. 69–88 |
article | DOI URL |

Abstract: Microfocal X-ray computed microtomography (CMT) is a novel technique that produces three-dimensional maps of the distribution of the linear attenuation coefficient inside an object. In contrast to the more conventional medical computerized tomography (CT) systems, a microfocal X-ray source is used. This enables a far better spatial resolution. The linear attenuation coefficient or tomodensity is dependent on the physical density and the mineralogy of the object to be imaged, and on the energy of the radiation used. Earlier work by Verhelst et al. (1996) presented the results of the correlation of the tomodensities obtained from three-dimensional CT scans with two-dimensional data on the composition of a coal sample acquired with color image analysis (CIA), a camera technique. This analysis assumed the linear proportionality of the tomodensity with the real physical bulk density, which is true only for certain energy ranges. In this paper, we use CMT-devices for a similar correlation. For sake of comparison, the same core sample was used. First, new CIA-data on the surface composition along two profiles were sampled with greater detail (100 micron). These data are subjected to a geostatistical analysis to quantify the spatial dependence between the measurements. Second, CMT-tomograms were made, yielding spatial resolutions twice as high as medical CT. A multivariate correlation was carried out, and two improved (geo)statistical methods are suggested. The different energy range of the microfocal X-ray source compared to medical CT, however, produces some bias in the correlation of the tomodensities with the surface percentages of the constituents. We therefore suggest that the linear attenuation coefficient be treated as a separate unit. No attempt was made to translate the linear attenuation coefficient to the physical bulk density of the different constituents of coal (e.g. macerals). |
|||||

BibTeX:
@article{Simons+97, author = { |
|||||

Frederik J. Simons, Zielhuis, A. and van der Hilst, R.D. |
The deep structure of the Australian continent from surface-wave tomography | 1999 | Lithos Vol. 48, pp. 17–43 |
article | DOI URL |

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. |
|||||

BibTeX:
@article{Simons+99, author = { |
|||||

Frederik J. Simons |
Slepian functions and their use in signal estimation and spectral analysis | 2010 | Handbook of Geomathematics, pp. 891–923 | incollection | DOI URL |

Abstract: It is a well-known fact that mathematical functions that are timelimited (or spacelimited) cannot be simultaneously bandlimited (in frequency). Yet the finite precision of measurement and computation unavoidably bandlimits our observation and modeling scientific data, and we often only have access to, or are only interested in, a study area that is temporally or spatially bounded. In the geosciences, we may be interested in spectrally modeling a time series defined only on a certain interval, or we may want to characterize a specific geographical area observed using an effectively bandlimited measurement device. It is clear that analyzing and representing data of this kind will be facilitated if a basis of functions can be found that are “spatiospectrally” concentrated, i.e., “localized” in both domains at the same time. Here, we give a theoretical overview of one particular approach to this “concentration” problem, as originally proposed for time series by Slepian and coworkers, in the 1960s. We show how this framework leads to practical algorithms and statistically performant methods for the analysis of signals and their power spectra in one and two dimensions, and on the surface of a sphere. |
|||||

BibTeX:
@incollection{Simons2010, author = { |
|||||

Frederik J. Simons |
On Foundations of Seismology: Bringing Idealizations Down to Earth, by James R. Brown and M. A. Slawinski [BibTeX] |
2018 | The Leading Edge Vol. 37(3), pp. 232 |
article | DOI URL |

BibTeX:
@article{Simons2018, author = { |
|||||

Frederik J. Simons |
Waves and Rays in Seismology: Answers to Unasked Questions, by Michael A. Slawinski [BibTeX] |
2019 | The Leading Edge Vol. 38(5), pp. 406-407 |
article | DOI URL |

BibTeX:
@article{Simons2019, author = { |
|||||

Slobbe, D.C., Frederik J. Simons and Klees, R. |
The spherical Slepian basis as a means to obtain spectral consistency between mean sea level and the geoid | 2012 | J. Geod. Vol. 86(8), pp. 609–628 |
article | DOI URL |

Abstract: The mean dynamic topography (MDT) can be computed as the difference between the mean sea level (MSL) and a gravimetric geoid. This requires that both data sets are spectrally consistent. In practice, it is quite common that the resolution of the geoid data is less than the resolution of the MSL data, hence, the latter need to be low-pass filtered before the MDT is computed. For this purpose conventional low-pass filters are inadequate, failing in coastal regions where they run into the undefined MSL signal on the continents. In this paper, we consider the use of a bandlimited, spatially concentrated Slepian basis to obtain a low-resolution approximation of the MSL signal. We compute Slepian functions for the oceans and parts of the oceans and compare the performance of calculating the MDT via this approach with other methods, in particular the iterative spherical harmonic approach in combination with Gaussian low-pass filtering, and various modifications. Based on the numerical experiments, we conclude that none of these methods provide a low-resolution MSL approximation at the sub-decimetre level. In particular, we show that Slepian functions are not appropriate basis functions for this problem, and a Slepian representation of the low-resolution MSL signal suffers from broadband leakage. We also show that a meaningful definition of a low-resolution MSL over incomplete spherical domains involves orthogonal basis functions with additional properties that Slepian functions do not possess. A low-resolution MSL signal, spectrally consistent with a given geoid model, is obtained by a suitable truncation of the expansions of the MSL signal in terms of these orthogonal basis functions. We compute one of these sets of orthogonal basis functions using the Gram–Schmidt orthogonalization for spherical harmonics. For the oceans, we could construct an orthogonal basis only for resolutions equivalent to a spherical harmonic degree 36. The computation of a basis with a higher resolution fails due to inherent instabilities. Regularization reduces the instabilities but destroys the orthogonality and, therefore, provides unrealistic low-resolution MSL approximations. More research is needed to solve the instability problem, perhaps by finding a different orthogonal basis that avoids it altogether. |
|||||

BibTeX:
@article{Slobbe+2012, author = {D. Cornelis Slobbe and |
|||||

Sukhovich, A., Irisson, J.-O., Frederik J. Simons, Ogé, A., Hello, Y.M., Deschamps, A. and Nolet, G. |
Automatic discrimination of underwater acoustic signals generated by teleseismic P-waves: A probabilistic approach | 2011 | Geophys. Res. Lett. Vol. 38, pp. L18605 |
article | DOI URL |

Abstract: We propose a new probabilistic scheme for the automatic recognition of underwater acoustic signals generated by teleseismic P‐waves recorded by hydrophones in the ocean. The recognition of a given signal is based on the relative distribution of its power among different frequency bands. The signal’s power distribution is compared with a statistical model developed by analyzing relative power distributions of many signals of the same origin and a numerical criterion is calculated, which can serve as a measure of the probability for the signal to belong to the statistical model. Our recognition scheme was applied to 6‐month‐long continuous records of seven ocean bottom hydrophones (OBH) deployed in the Ligurian Sea. A maximum of 94% of all detectable teleseismic P‐waves recorded during the deployment of the OBHs were recognized correctly with no false recognitions. The proposed recognition method will be implemented in autonomous underwater robots dedicated to detect and transmit acoustic signals generated by teleseismic P‐waves. |
|||||

BibTeX:
@article{Sukhovich+2011, author = {Alexey Sukhovich and Jean-Olivier Irisson and |
|||||

Sukhovich, A., Bonnieux, S., Hello, Y., Irisson, J.-O., Frederik J. Simons and Nolet, G. |
Seismic monitoring in the oceans by autonomous floats | 2015 | Nature Commun. Vol. 6, pp. 8027 |
article | DOI URL |

Abstract: Our understanding of the internal dynamics of the Earth is largely based on images of seismic velocity variations in the mantle obtained with global tomography. However, our ability to image the mantle is severely hampered by a lack of seismic data collected in marine areas. Here we report observations made under different noise conditions (in the Mediterranean Sea, the Indian and Pacific Oceans) by a submarine floating seismograph, and show that such floats are able to fill the oceanic data gap. Depending on the ambient noise level, the floats can record between 35 and 63% of distant earthquakes with a moment magnitude MZ6.5. Even magnitudes <6.0 can be successfully observed under favourable noise conditions. The serendipitous recording of an earthquake swarm near the Indian Ocean triple junction enabled us to establish a threshold magnitude between 2.7 and 3.4 for local earthquakes in the noisiest of the three environments. |
|||||

BibTeX:
@article{Sukhovich+2015, author = {Alexey Sukhovich and Sébastien Bonnieux and Yann Hello and Jean-Olivier Irisson and |
|||||

Wang, L., Shum, C.K., Frederik J. Simons, Tapley, B.D. and Dai, C. |
Coseismic and postseismic deformation of the 2011 Tohoku-Oki earthquake constrained by GRACE gravimetry | 2012 | Geophys. Res. Lett. Vol. 39, pp. L07301 |
article | DOI URL |

Abstract: Spaceborne gravimetry data from the Gravity Recovery And Climate Experiment (GRACE) are processed using spatio-spectral Slepian localization analysis enabling the high-resolution detection of permanent gravity change associated with both coseismic and postseismic deformation resulting from the great 11 March 2011 Mw 9.0 Tohoku-Oki earthquake. The GRACE observations are then used in a geophysical inversion to estimate a new slip model containing both coseismic slip and after-slip. The GRACE estimated moment for the total slip, up to the end of July 2011 is estimated as (4.59±0.49) 10^22 N m, equivalent to a composite Mw of 9.07±0.65. If the moment for the Tohoku-Oki main shock is assumed to be 3.8 x 10^22 N m, the contribution from the after-slip is estimated to be 3.0x10^21–12.8x10^21 N m, in good agreement with a postseismic slip model inverted from GPS data. We conclude that GRACE data provide an independent constraint to quantify co- and post-seismic deformation for the Tohoku-Oki event. |
|||||

BibTeX:
@article{Wang+2012a, author = {Lei Wang and C. K. Shum and |
|||||

Wang, L., Shum, C.K., Frederik J. Simons, Tassara, A., Erkan, K., Jekeli, C., Braun, A., Kuo, C., Lee, H. and Yuan, D.-N. |
Coseismic slip of the 2010 Mw 8.8 Great Maule, Chile, earthquake quantified by the inversion of GRACE observations | 2012 | Earth Planet. Sci. Lett. Vol. 335–336, pp. 167–179 |
article | DOI URL |

Abstract: The 27 February 2010 Mw 8.8 Maule, Chile, earthquake ruptured over 500 km along a mature seismic gap between 341 ◦S and 381 ◦S — the Concepcion-Constitucion gap, where no large megathrust earthquakes had occurred since the 1835 Mw 8.5 event. Notable discrepancies exist in slip distribution and moment magnitude estimated by various models inverted using traditional observations such as teleseismic networks, coastal/river markers, tsunami sensors, Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR). We conduct a spatio-spectral localization analysis, based on Slepian basis functions, of data from Gravity Recovery And Climate Experiment (GRACE) to extract coseismic gravity change signals of the Maule earthquake with improved spatial resolution (350 km half-wavelength). Our results reveal discernible differences in the average slip between the GRACE observation and predictions from various coseismic models. The sensitivity analysis reveals that GRACE observation is sensitive to the size of the fault, but unable to separate depth and slip. Here we assume the depth of the fault is known, and simultaneously invert for the fault-plane area and the average slip using the simulated annealing algorithm. Our GRACE-inverted fault plane length and width are 429±6 km, 146 ± 5km, respectively. The estimated slip is 8.1±1.2 m, indicating that most of the strain accumulated since 1835 in the Concepcion-Constitucion gap was released by the 2010 Maule earthquake. |
|||||

BibTeX:
@article{Wang+2012b, author = {Lei Wang and C. K. Shum and |
|||||

Wieczorek, M.A. and Frederik J. Simons |
Localized spectral analysis on the sphere | 2005 | Geophys. J. Int. Vol. 162(3), pp. 655–675 |
article | DOI URL |

Abstract: It is often advantageous to investigate the relationship between two geophysical data sets in the spectral domain by calculating admittance and coherence functions. While there exist powerful Cartesian windowing techniques to estimate spatially localized (cross-)spectral properties, the inherent sphericity of planetary bodies sometimes necessitates an approach based in spherical coordinates. Direct localized spectral estimates on the sphere can be obtained by tapering, or multiplying the data by a suitable windowing function, and expanding the resultant field in spherical harmonics. The localization of a window in space and its spectral bandlimitation jointly determine the quality of the spatiospectral estimation. Two kinds of axisymmetric windows are here constructed that are ideally suited to this purpose: bandlimited functions that maximize their spatial energy within a cap of angular radius θ_0, and spacelimited functions that maximize their spectral power within a spherical harmonic bandwidth L. Both concentration criteria yield an eigenvalue problem that is solved by an orthogonal family of data tapers, and the properties of these windows depend almost entirely upon the space–bandwidth product N_0 = (L + 1)θ_0/π. The first N_0−1 windows are near perfectly concentrated, and the best-concentrated window approaches a lower bound imposed by a spherical uncertainty principle. In order to make robust localized estimates of the admittance and coherence spectra between two fields on the sphere, we propose a method analogous to Cartesian multitaper spectral analysis that uses our optimally concentrated data tapers. We show that the expectation of localized (cross-)power spectra calculated using our data tapers is nearly unbiased for stochastic processes when the input spectrum is white and when averages are made over all possible realizations of the random variables. In physical situations, only one realization of such a process will be available, but in this case, a weighted average of the spectra obtained using multiple data tapers well approximates the expected spectrum. While developed primarily to solve problems in planetary science, our method has applications in all areas of science that investigate spatiospectral relationships between data fields defined on a sphere. |
|||||

BibTeX:
@article{Wieczorek+2005, author = {Mark A. Wieczorek and |
|||||

Wieczorek, M.A. and Frederik J. Simons |
Minimum-variance spectral analysis on the sphere | 2007 | J. Fourier Anal. Appl. Vol. 13(6), pp. 665–692 |
article | DOI URL |

Abstract: 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. |
|||||

BibTeX:
@article{Wieczorek+2007, author = {Mark A. Wieczorek and |
|||||

Wu, F.V., Borisov, D., Frederik J. Simons and Williamson, P. |
Waveform inversion for shear velocity and attenuation via the spectral-element adjoint method | 2021 | SEG Tech. Prog. Expanded Abstracts, pp. XXX | inproceedings | |

Abstract: Intrinsic attenuation affects both amplitude and phase behavior of the seismic wavefield. Its three-dimensional distribution in the subsurface needs to be robustly estimated to make earth models accurate. Full Waveform Inversion (FWI) has become the technique of choice for recovering high-resolution velocity structure, but its ability to reliably recover attenuation variations remains challenging. Multistage and hierarchical approaches offer promise to control the delicate trade-off between resolving elastic and anelastic structure in inversions for the shear velocity VS and quality factor Qμ, in regimes with strong surface waves prone to cycle skipping. We present a procedure for their joint reconstruction via spectral-element adjoint-based full-waveform modeling. Starting from a checkerboard example, we develop a frequency and offset continua- tion strategy for onshore active-source acquisition. In the early stages of our algorithm, the inversion targets the refinement of the shear-wave speed model by gradually widening the frequency band and increasing the offset range. The shear quality factor is of major consideration only in the final stage of the inversion. Our procedure is robust, competitive, and practical. |
|||||

BibTeX:
@inproceedings{Wu+2021, author = {Fan V. Wu and Dmitry Borisov and |
|||||

Yuan, Y.O. and Frederik J. Simons |
Multiscale adjoint waveform-difference tomography using wavelets | 2014 | Geophysics Vol. 79(3), pp. WA79–WA95 |
article | DOI URL |

Abstract: Full-waveform seismic inversions based on minimizing the distance between observed and predicted seismograms are, in principle, able to yield better-resolved earth models than those minimizing misfits derived from traveltimes alone. Adjoint-based methods provide an efficient way of calculating the gradient of the misfit function via a sequence of forward-modeling steps, which, using spectral-element codes, can be carried out in realistically complex media. Convergence and stability of full-waveform-difference adjoint schemes are greatly improved when data and synthetics are progressively presented to the algorithms in a constructive multiscale approximation using a (bi)orthogonal wavelet transform. Wavelets provide the nonredundant spectral decomposition that paves the way for the inversion to proceed successively from long-wavelength fitting to detailed exploration of the phases in the seismogram. The choice of wavelet class and type, the initial depth of the multiscale decomposition, and the minimization algorithms used at every level continue to play crucial roles in our procedure, but adequate choices can be made that test successfully on 2C elastic seismograms generated in toy models, as well as in the industry-standard 2D Marmousi model. Although for simplicity our inversion ignored surface waves by prior tapering and filtered removal, those also appeared to be very well matched in the final model. |
|||||

BibTeX:
@article{Yuan+2014a, author = {Yanhua O. Yuan and |
|||||

Yuan, Y.O., Frederik J. Simons and Bozdağ, E. |
Full-waveform adjoint tomography in a multiscale perspective | 2014 | SEG Technical Program Expanded Abstracts, pp. 1194–1199 | inproceedings | DOI URL |

Abstract: We discuss an algorithm to regularize elastic waveform inversions using wavelet-based constructive approximations of the data, synthetic and observed, in models that evolve as part of a gradient-based iterative scheme relying on forward and adjoint modeling carried out with a spectral-element method. For an elastic Marmousi model we show how our wavelet-based multiscale waveform inversion proceeds successively from large to small scales in the seismograms, with a progressive increase of the complexity of the resulting model. We explore the sensitivity of surface waves in imaging shallow structure. To circumvent cycle skipping we designed an envelope-misfit function within a wavelet-multiscale framework. We test our approach in a toy model in preparation for inversions at full complexity. |
|||||

BibTeX:
@inproceedings{Yuan+2014b, author = {Yanhua O. Yuan and |
|||||

Yuan, Y.O., Frederik J. Simons and Bozdağ, E. |
Multiscale adjoint tomography for surface and body waves | 2015 | Geophysics Vol. 80(5), pp. R281–R302 |
article | DOI URL |

Abstract: We have developed a wavelet-multiscale adjoint scheme for the elastic full-waveform inversion of seismic data, including body waves (BWs) and surface waves (SWs). We start the inversion on the SW portion of the seismograms. To avoid cycle skipping and reduce the dependence on the initial model of these dispersive waves, we commence by minimizing an envelope-based misfit function. Subsequently, we proceed to the minimization of a waveform-difference (WD) metric applied to the SWs only. After that, we fit BWs and SWs indiscriminately using a WD misfit metric. In each of these three steps, we guide the iterative inversion through a sequence of nested subspace projections in a wavelet basis. SW analysis preserves a wealth of near-surface features that would be lost in conventional BW tomography. We used a toy model to illustrate the dispersive and cycle-skipping behavior of the SWs, and to introduce the two ways by which we combat the nonlinearity of waveform inversions involving SWs. The first is the wavelet-based multiscale character of the method, and the second the envelope-based misfit function. Next, we used an industry synthetic model to perform realistic numerical experiments to further develop a strategy for SW and joint SW as well as BW tomography. The effect of incorrect density information on wave-speed inversions was also evaluated. We ultimately formalize a flexible scheme for full-waveform inversion based on adjoint methods that includes BWs and SWs, and also considers P- and S-wave speeds, as well as density. Our method is applicable to waveform inversion in exploration geophysics, geotechnical engineering, regional, and global seismology. |
|||||

BibTeX:
@article{Yuan+2015, author = {Yanhua O. Yuan and |
|||||

Yuan, Y.O., Frederik J. Simons and Tromp, J. |
Double-difference adjoint tomography | 2016 | Geophys. J. Int. Vol. 206(3), pp. 1599–1618 |
article | DOI URL |

Abstract: We introduce a `double-difference' method for the inversion for seismic wavespeed structure by adjoint tomography. Differences between seismic observations and model predictions at individual stations may arise from factors other than structural heterogeneity, such as errors in the assumed source-time function, inaccurate timings, and systematic uncertainties. To alleviate the corresponding nonuniqueness in the inverse problem, we construct differential measurements between stations, thereby reducing the influence of the source signature and systematic errors. We minimize the discrepancy between observations and simulations in terms of the differential measurements made on station pairs. We show how to implement the double-difference concept in adjoint tomography, both theoretically and in practice. We compare the sensitivities of absolute and differential measurements. The former provide absolute information on structure along the ray paths between stations and sources, whereas the latter explain relative (and thus higher-resolution) structural variations in areas close to the stations. Whereas in conventional tomography a measurement made on a single earthquake-station pair provides very limited structural information, in double-difference tomography one earthquake can actually resolve significant details of the structure. The double-difference methodology can be incorporated into the usual adjoint tomography workflow by simply pairing up all conventional measurements; the computational cost of the necessary adjoint simulations is largely unaffected. Rather than adding to the computational burden, the inversion of double-difference measurements merely modifies the construction of the adjoint sources for data assimilation. |
|||||

BibTeX:
@article{Yuan+2016, author = {Yanhua O. Yuan and |
|||||

Yuan, Y.O., Bozdağ, E., Ciardelli, C., Gao, F. and Frederik J. Simons |
The exponentiated phase measurement, and objective-function hybridization for adjoint waveform tomography | 2020 | Geophys. J. Int. Vol. 221, pp. 1145–1164 |
article | DOI URL |

Abstract: Seismic tomography has arrived at the threshold of the era of big data. However, how to extract information optimally from every available time series remains a challenge; one that is directly related to the objective function chosen as a distance metric between observed and synthetic data. Time-domain cross-correlation and frequency-dependent multitaper traveltime measurements are generally tied to window selection algorithms in order to balance the amplitude differences between seismic phases. Even then, such measurements naturally favor the dominant signals within the chosen windows. Hence, it is difficult to select all usable portions of seismograms with any sort of optimality. As a consequence, information ends up being lost, in particular from scattered waves. In contrast, measurements based on instantaneous phase allow extracting information uniformly over the seismic records without requiring their segmentation. And yet, measuring instantaneous phase, like any other phase measurement, is impeded by phase wrapping. In this paper, we address this limitation by using a complex-valued phase representation that we call `exponentiated phase'. We demonstrate that the xponentiated phase is a good substitute for instantaneous-phase measurements. To assimilate as much information as possible from every eismogram while tackling the nonlinearity of inversion problems, we discuss a flexible hybrid approach to combine various objective functions in adjoint seismic tomography. We focus on those based on the exponentiated phase, to take into account relatively small-magnitude scattered waves; on multitaper measurements of selected surface waves; and on cross-correlation measurements on specific windows to select distinct body-wave arrivals. Guided by synthetic experiments, we discuss how exponentiated-phase, multitaper, and cross-correlation measurements, and their hybridization, affect tomographic results. Despite their use of multiple measurements, the computational cost to evaluate gradient kernels for the objective functions is scarcely affected, allowing for issues with data quality and measurement challenges to be simultaneously addressed efficiently. |
|||||

BibTeX:
@article{Yuan+2020, author = {Yanhua O. Yuan and Ebru Bozdağ and Caio |