Search

We present a novel inverse method for discriminating regional deformation and long-term fault creep by inversion of GNSS velocities observed at the spatial scale of intraplate faults by exploiting the different spatial signatures of these two mechanisms. In doing so our method provides a refined estimate of the upper bound of the strain accumulation process. As case study, we apply this method to a six year GNSS campaign (2003–2008) set up in the southern portion of the Pollino Range over the Castrovillari and Pollino faults. We show that regional deformation alone cannot explain the observed deformation pattern and implies high geodetic strain rate, with a WSW-ENE extension of 86±41×10−9/yr. Allowing for the possibility of fault creep, the modelling of GNSS velocities is consistent with their uncertainties and they are mainly explained by a shallow creep over the Pollino fault, with a normal/strike-slip mechanism up to 5 mm/yr. The regional strain rate decrease by about 70 percent and is characterized by WNW-ESE extension of 24±28×10−9/yr. The large uncertainties affecting our estimate of regional strain rate do not allow infering whether the tectonic regime of the area is extensional or strike-slip, although the latter is slightly more likely.

To monitor the stability of a mountain slope in northern Italy, microseismic monitoring technique has been used since 2013. Locating microseismic events is a basic step of this technique. We performed a seismic tomographic survey on the mountain surface above the rock face to obtain a reliable velocity distribution in the rock mass for the localization procedure. Seismic travel-time inversion showed high heterogeneity of the rock mass with strong contrast in velocity distribution. Low velocities were found at shallow depth on the top of the rock cliff and intermediate velocities were observed in the most critical area of the rock face corresponding to a partially detached pillar. Using the 3D velocity model obtained from inversion, localization tests were performed based on the Equal Differential Time (EDT) localization method. The results showed hypocenter misfits to be around 15 m for the five geophones of the microseismic network and the error was significantly decreased compared to the results produced by a constant velocity model. Although the localization errors are relatively large, the accuracy is sufficient to distinguish microseismic events occurring in the most critical zone of the monitored rock mass from microseismic events generated far away. Thus, the 3D velocity model will be used in future studies to improve the classification of the recorded events.

A key initial step in geophysical imaging is to devise an effective means of mapping the sensitivity of an observation to the model parameters, that is to compute its Fréchet derivatives or sensitivity kernel. In the absence of any simplifying assumptions and when faced with a large number of free parameters, the adjoint method can be an effective and efficient approach to calculating Fréchet derivatives and requires just two numerical simulations. In the Glacial Isostatic Adjustment problem, these consist of a forward simulation driven by changes in ice mass and an adjoint simulation driven by fictitious loads that are applied at the observation sites. The theoretical basis for this approach has seen considerable development over the last decade. Here, we present the final elements needed to image 3-D mantle viscosity using a dataset of palaeo sea-level observations. Developments include the calculation of viscosity Fréchet derivatives (i.e. sensitivity kernels) for relative sea-level observations, a modification to the numerical implementation of the forward and adjoint problem that permits application to 3-D viscosity structure, and a recalibration of initial sea level that ensures the forward simulation honours present-day topography. In the process of addressing these items, we build intuition concerning how absolute sea-level and relative sea-level observations sense Earth's viscosity structure and the physical processes involved. We discuss examples for potential observations located in the near field (Andenes, Norway), far field (Seychelles), and edge of the forebulge of the Laurentide ice sheet (Barbados). Examination of these kernels: (1) reveals why 1-D estimates of mantle viscosity from far-field relative sea-level observations can be biased; (2) hints at why an appropriate differential relative sea-level observation can provide a better constraint on local mantle viscosity and (3) demonstrates that sea-level observations have non-negligible 3-D sensitivity to deep mantle viscosity structure, which is counter to the intuition gained from 1-D radial viscosity Fréchet derivatives. Finally, we explore the influence of lateral variations in viscosity on relative sea-level observations in the Amundsen Sea Embayment and at Barbados. These predictions are based on a new global 3-D viscosity inference derived from the shear-wave speeds of GLAD-M25 and an inverse calibration scheme that ensures compatibility with certain fundamental geophysical observations. Use of the 3-D viscosity inference leads to: (1) generally greater complexity within the kernel; (2) an increase in sensitivity and presence of shorter length-scale features within lower viscosity regions; (3) a zeroing out of the sensitivity kernel within high-viscosity regions where elastic deformation dominates and (4) shifting of sensitivity at a given depth towards distal regions of weaker viscosity. The tools and intuition built here provide the necessary framework to explore inversions for 3-D mantle viscosity based on palaeo sea-level data., (© The Author(s) 2023. Published by Oxford University Press on behalf of The Royal Astronomical Society.)

Abstract In Engquist et al. (2016) [8] , the Wasserstein metric was successfully introduced to the full waveform inversion. We apply this method to the earthquake location problem. For this problem, the seismic stations are far from each other. Thus, the trace by trace comparison (Yang et al. [47]) is a natural way to compare the earthquake signals. Under this framework, we have derived a concise analytic expression of the Frèchet gradient of the Wasserstein metric, which leads to a simple and efficient implementation of the adjoint method. We square and normalize the earthquake signals for comparison so that the convexity of the misfit function with respect to earthquake hypocenter and origin time can be realized and observed numerically. To reduce the impact of noise, which does not offset each other after the signals are squared, a new control parameter is introduced. Finally, the LMF (Levenberg–Marquardt–Fletcher) method is applied to solve the resulted optimization problem. According to the numerical experiments, only a few iterations are required to converge to the real earthquake hypocenter and origin time. Even for data with noise, we can obtain reasonable and convergent numerical results. Highlights • Apply the Wasserstein metric to the Earthquake location problem. • Derive a concise analytic expression of the Frèchet gradient. • Present numerical examples to show the performance of the method. [ABSTRACT FROM AUTHOR]

Most linear inverse problems require regularization to ensure that robust and meaningful solutions can be found. Typically, Tikhonov-style regularization is used, whereby a preference is expressed for models that are somehow 'small' and/or 'smooth'. The strength of such preferences is expressed through one or more damping parameters, which control the character of the solution, and which must be set by the user. However, identifying appropriate values is often regarded as a matter of art, guided by various heuristics. As a result, such choices have often been the source of controversy and concern. By treating these as hyperparameters within a hierarchical Bayesian framework, we are able to obtain solutions that encompass the range of permissible regularization parameters. Furthermore, we show that these solutions are often well-approximated by those obtained via standard analysis using certain regularization choices which are—in a certain sense—optimal. We obtain algorithms for determining these optimal values in various cases of common interest, and show that they generate solutions with a number of attractive properties. A reference implementation of these algorithms, written in Python, accompanies this paper. [ABSTRACT FROM AUTHOR]

