The Collaborative Seismic Earth Model: Generation 1.

We present a general concept for evolutionary, collaborative, multiscale inversion of geophysical data, specifically applied to the construction of a first-generation Collaborative Seismic Earth Model. This is intended to address the limited resources of individual researchers and the often limited use of previously accumulated knowledge. Model evolution rests on a Bayesian updating scheme, simplified into a deterministic method that honors today's computational restrictions. The scheme is able to harness distributed human and computing power. It furthermore handles conflicting updates, as well as variable parameterizations of different model refinements or different inversion techniques. The first-generation Collaborative Seismic Earth Model comprises 12 refinements from full seismic waveform inversion, ranging from regional crustal- to continental-scale models. A global full-waveform inversion ensures that regional refinements translate into whole-Earth structure.

Introduction

Evolving Multiscale Data and Inversion

Modern geophysical data are often characterized by their distribution across multiple spatial scales, meaning that high‐density local or regional deployments are complemented by coarser global‐scale recordings. Examples include data from seismic (e.g., Díaz et al., 2009; Meltzer et al., 1999), electromagnetic (e.g., Chulliat et al., 2017), gravity (e.g., International Gravimetric Bureau [http://bgi.omp.obs-mip.fr]), or Global Positioning System (GPS) (e.g., Kreemer et al., 2014) measurements. While data volumes increase steadily, the actually usable amount of information is limited by the available human and computing power of individual research groups. It therefore becomes increasingly difficult to assimilate all data consistently into one model using modern inversion techniques. The problem may be overcome by an inversion framework that supports collaborative and evolutionary inversion, where prior knowledge from earlier data can be incorporated.

Though being only one class of data, seismic recordings exemplify many of these challenges. Seismic data used for tomography are key to understand the Earth's structure and dynamics (e.g., Bunge et al., 2003; Colli et al., 2017; Koelemeijer et al., 2015; Trampert et al., 2004), explore and monitor reservoirs and volcanoes (e.g., Koulakov et al., 2013; Prieux et al., 2013), characterize earthquakes (e.g., Chen et al., 2007; Hejrani et al., 2017; Liu et al., 2004), and predict the ground motion they induce (e.g., Graves et al., 2010; Lee et al., 2014). Performing multiscale tomography is essential, as failure to resolve crustal details may induce artificial anisotropy (Bozdağ & Trampert, 2008; Ferreira et al., 2010; Fichtner, Kennett et al., 2013) that may be misinterpreted as mantle flow. While the need to jointly resolve crust and mantle is well recognized, prior knowledge with a level of detail known from crustal studies (e.g., Chen et al., 2007; Tape et al., 2010) is hardly incorporated in (global) tomography (e.g., Bozdağ et al., 2016; Koelemeijer et al., 2015; Schaeffer & Lebedev, 2013). The combination of resource limitations and insufficient use of prior knowledge prevent tomography from achieving the ideal of a multiscale model that assimilates the actually available wealth of data. A side effect are discrepancies between models, especially in parameters with strong dependence on small‐scale structure, for example, anisotropy and attenuation (Dalton et al., 2008; Ferreira et al., 2010).

Objectives and Outline

The main objective of this work is the development of an inversion framework that addresses the challenges of evolving multiscale data by meeting the following design criteria: (a) Enable a collaborative and evolutionary model construction to harness distributed human and computing power. (b) Consistently incorporate prior knowledge, accumulated during previous model updates. (c) Allow for variable parameterizations of different model refinements and of potentially different inversion techniques.

In section 2 we present the methodological developments, using the example of a Collaborative Seismic Earth Model (CSEM) for illustration. This includes the construction of an initial model, a Bayesian updating scheme for evolutionary multiscale refinement, pragmatic simplifications demanded by current computational limitations, and the handling of apparently conflicting updates. Following the general developments, we more specifically consider the first‐generation CSEM in section 3, which comprises 12 refinements from full seismic waveform inversion, ranging from regional crustal to global scales.

Evolutionary, Collaborative, Multiscale Inversion

In the following, we present a scheme for evolutionary, collaborative, multiscale inversion. It is designed to ensure computational efficiency, easy handling, completeness and flexibility of the parametrization, a clean refinement procedure, and the ability to handle both successive and simultaneous updates. Though being per se application independent, the general development is exemplified by the CSEM, which serves as a large‐scale use case that allows us to identify bottlenecks. In this context, we abandoned an earlier prototype based on a whole‐Earth refinable tetrahedral mesh (Afanasiev et al., 2016). While elegant in theory, it was too difficult to handle in practice, already requiring supercomputing resources for trivial tasks such as visualization and the addition of a refinement region.

Initial Model

We begin with the construction of an initial model. A physical parameter of the initial model as a function of position x, for example, an elastic parameter, can be expressed in terms of N 0 basis functions :

These may include polynomials, spherical harmonics, blocks, or a combination of these. The initial model vector is the maximum‐likelihood model of the prior probability density function (pdf), p 0(m 0), the shape of which will be specified in section 2.3.

In the specific context of the CSEM, the initial model is designed to be conservative in the sense of not containing structure that seismic data cannot modify, while still matching our longest‐period data to within half a cycle. It is built upon Preliminary Reference Earth Model (PREM; Dziewoński & Anderson, 1981). Since PREM's 220‐km discontinuity has been suggested to be only a regional feature (e.g., Gu et al., 2001), it is replaced by a linear gradient, as shown in Figures S2 and S3 in the supporting information. Superimposed are the S velocity variations of S20RTS (Ritsema et al., 1999), to which P velocity variations are scaled (Ritsema & van Heijst, 2002). The crust of PREM is replaced by the crustal model of (Meier et al., 2007), derived from surface wave inversion. Images of the initial model are presented in Figures S5 to S16.

Since bandlimited data can only constrain smooth, effective versions of discontinuities (Capdeville et al., 2013), we implemented the initial Moho as a gradient over ∼5 km, derived from the original crustal model by Backus averaging (Backus, 1962). This choice ensures that the effective sharpness of the Moho and thickness of the crust can be simply adjusted by not making an explicit distinction between crust and mantle during the inversion.

Model Parametrization and Bayesian Regional Updating

From the initial model the CSEM evolves through regional refinements that can be described within a rigorous Bayesian framework. While computational limitations currently enforce pragmatic simplifications (section 2.3), the idealistic Bayesian approach defines possible future improvements and offers solutions to emerging issues, such as simultaneous updates within overlapping regions.

The construction of a regionally refined version m 1(x) of m 0(x) starts with the addition of N 1 basis functions to the initial model representation,

The coefficient vector m 1 of the new basis functions may, for example, represent smaller‐scale variations relative to the initial model. They are initially constrained by the prior p 0(m 1|m 0), which is conditioned on the initial model itself because the probability of newly added variations depends on the structure that is already present. The prior contains, for instance, parameter correlations derived from mineral‐physics arguments or information on the allowable variations relative to the initial model. Through the assimilation of additional regional data, d 1, we obtain a conditional likelihood function L 0(m 1|m 0). Combining prior and likelihood via Bayes' theorem, yields the conditional posterior for the variations m 1, with normalization factor k. While both likelihood and posterior depend on d 1, this is omitted in the interest of a more succinct notation, especially in later developments. The joint posterior for the initial model coefficients m 0 and the variations m 1 follows from the conjunction of the conditional posterior p(m 1|m 0) with the initial model prior p 0(m 0),

This process can be iterated using p(m 0,m 1) as new prior, adding basis functions with coefficients m 2 for further refinement with a new data set d 2, and then using Bayes' theorem and the conjunction of pdfs to obtain the next posterior or next state of the CSEM, p(m 0,m 1,m 2). The basis functions of different refinements need not be the same.

The geometric parametrization of the CSEM differs strongly from the earlier prototype, defined on a large, regionally refined, tetrahedral mesh (Afanasiev et al., 2016). Each regional refinement retains its genuine basis functions used during its construction. The CSEM is thus a library of continuously defined updates that are assembled when the elastic properties of the Earth at some position x are queried. A set of positions may then be used for visualization or as grid points in a numerical wave propagation solver. The main benefits of the library versus the large, single mesh are greater flexibility in removing and adding refinement regions, as well as vastly reduced computational requirements for visualization and the extraction of submodels.

Pragmatic Simplifications

To avoid the computationally unfeasible sampling of high‐dimensional model spaces, the pure Bayesian approach requires simplifications. The first of these consists in the assumption that all pdfs are Gaussian, described by a mean and a covariance. The mean, or equivalently the maximum‐likelihood model of p(m 1|m 0), is then not computed via Bayes' theorem, as in equation (3), but approximated by minimizing a quadratic misfit functional χ, starting from the maximum‐likelihood initial model . Repeating the process leads to the representation of the current state's maximum‐likelihood model in terms of a sum of successive regional refinements, as illustrated in Figure 1.

Figure 1

Technical implementation and updating of the CSEM. (a) Schematic illustration of successive refinements, , , …, being added to the smooth maximum‐likelihood initial model . Each update is treated as a new item in an ordered library, keeping the parametrization used during its construction. Different updates may be parametrized differently, that is, have different types of basis functions. (b) Concrete example of SV velocity, v sv, at 120‐km depth in the initial model (top) and in the generation‐1 model of the CSEM (bottom). Refinement regions visible in this view include Europe, Turkey, the Sea of Marmara region, the North Atlantic, the South Atlantic, and North America. CSEM = Collaborative Seismic Earth Model.

The Hessian H of χ evaluated at equals the inverse covariance, of the posterior p(m 1|m 0). While H cannot be computed or stored explicitly, computationally manageable approximations can be estimated (e.g., Bui‐Thanh et al., 2013; Fichtner & Trampert, 2011a; Fichtner & van Leeuwen, 2015). Furthermore, arbitrary Hessian‐vector products can be evaluated efficiently using second‐order adjoints (e.g., Fichtner & Trampert, 2011b; Santosa & Symes, 1988), which is sufficient for uncertainty analysis and gradient‐based optimization schemes.

Apparent Inconsistencies and Dependent Updates

In the simplified scheme, the coefficients of already existing basis functions, m 0,…,m , remain unchanged when a new regional refinement with coefficients m is added. Instead of honoring potential dependencies of previous on newly added coefficients, the approximation of the maximum‐likelihood refinement is based on the minimization of a misfit functional with coefficients as fixed initial model and not by finding the actual maximum of the joint posterior p(m 0,…,m ,m ). The approximation error of the maximum‐likelihood model may lead to inconsistencies, meaning that previously assimilated data sets, d cannot be explained adequately anymore.

The occurrence of simplification‐related inconsistencies would in principle require reiterations of previously incorporated refinements. Alternatively, the Gaussian approximation allows us to replace repeated iterations by the solution of a least squares problem. To illustrate the concept, we combine two sets of model parameters, and , for instance from successive updates, into a joint parameter vector . We then consider two updates with conditional priors and , and with likelihood functions and , Bayes' theorem from equation (3) modifies to

Within the simplified framework of section 2.3, the maximum‐likelihood models, and , and posterior covariances, and , of the individual posteriors, and , are determined by independent misfit minimizations. It follows that the joint maximum‐likelihood model, , of p(m 1|m 0) is the minimum of

Noting that and are the Hessians and from the independent inversions; can be approximated iteratively from equation (6) using efficiently computable Hessian‐vector products (e.g., Fichtner & Trampert, 2011b; Santosa & Symes, 1988). Rearranging equation (6) and substituting the minimum yields the joint (inverse) posterior covariance,

Equation (7) indicates that the action of on a model (vector), needed for resolution analysis and further optimization, is simply the sum of actions of the individual Hessians on the model.

Equations (6) and (7) remove an apparent inconsistency by finding the solution that optimally agrees with the independent updates. While illustrated for updates of the initial model, the equations straightforwardly generalize to later‐stage updates of the CSEM, using the concept introduced in section 2.2.

In the current CSEM, described in detail in section 3, no regional refinement has so far led to an apparent inconsistency that would require the actual execution of the above‐derived scheme. This indicates that individual regional refinements indeed tend to improve the CSEM as a whole, at least at the scales that we currently consider. A more detailed discussion of this aspect can be found in section 4.

CSEM Generation 1

Data

The first‐generation CSEM includes 12 partially overlapping subregions, refined using three‐component waveform data at variable minimum periods, T. These are, in the order in which they have been incorporated, Australasia (T = 40 s; Fichtner et al., 2010), Europe (T = 50 s; Fichtner, Trampert et al., 2013), the North Atlantic (T = 25 s; Rickers et al., 2013), Turkey (T = 10 s; Fichtner, Saygin et al., 2013), the South Atlantic (T = 55 s; Colli et al., 2013), the Western Mediterranean (T = 15 s; Fichtner & Villaseñor, 2015), the wider Japanese Islands region (T = 20 s; Simute et al., 2016), the Iberian Peninsula (T = 8 s), Japan (T = 15 s), North America (T = 30 s), Southeast Asia (T = 20 s), and Western Turkey (T = 8 s; Çubuk Sabuncu et al., 2017). For the Iberian Peninsula and Southeast Asia, we incorporated vertical‐component ambient noise correlations, working under the simplifying assumption that correlations are proportional to interstation Rayleigh wave Green's functions (e.g., Lobkis & Weaver, 2001; Wapenaar & Fokkema, 2006).

To ensure consistency of regional refinements that do not overlap with previous ones, we incorporate three‐component waveforms from 223 events recorded globally by 149 high‐quality broadband stations at T = 55 s. The global data, assimilated between the updates in the Japanese Islands and the Western Mediterranean, translate the effects of regional refinements into improved whole‐Earth structure. As illustrated in the supporting information Figure S4, the global‐scale update changes S velocities by less than 1% within the uppermost mantle, mostly in response to large regional‐scale updates in Eurasia and Australasia.

A data summary is presented in Figure 2 in the form of a surface ray density plot. Since full‐waveform inversion used for the updates exploits a wide range of wave types with different finite‐frequency wave paths, ray coverage is only a rough, though useful, proxy for coverage. It illustrates, for instance, that western Europe, Japan, Australasia, and North America are already well covered, whereas the Pacific, central Africa, and central Asia remain to be improved. In total, the current CSEM incorporates waveform data from 948 events (including virtual events of noise correlations) recorded at 8,975 stations. The total number of source‐receiver pairs is 160,797.

Figure 2

Data coverage proxy. Shown are minor‐arc great circles connecting earthquake epicenters (red stars) and noise correlation virtual sources (yellow stars) to stations (green triangles). Gray scale represents the minimum period, T, of the waveform data used within a refinement region. Pale blue is used for the global data set with T = 55 s. Marked refinement regions are Australasia (Aus), Europe (Eu), Iberian Peninsula (Ib), wider Japanese Islands (Jp1), Japan (Jp2), North America (NAm), North Atlantic (NAt), South Atlantic (SAt), Southeast Asia (SEA), Turkey (Tu), Western Mediterranean (WMe), and Western Turkey (WTu). Refinement regions Ib, JP2, and WTu are marked by a white box for better visibility. The total number of distinct sources (earthquakes plus virtual sources from noise correlations) is 948. With 8,975 distinct stations , they provide 160,797 distinct source‐station pairs.

Inversion

To ensure that 3‐D wave propagation and finite‐frequency effects are properly taken into account, all updates to the CSEM have been performed using full‐waveform inversion, following the approach of Fichtner et al. (2010), to which we refer for technical aspects outside the focus of this work. We quantified differences between observed and synthetic seismograms using time‐frequency phase misfits (Fichtner et al., 2008), which we reduced iteratively until observations were explained to within the errors or the iteration stalled. To solve the forward and adjoint equations, we employed the spectral‐element solvers SPECFEM (Komatitsch & Tromp, 2002a, 2002b) and SES3D (Fichtner et al., 2009; Gokhberg & Fichtner, 2016) at the global and regional scales, respectively. Region‐specific inversion details and resolution analyses can be found in the references given in Figure 3.

Figure 3

Distribution of SV velocity, v sv, at depths of 15 km (top), 70 km (middle), and 200 km (bottom). At 15‐km depth, gray scale roughly corresponds to mantle velocities (4.0–4.8 km/s), and colors from red over white to blue roughly indicate crustal velocities (2.7–4.0 km/s). Key to marked features in regional close‐ups: (a) a.KM: Kazdağ Massif, a.KVP: Kula Volcanic Province, and a.MM: Menderes Massif; (b) b.JB: Japan Basin, b.SB: Shikoku Basin, and b.TB: Tsushima Basin; (c) c.PB: Pannonian Basin and c.TTL: Tornquist‐Teisseyre Line; (d) d.GC: Gawler Craton, d.KB: Kimberley Block, d.NFB: New England Fold Belt, d.PB: Pilbara Block, and d.YB: Yilgarn Block; (e) e.ACS: Apennine‐Calabria Slab, e.AS: Alboran Slab, and e.LPB: Liguro‐Provena̧l Basin; (f) f.PS: Pacific Slab, f.PSS: Philippine Sea Slab, and f.UA: Ulleung Anomaly.

Since the initial model is not the result of an inversion but an assembly of independent models, a consistent prior p 0(m) was not available and thus had to be designed. For simplicity, we set p 0(m) to be uniform. For all subsequent priors, p 0(m |m 0,…,m 1), we chose a covariance that enforces spatial correlation length scales approximately equal to the minimum wavelength. This common choice prevents the appearance of artifacts that are unlikely to be resolved. In all cases, we did not include correlations between physical model parameters.

Though the CSEM has flexibility to represent variations in any viscoelastic parameter, we restricted the inversions to parameters that are well resolved by the intermediate‐period waveform data, that is, P, SH, and SV velocities. Furthermore, we allowed for 3‐D density variations in order to avoid artifacts (Blom et al., 2017). Within the lower mantle, the CSEM is practically isotropic, though this has not been enforced explicitly.

Structural Model

A selection of slices through the SV velocity distribution of the current CSEM are presented in Figure 3. A more comprehensive gallery is shown in the supporting information Figures S5 to S16. Since a geologic interpretation is not among the objectives of this work, these are mostly intended to provide an impression of the multiscale nature of the CSEM. The images illustrate the large variation in length scales, ranging from ∼10 km in the crust of Western Turkey to ∼3,000 km in the Central Pacific upper mantle. Most of the significant deviations relative to the initial model are confined to the upper 500 km, that is, approximately the maximum penetration depth of regional long‐period surface and body waves.

To roughly separate crust and mantle visually, the 15‐km slice is shown in a mixed color scale, with gray scale above and colors below 4.0 km/s. The SV velocity distribution is marked by the outline of the continents and the smooth lateral variations of the initial crustal model (Meier et al., 2007). More detailed structure appears where regional‐scale shorter‐period data (<15 s) have been assimilated, that is, in western Europe, Turkey, and Japan.

From 70‐ to 200‐km depth, the model mostly represents upper‐mantle structure, an exception being the Tibetan Plateau. A general feature are anomalous SV velocities below −15%, compared to around −10% in S20RTS (Ritsema et al., 1999). In contrast, the largest positive SV variations of ∼10% are similar to those in S20RTS. While a clear causal relationship between anomalously low velocities, data, and inversion scheme is hard to establish, the iterative improvements of full‐waveform inversion and its ability to account for wavefront healing are likely to contribute (e.g., Igel & Gudmundsson, 1997; Malcolm & Trampert, 2011; Wielandt, 1987).

Discussion and Conclusions

The CSEM is not intended to be a long‐wavelength 3‐D reference for tomography or a community model in the sense of the Southern California Earthquake Center (SCEC) models (e.g., Small et al., 2017). Instead, it is designed as a framework for successive refinements on all scales and for the consistent use of prior knowledge. As such, it primarily tries to respond to current data and computing challenges, without claiming to be the only possible solution.

The updating scheme from section 2.2 is the foundation upon which the CSEM operates. Its generality makes it translatable to other domains, for example, the inversion of multiscale gravity, electromagnetic, or Global Positioning System (GPS) data. The current need to simplify suggests future improvements, especially in uncertainty quantification. Though the framework enables the incorporation of prior knowledge already at the level of the initial model, we decided not to do so, because previously constructed models usually have different parameterizations or lack sufficiently complete uncertainty information. Thus, the ability of the CSEM to take advantage of prior knowledge is to be understood in the sense of an update using information from all previous updates.

With the exception of the global update, the regional updates included in the current CSEM are sufficiently independent in coverage and frequency content to avoid the need for reiterations or the application of the scheme from section 2.3.1 in order to ensure consistency. Though this need not be the case for future refinements, our experience from global‐scale updates indicates that only one or two iterations are required to account for regional updates, thus keeping the computational costs at a reasonable level.

In this context, we note that the iterative refinement of the CSEM, including potential reiterations, constitutes an implementation of a stochastic gradient method. Instead of updating with a complete misfit function including all data (also from future deployments), only a quasi‐random subset of the data is used within an updating cycle. Since convergence of stochastic descent has been shown theoretically and empirically (e.g., Díaz & Guitton, 2011; Kiwiel, 2001), we expect that the CSEM should in fact converge towards a global optimum as new updates are added in a quasi‐random fashion, provided of course that cycle skipping issues continue to be carefully avoided.

Though all updates have so far been based on intermediate‐period waveform data, the generality of the framework allows for the incorporation of other data types, including derived data. The combination of the current full‐waveform updates with other seismic data types and inversion methods, for example, traveltime ray tomography, is work in progress.

To facilitate the participation of an increasing number of collaborators, the addition of refinements requires further automation. This, combined with the prospect of having CSEM hosted by the European Plate Observing System, will open the model to a wider community.

While multiscale tomography has been applied for at least two decades (e.g., Bijwaard et al., 1998; Karason & van der Hilst, 2000), the CSEM constitutes, to the best of our knowledge, the first collaborative and evolutionary framework for the construction of a multiscale seismic model. Despite imperfections of the early‐stage implementation, the CSEM is operational, and further refinements are in progress.

Supporting Information S1

Click here for additional data file.