Time gating is a commonly used approach in the pre-processing of nuclear magnetic resonance data before inversion. Gating suppresses spurious signals that can degrade recovered decay time distributions and therefore often stabilizes inversion. However, care must be taken in applying this technique when assumptions of uncorrelated Gaussianity break down and reduce stacking efficiency. If not properly accounted for, unreliable noise estimates introduce inversion artefacts through over- or underfitting of the data. Block-bootstrap resampling of noise realization proxies obtained through data phasing can be used to generate reliable noise estimates for the windows. Benefits of the approach are demonstrated through inversion of synthetic and borehole data. Analysis confirms that bootstrapped noise metrics are more reliable under variable noise conditions and result in more stable inversion results. [ABSTRACT FROM AUTHOR]

The characterization of seismic source mechanisms and geometries is of critical importance for understanding the underlying causes and physical processes of low-frequency seismicity at volcanoes. Because observational data are often limited in volcanic environments by logistical constraints, we use synthetic modelling to investigate the capability of seismic networks to properly resolve source mechanisms. For 16 synthetic networks with as many as 40 stations, and variable distance and azimuthal distributions, we perform nonlinear moment-tensor inversion for six input source models. Using a grid search for source type and constrained waveform inversions provides a quantitative measure of source mechanism reliability. If the source location is assumed to be correct, results suggest that complete azimuthal coverage, with stations located at different distance ranges from the source will allow for a higher resolution recovery of the source-time function. In general, a similar degree of uncertainty characterizes configurations with as many as 40 stations and as few as eight stations. Although the level of uncertainty in the source-time function increases when fewer than eight stations are used in the inversions, sources are still recoverable when as few as four stations are used. Deviations from this general trend are present across the different input source models, with performance of configurations with fewer than eight stations showing a strong dependence on the source type. [ABSTRACT FROM AUTHOR]

While single point source is an oversimplified representation of medium to large earthquakes, finite fault models in many cases over parameterize the inversion due to the lack of sufficient near field data. Multiple point source solutions can fill in the gap in-between these two representations. Here, we propose a Markov-Chain-Monte-Carlo multiple point source inversion scheme, in combination with the advantage of the cut-and-paste technique, which cuts the seismogram into Pnl and surface portions and allows different time-shifts for each segment to align data with the synthetics. We apply the approach to the M w 6.2 foreshock in the 2016 Kumamoto earthquake sequence by using the strong-motion observations within 100 km. We are able to perform the inversion at relatively high-frequency ranges (0.02–0.28 Hz for Pnl and 0.02–0.20 Hz for surface waves) with confidence on the velocity model built on the M w 5.4 path calibration event. Our results show that the rupture was mainly composed of three subevents, with total duration of about 12 s and total M w of 6.2. The strikes of three subevents agree well with surface fault mapping where the Futagawa Fault intersects the Hinagu Fault with ∼30° difference in strike. Our solution shows that the first subevent dips to the southeast, while the second and the third subevents, located ∼3 km the north and ∼4 km to the southwest of the first subevent, dip to the northwest. The focal mechanism of the first subevent shows remarkable agreement with the first-motion solution. The fault geometry also shows well consistence with the relocated aftershocks, which delineate a SE-dipping fault around the first subevent and two NW-dipping faults to the north and south, respectively, corresponding to the second and the third subevent. The sum of moment tensor of subevents that have not only different geometry but also rake angles shows strong Compensated Linear Vector Dipole (CLVD) component (48–50 per cent). With a local 1D crustal model, a full-moment-tensor inversion using regional long-period waveform data also detects strong CLVD component of this earthquake. In contrast, using the PREM model results in an almost pure double-couple solution. In short, we have precisely resolved the rupture process, the intricate fault geometry, and the strong CLVD component with strong-motion data. This highlights the importance of extracting the relatively high-frequency information from the waveform data and with accurate velocity model in seismic source analyses of large earthquakes. [ABSTRACT FROM AUTHOR]

Automatic arrival picking of certain seismic or microseismic phases has been studied for decades. However, automatic detection of continuous signal waveforms has been seldom addressed. In this paper, I propose a novel approach for automatically detecting the waveforms in the microseismic data. The waveform detection can be formulated into a classification-based machine learning (ML) problem, that is each data point in the microseismic record needs to be classified as either waveform or non-waveform. I use the classic K-means clustering based unsupervised machine learning algorithm to solve this problem. I use mean, power, and spectral centroid as the three features to help the machine to characterize each data point. Both synthetic and real microseismic data examples are used to demonstrate the feasibility of the proposed algorithm. Results show that the algorithm can help detect the dominant waveforms in the data in an effective and efficient manner. The automatically detected waveforms can help quickly obtain the microseismic imaging result using an amplitude-based reverse time migration method. [ABSTRACT FROM AUTHOR]

We formulate a full wavefield inversion (FWFI) for ambient seismic noise recorded by large and dense seismograph arrays. FWFI exploits the constraints on the gradients of the wavefield that array data inherently possess. We pose FWFI as a partial differential equation (PDE) constrained inverse problem resulting in a joint estimation of a reconstructed wavefield and the medium parameters. The inverse problem is solved by variable projection. We explicitly allow for non-unique solutions to the PDE system that is imposed as a constraint. The boundary conditions of the wavefield do not need to be specified, and can remain unknown. This makes the algorithm suitable for inverting observations of ambient seismic noise by dense arrays. The result is that the inverse problem for subsurface properties becomes insensitive to the character and distribution of the noise sources that excited the seismic wavefield. In principle the formulation holds equally for ambient noise wavefields and for wavefields excited by controlled sources. The theory is supported with examples in one dimension in the time domain, and in two dimensions in the frequency domain. The latter are of interest in the inversion of surface wave ambient noise for phase velocity maps. [ABSTRACT FROM AUTHOR]

A big potential lies in the quantitative integration of multiple geophysical measurements for what concerns more unique and robust inversion results, complementary sensitivity to geological features and enhanced resolution. The mechanisms that enable such integration typically rely on statistics where the multiphysics parameter values are related to each other through regression functions (rock physics) or through the shape of the parameter spatial distributions (structure). Such coupling operators are jointly minimized with the data misfit to obtain coupled parameter distributions. We explore the application of various coupling mechanisms to synthetic and real data comprising seismic and electromagnetic measurements acquired in complex geological conditions such as salt geology and complex near surface in desert environment. Two types of structure operators, consisting of the standard cross-gradient and a newly introduced summative gradient, together with rock-physics operators are tested and combined for velocity model reconstruction of salt overburden in a marine environment and in a complex near surface case. Results of seismic-EM joint minimization suggest that the summative gradient operator with the introduction of the sign of the gradient correlation provides a strong coupling mechanism that can become advantageous in the presence of noise-affected seismic data. The rock-physics coupling mechanism is extremely strong and its effectiveness depends primarily on the reliability of the rock-physics relation. A Bayesian approach for the rock-physics operator is introduced to balance the related uncertainties and it is successfully tested in a joint minimization scheme. The combined use of rock-physics and structure coupling operators provides the best results in synthetic and real data applications. The use and combination of various regularization operators, as described, provide a formidable toolbox for solving a wide variety of ill-posed and non-unique geophysical inverse problems. [ABSTRACT FROM AUTHOR]

Inversion of surface wave data for crustal and upper-mantle structure is a staple of passive seismology, particularly since the advent of techniques enabling surface wave dispersion (SWD) and Rayleigh wave ellipticity measurements from ambient noise. Recent development and application of transdimensional Bayesian (TB) seismic inversion offers an approach to quantify model parameter uncertainties and trade-offs with fewer assumptions than traditional methods. Using synthetic tests, we investigate choices in the implementation of TB for the inversion of SWD and Rayleigh wave ellipticity to constrain the structure of Earth's continental lithosphere. We focus on three aspects of the inversion: limitation of data sensitivity, assumed scaling among parameters (compressional wave speed, Vp, shear wave speed, Vs, density and radial anisotropy) and parametrization choices. We show that while surface wave data provide relatively strong constraints on the posterior distribution of Vs and, to a lesser extent, Vp, common parametrization choices can potentially bias structure estimates. This is particularly the case for radial anisotropy (ξ), due to the inability to distinguish variations of Vp and density from those of ξ. Inferred results therefore depend substantially on the parametrization and scaling choices. We illustrate how layered parametrizations can, in the TB framework, recover smoothly varying profiles, and quantify the number of layers recoverable at different levels of measurement uncertainty. Finally, we propose two types of model parametrization for TB inversion involving multiple types of parameters. We demonstrate that by implementing an independent parametrization for different physical quantities, we can avoid imposing identical model geometry across multiple types of model parameters, and obtain better model estimates with reduced trade-offs. We advocate for such a parametrization in TB inversion of radial anisotropy using surface wave data, and when targeting disparate Vp and Vs structures such as those associated with |$\alpha $| - |$\beta $| quartz transtion. [ABSTRACT FROM AUTHOR]

Natural remanent magnetization and self-demagnetization on high-susceptibility bodies are two important factors affecting magnetic data inversion. We propose a framework for the inversion and interpretation of magnetic anomalies, in which significant remanent magnetization and self-demagnetization are present simultaneously. The framework is based on the assumptions that the external applied field and internal self-demagnetization field are uniform and the deflection of self-demagnetization in the total magnetization direction is negligible. First, the magnetization vector distributions are obtained from magnetic data by estimating the magnetization direction, then inverting for the magnetization intensity distribution, using the inferred magnetization direction as a constraint. Based on a priori information about the Koenigsberger ratio derived from petrophysical measurements, the direction and intensity of the remanent magnetization are obtained. The self-demagnetization factor is then computed using the finite volume method. Finally, the true-susceptibility distribution is achieved by correcting for the self-demagnetization effect. The method is first applied to synthetic magnetic data produced by a prism-shaped source model that has significant remanent magnetization and high susceptibility. In a case study of the Daye iron-ore deposit, Hubei province, China, the true susceptibility and remanent magnetization are reconstructed. The remanence direction information reveals that local geological activities such as synclines and faults lead to changes in the remanence directions at different local deposits. [ABSTRACT FROM AUTHOR]

A spectral element method for modelling surface wave dispersion of Love and Rayleigh waves is presented. This method uses standard Gauss–Lobatto–Legendre polynomials coupled with a Gauss–Laguerre–Legendre element to represent a half-space, and in doing so, improves both the efficiency and accuracy of the calculation of phase and group velocities, particularly at lower frequencies. It is demonstrated that this method is able to directly represent 1-D earth models with smoothly varying structure, layered structure or combinations thereof. The method is both efficient and accurate, while the solution error can be tuned across a range of frequencies to balance the trade-off with computational cost. In addition, the adjoint technique is used to develop an efficient algorithm for the calculation of the gradient of arbitrary misfit functions with respect to earth model parameters, which is of prime interest in the inversion of surface waves. It is demonstrated that accurate gradients can be computed for misfit functions based on phase velocity, group velocity, and Rayleigh wave ellipticity. To demonstrate both the spectral element method and the adjoint method developed in this paper, two simulated inverse problems are presented using Love wave phase velocity observations and Rayleigh wave ellipticity observations. [ABSTRACT FROM AUTHOR]

Elastic reflection waveform inversion (ERWI) utilizes reflections to update the low and intermediate wavenumbers in the deeper part of elastic models and can provide good initial models for elastic full waveform inversion (EFWI). Although ERWI aims to mitigate the nonlinearity of inversion when starting from a poor initial model, it suffers from the cycle-skipping problem due to the objective function of waveform fitting. Building initial P - and S -wave velocity models for EFWI through elastic wave-equation reflection traveltime inversion (ERTI) would be effective and robust since traveltime information relates to the background model more linearly. However, the current implementations of acoustic traveltime inversion is not straightforward in elastic media due to the existence of S -wavefields. Wave mode decomposition, both on the recording surface and in the extrapolated wavefields, is important for ERTI. First, for seismic data with P -wave sources, the P / S separation of multicomponent seismograms isolates the PP and PS reflection events and thus make it possible to extract the event-to-event time-shifts of these isolated reflections through dynamic image warping (DIW). Then, we can use the traveltime residuals of PP and PS reflections to build the objective function for ERTI. Second, based on the investigation of the complicated reflection kernels in an elastic medium, we demonstrate the necessity of wave mode decomposition applied on the extrapolated elastic wavefields, to suppress the artefacts induced by the undesirable cross-correlations of the components in forward and back-propagated wavefields. Therefore, the decomposition of surface recording data and extrapolated wavefields guarantees the dominate contribution of the traveltime is included during the ERTI. Accordingly, we propose a two-stage method to first build the P -wave background velocity using the separated PP reflections and then build the S -wave background velocity using the separated PS reflections based on the well-recovered P -wave velocity model. A numerical example of the Sigsbee2A model shows the effectiveness of the proposed ERTI approach. [ABSTRACT FROM AUTHOR]

We jointly invert for magnetic and velocity fields at the core surface over the period 1997–2017, directly using ground-based observatory time-series and measurements from the CHAMP and Swarm satellites. Satellite data are reduced to the form of virtual observatory time-series distributed on a regular grid in space. Such a sequential storage helps incorporate voluminous modern magnetic data into a stochastic Kalman filter, whereby spatial constraints are incorporated based on a norm derived from statistics of a numerical geodynamo model. Our algorithm produces consistent solutions both in terms of the misfit to the data and the estimated posterior model uncertainties. We retrieve core flow features previously documented from the analysis of spherical harmonic field models, such as the eccentric anticyclonic gyre. We find enhanced diffusion patterns under both Indonesia and Africa. In contrast to a steady flow that is strong under the Atlantic hemisphere but very weak below the Pacific, interannual motions appear evenly distributed over the two hemispheres. Recovered interannual to decadal flow changes are predominantly symmetrical with respect to the equator outside the tangent cylinder. In contrast, under the Northern Pacific we find an intensification of a high latitude jet, but see no evidence for a corresponding feature in the Southern hemisphere. The largest flow accelerations that we isolate over the studied era are associated with meanders, attached to the equatorward meridional branch of the planetary gyre in the Eastern hemisphere, that are linked to the appearance of an eastward equatorial jet below the Western Pacific. [ABSTRACT FROM AUTHOR]

We propose two ensembles of geomagnetic field models spanning the last three millennia: COV-ARCH is calculated using all available archeological artefacts and volcanic lava flows; lake and marine sediment records are added to this data set to build COV-LAKE. Given the sparse distribution of archeomagnetic observations and their associated large uncertainties, the recovery of magnetic field changes from such data is an ill-posed inverse problem that requires assuming some prior information. This is usually performed by imposing arbitrary regularizations in space and time. Instead, we construct the prior knowledge entering the objective function to be minimized from spatial and temporal statistics of the geomagnetic field, as available from satellites, ground-based observatories and paleomagnetic measurements, and validated by numerical geodynamo simulations. Our approach relies on the projection of model coefficients onto temporal cross-covariance functions. We show with synthetic experiments that the dispersion within the ensemble of solutions provides a coherent measure of the model uncertainties. Gauss coefficients inverted from geophysical records compare satisfactorily with those deduced from the independent database built upon historical and observatory records. A posteriori model errors are reduced when incorporating sediment records; they logically increase towards decreasing length-scales, indicating that a partial information is available up to a spherical degree 4–5. Such models and their associated uncertainties are suited to be used as observations in geomagnetic data assimilation studies. Our results advocate for an approximately constant dipole decay since ≈1700 AD, preceded by an era (going back to 1000 AD) where the dipole trend is weak, possibly slightly positive. We observe in both hemispheres, at both low- and high-latitude, persistent patches over the past 3000 yr. We also confirm a westward drift of flux patches at the core–mantle boundary at a speed of about 0.20 to 0.25° yr−1. Despite the sparse data distribution in the southern hemisphere, the South Atlantic Anomaly appears in both ensembles of models around 1800 AD. A similar low-intensity event seems to have appeared below the Indian Ocean over 600–1400 AD. Both global models are in general good agreement with regional master curves, though filtering out some of the centennial oscillations. [ABSTRACT FROM AUTHOR]

The selection of earthquake focal mechanisms (FMs) for stress tensor inversion (STI) is commonly done on a spatial basis, that is, hypocentres. However, this selection approach may include data that are undesired, for example, by mixing events that are caused by different stress tensors when for the STI a single stress tensor is assumed. Due to the significant increase of FM data in the past decades, objective data-driven data selection is feasible, allowing more refined FM catalogues that avoid these issues and provide data weights for the STI routines. We present the application of angular classification with expectation-maximization (ACE) as a tool for data selection. ACE identifies clusters of FM without a priori information. The identified clusters can be used for the classification of the style-of-faulting and as weights of the FM data. We demonstrate that ACE effectively selects data that can be associated with a single stress tensor. Two application examples are given for weighted STI from South America. We use the resulting clusters and weights as a priori information for an STI for these regions and show that uncertainties of the stress tensor estimates are reduced significantly. [ABSTRACT FROM AUTHOR]

The goal of an investigation, scientific or otherwise, is usually to find answers to some specific set of questions about the state of nature: what is the seismic velocity structure? How likely is this volcano to erupt within a certain period? Does a subsurface reservoir contain resources of interest? Background research may reveal the existence of pertinent knowledge and information discovered previously, new data are normally acquired, and an inference problem is solved in order to answer the questions taking both all of the a priori information and the new data into account. Inverse theory, decision theory and the theory of experimental design provide methods to optimize the design of the investigation and to estimate results. However, those theories are normally set in the context of a particular model of the universe, with its particular parametrization. This requires the investigator to specify a priori a coherent utility (a function that describes the risks and rewards) of all possible outcomes under that parametrization. Quite commonly, the investigator may not be able to do this. Ideally an investigator would be able merely to pose a set of questions, define a set of constraints on the data types, acquisition costs and logistics, and provide a functional to relate the questions to any particular parameter space. Theory and methodology would then semi-autonomously drive the interrogation of the state of nature by optimally selecting one or more relevant models and parameter spaces, and designing, acquiring and analysing data, in order to best answer the questions. If necessary this could be done in a sequential or iterative manner, which potentially then involves changing the questions posed in each iteration based on both previous results and inspiration from the investigator. We present such a theory of interrogation in this paper. We review the relevant aspects of decision and design theory, and cast them in a framework where the investigator specifies a utility only at the level required by the general questions to be posed. Each model under consideration is then mapped into this utility space of possible answers. We then extend this framework to sequential investigations, where the outcome of each step may affect all aspects of the problem: the models entertained, the utilities and even the questions themselves. A variety of examples illustrates the generality of this method: an asset team investigating how best to exploit a subsurface reservoir, Monte Carlo sampling to estimate the Bayesian evidence for geophysical models, discriminating between different rock physics models of strain in laboratory deformation experiments, an organization sequentially assessing the effectiveness of its methods to evaluate subsurface assets, assessing whether subsurface CO2 storage should be promoted for climate change mitigation, and examples running through the text of seismic tomography, earthquake characterization and autonomous interplanetary robotic exploration. [ABSTRACT FROM AUTHOR]

Characterization of the top 10–50 m of the subseabed is key for landslide hazard assessment, offshore structure engineering design and underground gas-storage monitoring. In this paper, we present a methodology for the stochastic inversion of ultra-high-frequency (UHF, 0.2–4.0 kHz) pre-stack seismic reflection waveforms, designed to obtain a decimetric-resolution remote elastic characterization of the shallow sediments with minimal pre-processing and little a priori information. We use a genetic algorithm in which the space of possible solutions is sampled by explicitly decoupling the short and long wavelengths of the P -wave velocity model. This approach, combined with an objective function robust to cycle skipping, outperforms a conventional model parametrization when the ground-truth is offset from the centre of the search domain. The robust P -wave velocity model is used to precondition the width of the search range of the multiparameter elastic inversion, thereby improving the efficiency in high-dimensional parametrizations. Multiple independent runs provide a set of independent results from which the reproducibility of the solution can be estimated. In a real data set acquired in Finneidfjord, Norway, we also demonstrate the sensitivity of UHF seismic inversion to shallow subseabed anomalies that play a role in submarine slope stability. Thus, the methodology has the potential to become an important practical tool for marine ground model building in spatially heterogeneous areas, reducing the reliance on expensive and time-consuming coring campaigns for geohazard mitigation in marine areas. [ABSTRACT FROM AUTHOR]

Abstract: Surface waves in seismic noise correlation functions are routinely exploited for tomographic studies. However, the observation and proper interpretation of body waves is considerably more difficult. To understand body wave generation in particular and the emergence of noise correlation signals in general, we investigate the sensitivity of correlation functions with respect to global noise sources. Treating noise correlation functions as self‐consistent observables, without the enforcement of Green function retrieval, enables us to investigate the quantitative contribution of different noise sources in space to a signal in a specific time window. The source sensitivity kernel for a surface wave reveals the expected dominant and broad region of stationary phase and thereby confirms the robustness of noise tomography in its current state. In contrast, kernels for potential body waves are highly oscillatory and indicate significant contributions from cross terms, for example, from body wave‐surface wave interactions. While the short‐wavelength nature of the body wave kernel may in the future allow us to infer noise sources with great detail, the importance of cross terms combined with the natural heterogeneity of noise sources renders an interpretation of such an arrival as a pure body wave almost impossible. Furthermore, we show that even long‐wavelength Earth structure significantly distorts sensitivity kernels. It is therefore essential to incorporate 3‐D heterogeneities into high‐resolution imaging of ambient noise sources, especially at higher frequencies. [ABSTRACT FROM AUTHOR]

We introduce a new Bayesian inversion method that estimates the spatial distribution of geological facies from attributes of seismic data, by showing how the usual probabilistic inverse problem can be solved using an optimization framework while still providing full probabilistic results. Our mathematical model consists of seismic attributes as observed data, which are assumed to have been generated by the geological facies. The method infers the post-inversion (posterior) probability density of the facies plus some other unknown model parameters, from both the seismic attributes and geological prior information. Most previous research in this domain is based on the localized likelihoods assumption, whereby the seismic attributes at a location are assumed to depend on the facies only at that location. Such an assumption is unrealistic because of imperfect seismic data acquisition and processing, and fundamental limitations of seismic imaging methods. In this paper, we relax this assumption: we allow probabilistic dependence between seismic attributes at a location and the facies in any neighbourhood of that location through a spatial filter. We term such likelihoods quasi-localized. Exact Bayesian inference is impractical because it requires normalization of the posterior distribution which is intractable for large models and must be approximated. Stochastic sampling (e.g. by using Markov chain Monte Carlo) is the most commonly used approximate inference method but it is computationally expensive and detection of its convergence is often subjective and unreliable. We use the variational Bayes method which is a more efficient alternative that offers reliable detection of convergence. It achieves this by replacing the intractable posterior distribution by a tractable approximation. Inference can then be performed using the approximate distribution in an optimization framework, thus circumventing the need for sampling, while still providing probabilistic results. We show in a noisy synthetic example that the new method recovered the coefficients of the spatial filter with reasonable accuracy, and recovered the correct facies distribution. We also show that our method is robust against weak prior information and non-localized likelihoods, and that it outperforms previous methods which require likelihoods to be localized. Our method is computationally efficient, and is expected to be applicable to 3-D models of realistic size on modern computers without incurring any significant computational limitations. [ABSTRACT FROM AUTHOR]

Four solutions of the spherical gravitational curvature boundary-value problems can be exploited for the determination of the Earth's gravitational potential. In this paper we discuss the combination of simulated satellite gravitational curvatures, that is, components of the third-order gravitational tensor, by merging these solutions using the spectral combination method. For this purpose, integral estimators of biased- and unbiased-types are derived. In numerical studies, we investigate the performance of the developed mathematical models for the gravitational field modelling in the area of Central Europe based on simulated satellite measurements. First, we verify the correctness of the integral estimators for the spectral downward continuation by a closed-loop test. Estimated errors of the combined solution are about eight orders smaller than those from the individual solutions. Second, we perform a numerical experiment by considering the Gaussian noise with the standard deviation of 6.5 × 10−17 m−1 s−2 in the input data at the satellite altitude of 250 km above the mean Earth sphere. This value of standard deviation is equivalent to a signal-to-noise ratio of 10. Superior results with respect to the global geopotential model TIM-r5 (Brockmann et al. 2014) are obtained by the spectral downward continuation of the vertical–vertical–vertical component with the standard deviation of 2.104 m2 s−2, but the root mean square error is the largest and reaches 9.734 m2 s−2. Using the spectral combination of all gravitational curvatures the root mean square error is more than 400 times smaller but the standard deviation reaches 17.234 m2 s−2. The combination of more components decreases the root mean square error of the corresponding solutions while the standard deviations of the combined solutions do not improve as compared to the solution from the vertical–vertical–vertical component. The presented method represents a weight mean in the spectral domain that minimizes the root mean square error of the combined solutions and improves standard deviation of the solution based only on the least accurate components. [ABSTRACT FROM AUTHOR]

Source encoding techniques alleviate the computational burden of sequential-source full waveform inversion (FWI) by considering multiple sources simultaneously rather than independently. The reduced data volume requires fewer forward/adjoint simulations per nonlinear iteration. Applications of source-encoded FWI (SEFWI) have thus far focused on monoparameter acoustic inversion. We extend SEFWI to the multiparameter case with applications presented for elastic isotropic inversion. Estimating multiple parameters can be challenging as perturbations in different parameters can prompt similar responses in the data. We investigate the relationship between source encoding and parameter trade-off by examining the multiparameter source-encoded Hessian. Probing of the Hessian demonstrates the convergence of the expected source-encoded Hessian, to that of conventional FWI. The convergence implies that the parameter trade-off in SEFWI is comparable to that observed in FWI. A series of synthetic inversions are conducted to establish the feasibility of source-encoded multiparameter FWI. We demonstrate that SEFWI requires fewer overall simulations than FWI to achieve a target model error for a range of first-order optimization methods. An inversion for spatially inconsistent P -wave (α) and S -wave (β) velocity models corroborates the expectation of comparable parameter trade-off in SEFWI and FWI. The final example demonstrates a shortcoming of SEFWI when confronted with time-windowing in data-driven inversion schemes. The limitation is a consequence of the implicit fixed-spread acquisition assumption in SEFWI. Alternative objective functions, namely the normalized cross-correlation and L1 waveform misfit, do not enable SEFWI to overcome this limitation. [ABSTRACT FROM AUTHOR]

We present a method for calculating the derivatives of measurements of glacial isostatic adjustment (GIA) with respect to the viscosity structure of the Earth and the ice-sheet history. These derivatives, or kernels, quantify the linearized sensitivity of measurements to the underlying model parameters. The adjoint method is used to enable efficient calculation of theoretically exact sensitivity kernels within laterally heterogeneous earth models that can have a range of linear or nonlinear viscoelastic rheologies. We first present a new approach to calculate GIA in the time domain, which, in contrast to the more usual formulation in the Laplace domain, is well suited to continuously varying earth models and to the use of the adjoint method. Benchmarking results show excellent agreement between our formulation and previous methods. We illustrate the potential applications of the kernels calculated in this way through a range of numerical calculations relative to a spherically symmetric background model. The complex spatial patterns of the sensitivities are not intuitive, and this is the first time that such effects are quantified in an efficient and accurate manner. [ABSTRACT FROM AUTHOR]

The finite element method using an unstructured tetrahedral mesh is one of the most effective methods for the 3-D modelling of marine magnetotelluric data that are strongly affected by bathymetry, because it enables us to incorporate both small-scale and regional-scale bathymetry into a computational mesh with a practical number of elements. The authors applied a 3-D inversion scheme using mesh of this type to marine magnetotelluric problems for the first time and verified its applicability. Forward calculations for two bathymetry models demonstrated that the results obtained with an unstructured tetrahedral mesh are close to the reference solutions. To evaluate the forward calculation results, we developed a general TM-mode analytical formulation for a 2-D sinusoidal topography. Moreover, synthetic inversion test results confirmed that a 3-D inversion scheme with an unstructured tetrahedral mesh enables us to recover subseafloor resistivity structure properly even for a model including a land–sea boundary as well as seafloor undulations. The verified inversion scheme was subsequently applied to a set of marine magnetotelluric data observed around the Iheya North Knoll, the middle Okinawa Trough. 3-D modelling using a mesh with precise bathymetry demonstrated that the data observed around the Iheya North Knoll are strongly affected by bathymetry, especially by the sea-depth differences between the depression of the trough and the shallow East China Sea. The estimated resistivity structure under the knoll is characterized by a conductive surface layer underlain by a resistive layer. The conductive layer implies permeable pelagic/hemipelagic sediments, which are consistent with a previous seismological study. Furthermore, the conductive layer has a resistive part immediately below the knoll, which is regarded as the consolidated magma intrusion that formed the knoll. Furthermore, at a depth of 10 km, we found that the resistor underneath the knoll extends to the southeast, implying that subseafloor resistivity under the Volcanic Arc Migration Phenomenon area is more resistive than the surroundings due to the presence of consolidated magma. [ABSTRACT FROM AUTHOR]

An earthquake rupture process can be kinematically described by rupture velocity, duration and spatial extent. These key kinematic source parameters provide important constraints on earthquake physics and rupture dynamics. In particular, core questions in earthquake science can be addressed once these properties of small earthquakes are well resolved. However, these parameters of small earthquakes are poorly understood, often limited by available data sets and methodologies. The Incorporated Research Institutions for Seismology Community Wavefield Experiment in Oklahoma deployed ∼350 three-component nodal stations within 40 km2 for a month, offering an unprecedented opportunity to test new methodologies for resolving small earthquake finite source properties in high resolution. In this study, we demonstrate the power of the nodal data set to resolve the variations in the seismic wavefield over the focal sphere due to the finite source attributes of an M2 earthquake within the array. The dense coverage allows us to tightly constrain rupture area using the second moment method even for such a small earthquake. The M2 earthquake was a strike-slip event and unilaterally propagated towards the surface at 90 per cent local S -wave speed (2.93 km s−1). The earthquake lasted ∼0.019 s and ruptured L c ∼70 m and W c ∼45 m. With the resolved rupture area, the stress-drop of the earthquake is estimated as 7.3 MPa for M w 2.3. We demonstrate that the maximum and minimum bounds on rupture area are within a factor of two, much lower than typical stress-drop uncertainty, despite a suboptimal station distribution. The rupture properties suggest that there is little difference between the M2 Oklahoma earthquake and typical large earthquakes. The new three-component nodal systems have great potential for improving the resolution of studies of earthquake source properties. [ABSTRACT FROM AUTHOR]

In this paper we present a comparison study of the time-domain (TD) and frequency-domain (FD) spectral induced polarization (IP) methods in terms of acquisition time, data quality and spectral information retrieved from inversion. We collected TDIP and FDIP surface measurements on three profiles with identical electrode setups at two different field sites with different lithology. In addition, TDIP data were collected in two boreholes using the El-Log drilling technique in which apparent formation resistivity and chargeability values are measured during drilling using electrodes integrated within the stem auger. The TD and FD data were processed and inverted with similar approaches. The IP voltage decays and complex impedance spectra were inspected quadrupole by quadrupole for outliers. In the inversion, the apparent resistivity values and the full IP decays (in TD) or the full complex impedance spectra (in FD) for all quadrupoles were inverted simultaneously in 2-D using a modified parametrization of the Cole-Cole model. Furthermore, inversions of synthetic models mimicking the field situations were performed. The comparisons reveal that TD and FD results are comparable not only qualitatively, but also quantitatively. Furthermore, the surface inversions are in agreement with the borehole results and lithology, within the resolution capability of the surface measurements. In particular, comparable spectral IP parameters were retrieved for both sites with both measurement approaches, with relaxation times between 10−2 s and 1 s, suggesting a spectral coverage at least up to 100 Hz for both TD and FD measurements. However, for the employed data acquisition procedures and instrumentations, the TD measurements had an advantage in terms of acquisition speed and usable spectral acquisition range compared to the FD measurements. We conclude that the TD method, as the FD method, is a suitable tool to recover spectral IP information in the field, provided that the measurement procedures are planned accordingly and the subsurface IP response is in the spectral range covered by the measurements. [ABSTRACT FROM AUTHOR]

Mesozoic to Cenozoic continental rifting, breakup and spreading between North America and Greenland led to the opening, from south to north, of the Labrador Sea and eventually Baffin Bay between Baffin Island, northeast Canada and northwest Greenland. Baffin Bay lies at the northern limit of this extinct rift, transform and spreading system and remains largely underexplored. With the sparsity of existing crustal-scale geophysical investigations of Baffin Bay, regional potential field methods and quantitative deformation assessments based on plate reconstructions provide two means of examining Baffin Bay at the regional scale and drawing conclusions about its crustal structure, its rifting history and the role of pre-existing structures in its evolution. Despite the identification of extinct spreading axes and fracture zones based on gravity data, insights into the nature and structure of the underlying crust have only been gleaned from limited deep seismic experiments, mostly concentrated in the north and east where the continental shelf is shallower and wider. Baffin Bay is partially underlain by oceanic crust with zones of variable width of extended continental crust along its margins. 3-D gravity inversions, constrained by bathymetric and depth to basement constraints, have generated a range of 3-D crustal density models that collectively reveal an asymmetric distribution of extended continental crust, approximately 25–30 km thick, along the margins of Baffin Bay, with a wider zone on the Greenland margin. A zone of 5–13 km thick crust lies at the centre of Baffin Bay, with the thinnest crust (5 km thick) clearly aligning with Eocene spreading centres. The resolved crustal thicknesses are generally in agreement with available seismic constraints, with discrepancies mostly corresponding to zones of higher density lower crust along the Greenland margin and Nares Strait. Deformation modelling from independent plate reconstructions using GPlates of the rifted margins of Baffin Bay was performed to gauge the influence of original crustal thickness and the width of the deformation zone on the crustal thicknesses obtained from the gravity inversions. These results show the best match with the results from the gravity inversions for an original unstretched crustal thickness of 34–36 km, consistent with present-day crustal thicknesses derived from teleseismic studies beyond the likely continentward limits of rifting around the margins of Baffin Bay. The width of the deformation zone has only a minimal influence on the modelled crustal thicknesses if the zone is of sufficient width that edge effects do not interfere with the main modelled domain. [ABSTRACT FROM AUTHOR]

The ill-posed nature of earthquake source estimation derives from several factors including the quality and quantity of available observations and the fidelity of our forward theory. Observational errors are usually accounted for in the inversion process. Epistemic errors, which stem from our simplified description of the forward problem, are rarely dealt with despite their potential to bias the estimate of a source model. In this study, we explore the impact of uncertainties related to the choice of a fault geometry in source inversion problems. The geometry of a fault structure is generally reduced to a set of parameters, such as position, strike and dip, for one or a few planar fault segments. While some of these parameters can be solved for, more often they are fixed to an uncertain value. We propose a practical framework to address this limitation by following a previously implemented method exploring the impact of uncertainties on the elastic properties of our models. We develop a sensitivity analysis to small perturbations of fault dip and position. The uncertainties of our fixed fault geometry are included in the inverse problem under the formulation of the misfit covariance matrix that combines both prediction and observation uncertainties. We validate this approach with the simplified case of a fault that extends infinitely along strike, using both Bayesian and optimization formulations of a static slip inversion. If epistemic errors are ignored, predictions are overconfident in the data and slip parameters are not reliably estimated. In contrast, inclusion of uncertainties in fault geometry allows us to infer a robust posterior slip model. Epistemic uncertainties can be many orders of magnitude larger than observational errors for great earthquakes (M w > 8). Not accounting for uncertainties in fault geometry may partly explain observed shallow slip deficits for continental earthquakes. Similarly, ignoring the impact of epistemic errors can also bias estimates of near-surface slip and predictions of tsunamis induced by megathrust earthquakes. [ABSTRACT FROM AUTHOR]

With a finite recording aperture, a limited source spectrum and unbalanced illumination, traditional imaging methods are insufficient to generate satisfactory depth profiles with high resolution and high amplitude fidelity. This is because traditional migration uses the adjoint operator of the forward modelling rather than the inverse operator. We propose a least-squares migration approach based on the time-domain Gaussian beam summation, which helps to balance subsurface illumination and improve image resolution. Based on the Born approximation for the isotropic acoustic wave equation, we derive a linear time-domain Gaussian beam modelling operator, which significantly reduces computational costs in comparison with the spectral method. Then, we formulate the corresponding adjoint Gaussian beam migration, as the gradient of an L2-norm waveform misfit function. An L1-norm regularization is introduced to the inversion to enhance the robustness of least-squares migration, and an approximated diagonal Hessian is used as a pre-conditioner to speed convergence. Synthetic and field data examples demonstrate that the proposed approach improves imaging resolution and amplitude fidelity in comparison with traditional Gaussian beam migration. [ABSTRACT FROM AUTHOR]

Full wavefield inversion (FWI) is a powerful and elegant approach for seismic imaging that is on the way to becoming the method of choice when processing exploration or global seismic data. In the case of processing marine survey data, one may be tempted to assume that acoustic FWI is sufficient given that only pressure waves exist in the water layer. In this paper, we pose the question as to whether or not in theory-at least for a hard waterbottom case-it should be possible to resolve the shear modulus or S-wave velocity in a marine setting using large offset data. We, therefore, conduct numerical experiments with idealized marine data calculated with the elastic wave equation. We study two cases, FWI of data due to a diffractor model, and FWI of data due to a fault model.We find that at least in idealized situation, elastic FWI of hard waterbottom data is capable of resolving between the two Lamé parameters λ and µ. Another numerical experiment with a soft waterbottom layer gives the same result. In contrast, acoustic FWI of the synthetic elastic data results in a single image of the first Lam´e parameter λ which contains severe artefacts for diffraction data and notable artefacts for layer reflection data. Based on these results, it would appear that at least the inversions of large offset marine data should be fully elastic rather than acoustic, unless it has been demonstrated that for the specific case in question (offsets, model and water depth, practical issues such as soft sediment attenuation of shear waves or computational time), an acoustic-only inversion provides a reasonably good quality of image comparable to that of an elastic inversion. Further research with real data is required to determine the degree to which practical issues such as shear wave attenuation in soft sediments may affect this result. [ABSTRACT FROM AUTHOR]

Despite surface displacements observed by geodesy are linear combinations of slip at faults in an elastic medium, determining the spatial distribution of fault slip remains a ill-posed inverse problem. A widely used approach to circumvent the illness of the inversion is to add regularization constraints in terms of smoothing and/or damping so that the linear system becomes invertible. However, the choice of regularization parameters is often arbitrary, and sometimes leads to significantly different results. Furthermore, the resolution analysis is usually empirical and cannot be made independently of the regularization. The stochastic approach of inverse problems provides a rigorous framework where the a priori information about the searched parameters is combined with the observations in order to derive posterior probabilities of the unkown parameters. Here, I investigate an approach where the prior probability density function (pdf) is a multivariate Gaussian function, with single truncation to impose positivity of slip or double truncation to impose positivity and upper bounds on slip for interseismic modelling. I show that the joint posterior pdf is similar to the linear untruncated Gaussian case and can be expressed as a truncated multivariate normal (TMVN) distribution. The TMVN form can then be used to obtain semi-analytical formulae for the single, 2-D or n-D marginal pdf. The semi-analytical formula involves the product of a Gaussian by an integral term that can be evaluated using recent developments in TMVN probabilities calculations. Posterior mean and covariance can also be efficiently derived. I show that the maximum posterior (MAP) can be obtained using a non-negative least-squares algorithm for the single truncated case or using the bounded-variable least-squares algorithm for the double truncated case. I show that the case of independent uniform priors can be approximated using TMVN. The numerical equivalence to Bayesian inversions using Monte Carlo Markov chain (MCMC) sampling is shown for a synthetic example and a real case for interseismic modelling in Central Peru. The TMVN method overcomes several limitations of the Bayesian approach using MCMC sampling. First, the need of computer power is largely reduced. Second, unlike Bayesian MCMC-based approach, marginal pdf, mean, variance or covariance are obtained independently one from each other. Third, the probability and cumulative density functions can be obtained with any density of points. Finally, determining the MAP is extremely fast. [ABSTRACT FROM AUTHOR]

Accounting for lateral heterogeneities in the 3-D velocity structure of the crust is known to improve earthquake source inversion, compared to results based on 1-D velocity models which are routinely assumed to derive finite-fault slip models. The conventional approach to include known 3-D heterogeneity in source inversion involves pre-computing 3-D Green's functions, which requires a number of 3-D wave propagation simulations proportional to the number of stations or to the number of fault cells. The computational cost of such an approach is prohibitive for the dense data sets that could be provided by future earthquake observation systems. Here, we propose an adjoint-based optimization technique to invert for the spatiotemporal evolution of slip velocity. The approach does not require pre-computed Green's functions. The adjoint method provides the gradient of the cost function, which is used to improve the model iteratively employing an iterative gradient-based minimization method. The adjoint approach is shown to be computationally more efficient than the conventional approach based on pre-computed Green's functions in a broad range of situations. We consider data up to 1 Hz from a Haskell source scenario (a steady pulse-like rupture) on a vertical strike-slip fault embedded in an elastic 3-D heterogeneous velocity model. The velocity model comprises a uniform background and a 3-D stochastic perturbation with the von Karman correlation function. Source inversions based on the 3-D velocity model are performed for two different station configurations, a dense and a sparse network with 1 and 20 km station spacing, respectively. These reference inversions show that our inversion scheme adequately retrieves the rise time when the velocity model is exactly known, and illustrates how dense coverage improves the inference of peak-slip velocities. We investigate the effects of uncertainties in the velocity model by performing source inversions based on an incorrect, homogeneous velocity model. We find that, for velocity uncertainties that have standard deviation and correlation length typical of available 3-D crustal models, the inverted sources can be severely contaminated by spurious features even if the station density is high. When data from thousand or more receivers is used in source inversions in 3-D heterogeneous media, the computational cost of the method proposed in this work is at least two orders of magnitude lower than source inversion based on pre-computed Green's functions. [ABSTRACT FROM AUTHOR]

We present a methodology for 1-D imaging of upper-mantle structure using a Bayesian approach that incorporates a novel combination of seismic data types and an adaptive parametrization based on piecewise discontinuous splines. Our inversion algorithm lays the groundwork for improved seismic velocity models of the lithosphere and asthenosphere by harnessing the recent expansion of large seismic arrays and computational power alongside sophisticated data analysis. Careful processing of P- and S-wave arrivals isolates converted phases generated at velocity gradients between the mid-crust and 300 km depth. This data is allied with ambient noise and earthquake Rayleigh wave phase velocities to obtain detailed VS and VP velocity models. Synthetic tests demonstrate that converted phases are necessary to accurately constrain velocity gradients, and S--p phases are particularly important for resolving mantle structure, while surface waves are necessary for capturing absolute velocities. We apply the method to several stations in the northwest and north-central United States, finding that the imaged structure improves upon existing models by sharpening the vertical resolution of absolute velocity profiles, offering robust uncertainty estimates, and revealing mid-lithospheric velocity gradients indicative of thermochemical cratonic layering. This flexible method holds promise for increasingly detailed understanding of the upper mantle. [ABSTRACT FROM AUTHOR]

The empirical orthogonal function (EOF) reconstruction (EOFR) combines the dense observations, for example, satellite altimetry data or model output covering a short period, and sparse tide gauge records (TGRs) spanning a long period to infer the past sea level variability (SLV). Given that the number of the TGRs reduces significantly backward over time, it is necessary to assess how the sparse TGRs affect the reconstructed SLV. Meanwhile, EOFR involves two techniques, that is, using error matrix or not, which has not yet been fully assessed. We find that error matrix plays an important role in the EOFR. Using error matrix produces better reconstructed SLV than those without error matrix, especially when the TGRs are sparse (e.g. <100). If the error matrix is not included, the SLV reconstructed with sparse TGRs tends to be seriously overestimated. It is also found that global mean sea level variability reconstructed with sparse TGRs is not reliable even the error matrix is used. However, the recovery of first EOF pattern is shown to be robust through the entire 20th century (correlation >0.7). Since the first EOF is highly related to El Niño--Southern Oscillation (ENSO), this suggests that EOFR is useful for revealing the evolution of ENSO variability. [ABSTRACT FROM AUTHOR]