Tomographic analysis of Pn arrivals times—the guided P-wave propagating within the lithospheric mantle—is ideal for studying the structure of the uppermost mantle. While plate-scale seismic images of Pn wave speeds are common beneath the continents, similar scale studies have not been possible within ocean basins due to the sparsity of seismic stations. The Cascadia Initiative (CI) dataset provides the first opportunity to image spatial variations in lithospheric structure across an entire oceanic plate. We measure 2,862 Pn arrivals from local earthquakes recorded by the CI array. Our dataset provides complete coverage of both the Juan de Fuca (JdF) and deforming Gorda plates. We invert the measured arrival times for 3D variations in anisotropic P-wave velocity and hypocentral parameters. Despite surficial evidence of extensive active faulting, the velocity structure of the Gorda uppermost mantle is remarkably consistent with predictions from a conductive cooling model (attached figure). Limited deformation at mantle depths is supported by seismic anisotropy measurements that show the fast-direction of P-wave propagation rotates in concert with the magnetic anomaly lineations. This rotation may be explained by local plate kinematics without internal deformation and hydration of the shallow mantle. In contrast to Gorda, the seismic velocity structure of the JdF plate does not exhibit a clear age dependence. Three pronounced mantle low-velocity zones are found along the southern edge of the JdF plate near the termini of large pseudofaults that contributed to the formation of the Blanco Transform fault. We attribute these velocity reductions to mantle alteration by seawater. We note that within the interior of the JdF plate pseudofaults do not appear as uniformly slow features in our seismic images. Beneath the central and northern JdF plate, P-wave speeds are ~7.7 km/s out to ~4-5 Myr before abruptly increasing to ~7.9 km/s. Curiously, this transition occurs near the onset of mantle downwelling inferred from teleseismic body wave tomography and attenuation suggesting that mantle flow dynamics may influence the structure of young oceanic lithosphere. Lastly, we note that our results do not suggest a relationship between the structure of the uppermost slab mantle and segmentation of the Cascadia megathrust.
The 2018 eruption of Kīlauea Volcano included caldera collapse at the summit of the volcano that was well recorded by a network of both permanent and temporary seismometers and infrasound microphones. Volcanic activity prior to the start of the eruption was elevated, including high lava lake levels and increased seismicity. Deflation of the summit began shortly after the eruption in the lower East Rift Zone commenced and was accompanied by a drop in the summit lava lake, which eventually disappeared from view after lowering by several hundred meters. Continued volume loss from beneath the summit eventually led to caldera collapse, which was accompanied by increases in earthquake rates and tremor amplitudes. Infrasonic tremor that originally rose from spattering at the surface of the lava lake was replaced by discrete infrasonic events arising from rockfalls and slumps. On May 16, the first of tens of collapse events occurred, beginning a cycle of increasing earthquake rates and energy release in the presence of deflation. After an initial increase in rate, seismicity rates remained constant for several hours prior to the next collapse event. Earthquakes during these events were typically part of repeating earthquake families and occurred on one of several circumferential cracks that marked active collapsing blocks. Through time, earthquake activity has mirrored the morphology of the collapse, migrating primarily north and east as caldera down drop extended in those directions. Collapse-related infrasound arrivals were initially down, suggesting downdropping of the surface, followed several seconds later by higher frequency infrasound, which we interpret to reflect the explosive or expulsion phase reaching the surface. Infrasound signals were 2-3 times more energetic below 0.5 Hz than above. Especially after May 29, when caldera collapse became much larger in surface area, infrasound signals were highly repetitive with a strong downward first motion. The collapse events were also deficient at seismic frequencies > 0.5 Hz, compared to typical tectonic earthquake sources. Despite similarities in waveforms at low frequencies (<0.1 Hz) and in infrasound, the seismic waveforms at high frequencies were not similar, reflecting either a unique source in each collapse event or a change in location.
We develop a Hamiltonian Monte Carlo (HMC) sampler which solves a multi-parameter elastic full-waveform inversion (FWI) in a probabilistic setting for the first time. This gives novel access to the full posterior distribution for this type of highly non-linear inverse problem. Typically, FWI has focused on using gradient descent methods with proper regularization to iteratively update models to a minimum misfit value. Non-uniqueness and uncertainties are mostly in this approach. Bayesian inversions offer an alternative by assigning a probability to each model in model space given some data and prior constraints. The drawback is the need to evaluate a very large number of models. Random walks from Markov chains counter this effect by only exploring regions of model space where probability is significant. HMC method additionally incorporates gradient information, i.e. local structure, typically available for numerical waveform tomography experiments. So far, HMC has only been implented for acoustic FWI. We implement HMC for multiple 2D elastic FWI set-ups. Using parallelized wave propagation code, wavefields and kernels are computed on an regular numerical grid and projected onto basis functions. These gradients are subsequently used to explore the posterior space of different target models using HMC. The free parameters in these experiments are P and S velocity, and density. Although simulating Hamiltonian dynamics in the resulting phase space is approximated numerically, the results of the Markov chain are nevertheless very insightful. No prior tuning of kernels, data or model space is required, under the constraint that the sampler is properly tuned. After a burn in phase during which the mass matrix is iteratively optimized, the Markov chain is run on multiple nodes. After approximately 100,000 samples (combined from all nodes) the Markov chain mixes well. The resulting samples give access to the full posterior distribution, including the mean and maximum-likelihood models, conditional probabilities, inter-parameter correlations and marginal distributions.
The Earth’s crystalline inner core (IC) solidifies from the liquid Fe alloy of the outer core (OC), which releases latent heat and light elements sustaining the geodynamo. Variability in solidification regime at the inner core boundary (ICB) may result in compositional and thermal multi-scale mosaic of the IC surface and dissimilarity of its hemispheres. Both the mosaic and hemisphericity are poorly constrained, not least due to a lack of available sampling by short-period reflected waves. Measured amplitude ratio of seismic phases of PKiKP and PcP reflected, respectively, off the inner and outer boundary of the liquid core, yields direct estimate of the ICB density jump. This parameter is capable of constraining the inner – outer core compositional difference and latent energy release, but is not well known (0.2–1.2 g/cub cm), and its distribution is obscure. Travel time measurements of PKiKP and PcP waveforms can be useful in terms of getting an insight into fine structure of ICB and its topography. We analyse a new representative sample of pre-critical PKiKP/PcP differential travel times and amplitude ratios that probes the core’s spots under Southeastern Asia and South America . We observe a statistically significant systematic bias between the measurements collected in western and eastern hemispheres, and carefully examine its origin. Separating the effects of core-mantle boundary (CMB) and ICB on the measured differentials is particularly challenging and we acknowledge that a whole class of physically valid models involving D” heterogeneities and lateral variation in lower mantle attenuation can be addressed to account for the observed hemisphericity. However, we find that variance in PKiKP-PcP differential travel times measured above the epicentral distance of 16 degrees is essentially due to mantle heterogeneities. Analysis of data below this distance indicates the ICB density jump under Southeastern Asia can be about 0.3 g/cub. cm, which is three times as small as under South America where also the thickness of the liquid core can be by 1-3 km in excess of the one in the East. The findings are interpretable as evidence for IC hemispherical asymmetry whereby crystallization dominates in the West and melting in the East (and not vice versa) or in terms of two disconnected mosaic patches with contrasting properties.
Timing errors are a notorious problem in seismic data acquisition and processing. We present a technique that allows such time shifts to be detected and corrected in a systematic fashion. The technique relies on virtual-source surface-wave responses retrieved through the crosscorrelation of ambient seismic noise. In particular, it relies on the theoretical time-symmetry of these time-averaged receiver-receiver crosscorrelations. By comparing the arrival time of the surface waves at positive time to the arrival time of the surface waves at negative time for a large a number of receiver-receiver pairs, relative timing errors can be determined in a least-squared sense. The time-symmetry of the receiver-receiver crosscorrelations, however, is contingent on a uniform surface-wave (noise) illumination pattern. In practice, the illumination pattern is often not uniform. We therefore show that weighting different receiver-receiver pairs differently in the inversion allows timing errors to be determined more accurately. The weights are based on the susceptibility of different receiver pairs to illumination-related travel-time errors. The proposed methodology is validated using both synthetic data and field data. The field data consists of recordings of ambient seismic noise by an array of stations centered around the tip of the Reykjanes peninsula, southwest Iceland (some of these stations exhibit time shifts of an unknown nature).
Since 2008, the rate of seismic events within the Central United States has dramatically increased, which is likely associated with wastewater injection from nearby oil and gas operations. Surface deformation measurements derived from spaceborne interferometric synthetic aperture radar (InSAR) data can be used to quantify the magnitude and spatial extent of the injection-related stress perturbation, which are critical for understanding the complex interaction between the injected fluid and the earth’s subsurface. In this study, we processed Sentinel-1 InSAR data over Central and West Texas using a recently developed processing framework that performs topography/geometry phase corrections prior to the interferogram formation (Zebker 2017). We streamlined the creation of upsampled digital elevation maps (DEMs) from NASA Shuttle Radar Topographic Mission (SRTM) data, as well as the collection of Sentinel-1 precise orbit data. We developed a tool for InSAR time-series analysis and data visualization. To detect unknown deformation signatures from large volumes of InSAR data, we employed computer vision ideas for feature detection independent of scale, well known through their success in the Scale Invariant Feature Transform (SIFT). We used multi-scale Laplacian-of-Gaussian (LoG) filters to find local maxima and minima in a coarse deformation solution, corresponding to “bowls” of uplift and subsidence, respectively. This allowed us to drastically cut down processing time of high-resolution InSAR products. As a validation, our method successfully detected all sinkhole locations, injection-related uplift signals and production-related subsidence signals as reported in Kim and Lu (2017) over a 100 km x 100km search area without the need for manual inspection. We then examined the Dallas Fort Worth Basin area for evidence of deformation near wastewater injection and oil/gas production sites. We begin to quantify the uncertainty from common noise sources to produce more confident time-series results.
Measurements and models of global geophysical parameters such as potential fields, seismic velocity models and dynamic / residual topography are well represented as 2D coloured / contoured maps. However, as teaching aids and for outreach, they offer little impact. Modern 3D-printing techniques help to visualize these and other concepts that are difficult to grasp, such as the intangible structures in the deep Earth. We developed a simple method for portraying scalar fields by 3D printing modified globes of surface topography, representing the parameter of interest as additional, exaggerated topography. This is particularly effective for long-wavelength (>500 km) fields. The workflow uses only open source and free-to-use software, and the resulting models print easily and effectively on a cheap (<300 GBP, 400USD) desktop 3D-printer. We have printed 3D representations of different scalar fields, including models of the surface topography of rocky planets, which can be used in outreach events. These objects are powerful to explain the importance of plate tectonics in shaping a planet. The workflow was extended to 3D scalar fields by analogy to Russian nesting dolls, where the audience can remove shallower layers to see how structures change with depth. We applied this to global seismic tomography models resulting in prototypes of “seismic matryoshkas” (see Figure). The tactile nature of these objects ensures that anyone, including the visually impaired, can explore the structures present deep within our planet
The laboratory approach brings a significant simplification compared to the actual conditions in situ. Laboratory experiments represent the only chance to control the physical conditions under which the investigated physical phenomena occur. Acoustic emission (AE) is the process accompanying the brittle fracturing of solid body and simultaneously an indispensable tool for its study. Laboratory experiments under controlled loading conditions make it possible to differentiate the effect of important factors like material structure, stress field, crack occurrence, etc. on fracture initiation and development, and allow simulate the nature in situ. Microearthquakes detected during AE can be analyzed by methods developed in earthquake seismology. We apply the shear-tensile crack (STC) to describe the source mechanism of the AE events with the aim to detect the mode of rock fracturing, in particular to distinguish between a shear slip and tensile crack, the latter both in the phase of its opening and closing. The benefit of discerning between shear and tensile fracturing is an insight into changes of the permeability of the rock massif both in space and time. By contrast to natural seismology, tens of thousands AE events occur in laboratory during the experiment. Expecting to process large volumes of data, an urgent demand was to make the non-linear STC search together with the estimate of the errors involved as fast as possible. To assess the reliability of the STC solution, the confidence regions of source model parameters are constructed. The misfit function is converted into the probability density function which is integrated over a trial volume of low misfit until requested probability content is achieved. For individual microearthquakes, we display confidence regions both for the mechanism orientation and its decomposition. Aiming to process a large bulk of AE data, a method of assessing of these zones needs to be proposed, which describes them by estimates of their extreme size. This allows us select for subsequent interpretation from all solutions only those that are stable and reliable. We have applied this approach to the experimental data obtained from a couple of uniaxial loading tests performed on a Westerly Granite and Liberec Granite specimen using a 14 channel AE monitoring system.
Ultrasound Computer Tomography (USCT) is an emerging technique for breast cancer screening. Ultrasonic waves are propagated through the tissue and recorded by a set of transducers that are surrounding the breast. The experiment collects transmission and reflection data, which are used to obtain quantitative images of acoustic properties of the tissue (see Figure). This information is useful to characterize the breast tissue, and improves the specificity of standard imaging modalities. However, providing a diagnostic tool with high accuracy and clinically affordable time-to-solution (goal: ~15 min/patient) still remains a challenge. The goal of this work is to show that, despite the vast scale differences, experiments in seismology and USCT share many similarities. In both fields, the relative wave speed variations are comparable and the number of propagated wavelengths in the domain has the same order of magnitude. Because the wave equation is scale invariant, the cross-fertilization between both fields will benefit imaging methods on all scales. In this study, we present methods from seismic tomography that we have recently introduced to USCT: 1) We employ a linearized finite-frequency traveltime tomography approach for speed-of-sound reconstruction. Using the cross-correlation traveltime misfit functional, we compute analytically the sensitivity kernels using adjoint techniques. Our method can operate almost in real time while still including finite-frequency effects. It also can retrieve useful 3D information from 2D acquisition systems. 2) Similar to exploration geophysics, both speed-of-sound and reflectivity images are important for the interpretation. Here, we suggest a framework that combines full-waveform inversion for speed-of-sound and reverse time migration for reflectivity. 3) We apply the Sequential Optimal Experimental Design (SOED) method to optimize the position and number of transducers, in terms of accuracy and cost, to image both reflection and transmission information. Using the Bayesian approach, we define the quality of a design as the average of the posterior variances of the parameters. SOED provides cost-benefit curves that quantifies the information gain versus the computational cost. These are useful to control the trade-off between accuracy and practicality.
We perform physics-based simulations of earthquake rupture propagation on geometrically complex strike-slip faults. We consider many different realization of the fault roughness and obtain heterogeneous stress fields by performing dynamic rupture simulation of large earthquakes. We calculate the Coulomb failure function (CFF) for all these realizations so that we can quantify zones of stress increase/shadows surrounding the main fault and compare our results to seismic catalogs. To do this comparison, we use relocated earthquake catalogs from Northern and Southern California. We specify the range of fault roughness parameters based on past observational studies. The Hurst exponent (H) varies in range from 0.5 to 1 and RMS height to wavelength ratio ( RMS deviation of a fault profile from planarity) has values between 10-2 to 10-3. For any realization of fault roughness, the Probability density function (PDF) values relative to the mean CFF change show a wider spread near the fault and this spread squeezes into a narrow band as we move away from fault. For lower value of RMS ratio ( 10-3), we see bigger zones of stress change near the hypocenter and for higher value of RMS ratio ( 10-2), we see alternate zones of stress increase/decrease surrounding the fault to have comparable lengths. We also couple short-term dynamic rupture simulation with long-term tectonic modelling. We do this by giving the stress output from one of the dynamic rupture simulation (of a single realization of fault roughness) to long term tectonic model (LTM) as initial condition and then run LTM over duration of seismic cycle. This short term and long term coupling enables us to understand how heterogeneous stresses due to fault geometry influence the dynamics of strain accumulation in the post-seismic and inter-seismic phase of seismic cycle.
The Superior Province Rifting Earthscope Experiment (SPREE) deployed seismic stations in 2011-2013 throughout Wisconsin, Minnesota, and Ontario. To protect equipment from groundwater damage, SPREE stations were buried at unusually shallow depths, increasing the power of long period noise and facilitating an investigation into the regional effects of atmospheric tides and soil properties (Wolin et al., 2015). Here we utilize the SPREE array to study the effects of solid-earth tides and meteorological conditions, on very long-period seismic noise in the U.S. midcontinent. Continuous seismic data was collected from SPREE and Transportable Array (TA) stations located in Wisconsin and Minnesota (WIMN) between July 2011 and September 2013. This data was “cleaned”, filtered, and averaged to produce a monthly representation of the very-long period signals recorded by the SPREE stations. The signals showed diurnal (24 hr) and semidiurnal (12 hr) periodicities, whose magnitudes and dominance vary seasonally. Using cross correlations, we compare our very-long period observations with theoretical solid-earth tides (Milbert, 2018) as well as meteorological components in the WIMN region. Meteorological data, specifically temperature and pressure, was obtained from the National Oceanic and Atmospheric Administration’s (NOAA) National Center for Environmental Information (NCEI). Solid-earth tides result from the gravitational pull of the moon and sun, and have previously been documented in seismic data (e.g. Pillet et al.,1994; Lambotte et al., 2005). We observe a distinct correlation between theoretical solid-earth tides and very-long period signals in seismic data from SPREE and TA stations in the WIMN region, where one frequency component is correlated while the other appears delayed. In addition, we observe a remarkable seasonal change in SPREE recordings of these signals, but not in TA recordings. We will report our findings from testing the hypothesis that the observed very-long period signals in SPREE data are a combination of both tidal and thermal effects and that these cumulative effects are the result of the unusual burial depth of SPREE stations.
We are engaging citizen scientists in an experiment to test if many human ears can replace the process of a professional seismologist in identifying dynamically triggered seismic events. Ordinarily, this process involves interactive data processing and visualization efforts on a volume of earthquake recordings (seismograms) that exploded during the recent big-data revolution, for example through EarthScope. In this citizen seismology project, we ask citizens to listen to relevant sections of seismograms that are accelerated to audible frequencies. This approach has five advantages: 1) The human ear implicitly performs a time-frequency analysis and is capable of discerning a wide range of different signals, 2) Many human ears listening to the same data provides statistics that rank seismograms in order of their likelihood to contain a recording of a triggered event, which is helpful to researchers’ analysis of this data as well as to 3) the ability of a deep-learning algorithm to model the boolean identifications or bulk statistics of the analyses, 4) the project has the potential to enhance informal learning because of the online platform that hosts the project, Zooniverse, is available to people of all identities and hosts many other citizen science projects, and 5) it helps prepare our team for diverse post-graduation careers as part of IDEAS, an NRT program at Northwestern University. The events we are asking citizens to help identify via listening are small seismic events such as local earthquakes and tectonic tremor, that occur in response to transient stresses from passing seismic surface waves from a large, distant earthquake. While much research progress has been made in understanding how these events are triggered, there is no reliable deterministic recipe for their occurrence. The aim of our project is to enlist the help of citizens to increase the data set of known triggered seismic events and known absences of triggered events in order to help researchers unravel key aspects of that recipe. A better understanding of triggered seismic events is expected to provide important clues towards a fundamental understanding of all seismic activity, including damaging earthquakes.
In the field geophysical instruments, various seismic sensors (an accelerometer, a velocimeter) and auxiliary equipment in the form of recorders (loggers) have been designed, which provide wireless information transfer to the data collecting and processing center. On the basis of these, instrumental earthquake observation systems have been developed to ensure the integrated security of critical facilities. These systems can be transformed to control global seismic monitoring of the territories, dams, also facilities of critical importance such as nuclear power plants, bridges, tunnels and technological pipelines for chemical and oil-processing plants. The report presents the technical characteristics of the above mentioned instruments and the possibilities of their wide application in different countries. Additionally, monitoring of seismicity of Armenian regions is implemented for the earthquake research.
Coseismic and rainfall-triggered landslides are a common hazard in many terrains, and the risk associated with them can be quantified, usually by probabilistic modelling. These events are well-documented as a special case of a cascading hazard chain, and the assessment is commonly done via spatial modelling of susceptibility (suppressing temporal dependence) or tailoring models to specific areas and events. The interaction between Earthquakes and rainfall is not usually implemented in a model, as it is considered coincidental. However, because landslides have multiple triggering factors, there is a need for a statistical model that incorporates both features, in a manner such that the separate and joint effects can be estimated. This helps with understanding the interactions between primary events in the triggering of a single secondary hazard type that is crucial for generally applicable multi-hazard methodologies. The presented work aims at the apportioning of the relative and combined effect on landslide triggering by earthquakes and rainfall using a discrete approximation to a multivariate hierarchical point process. Doing so provides a building block in a general framework with the potential to be extended to other chains of events. A case study on the Italian region of Emilia-Romagna is included, using one of the longest and most complete landslide data sets known. Multiple models for the rainfall-earthquake interaction in landslide triggering are trialed, with the best explanation being additive effects from rainfall intensity, rainfall duration and coseismic triggering.
The NASA InSight mission, landing in late November 2018, promises to revolutionize our understanding of Martian interior structure via analysis of seismic data returned by the SEIS instrument. The extent to which the mission’s potential is realized will depend on the number of detectable seismic events that occur during the period of operation. Here I estimate the rate of detectable events generated by volcano-tectonic activity on Mars based on extrapolation of Hawai’i’s seismic record. I use a catalog of 5603 earthquakes spanning 48 years [IRIS, 2018], with moment magnitudes MW ranging from 3.0 to 6.9, to derive the Gutenberg-Richter (G-R) frequency-magnitude relation for the Island of Hawai’i, expressed as log(N) = a – b MW, where N is the number of earthquakes with magnitudes greater than or equal to MW, and a and b are constants. By this analysis, one earthquake with MW 5.1 or greater can be expected every year at Hawai’i. Under the assumption that the mechanisms of seismicity associated with edifice building are similar at Hawai’i and Olympus Mons (supported by observation of decollement-based volcanic spreading at both), I use the same b for both settings and scale a according to estimates of magmatic volume flux rates dV/dt at both settings. Over the 80 Myr history of the Hawaiian-Emperor volcanic chain, dV/dt ≈ 1.7 x 10-2 km3/yr. An estimate of dV/dt for the Olympus Mons volcano on Mars was derived from paleotopographic analysis of a set of lava flows south of Olympus Mons with discordant topography. Given the mean flow age from crater counts (210 Ma), an estimate of the amount of volcanic material needed to cause deflections of flow orientations of the required magnitudes yields estimates of dV/dt over this timespan ranging from 6.33 x 10-4 to 6.43 x 10-3 km3/s. Taking the mean of these values and scaling a by the ratio of dV/dt values for Mars and Hawai’i yields a rate of at least 1 quake of MW = 4.4 or greater per year (Figure 1). Thus, under several assumptions (including a steady recent magma supply rate for Olympus Mons), we can expect ≈ 2 volcano-tectonically driven quakes of magnitude MW > 4.4 from the vicinity of Olympus Mons during the nominal 2 Earth-year InSight prime mission. This is a conservative lower bound that does not consider contributions from numerous potential volcano-tectonic sources in Tharsis and elsewhere on Mars.
Seismograms recorded at an array of 10 closely spaced seismometers in the western United States are stacked about the ScP phase and display a complex set of phase arrivals in addition to ScP. One dimensional forward modeling with ultra-low velocity zone structures are unable to match the observed waveform. Forward modeling with core-rigidity zone structures, where S-waves can propagate in a thin layer below the core-mantle boundary, provide reasonable waveform fits to the largest amplitude phases. The basic core-rigidity zone model has thickness of 1.1 km, shear wave velocity of 1.8 km/s, and density decrease of 10%, although tradeoffs in model parameters are observed. Waveform fits are insensitive to small changes in P-wave velocity in the core-rigidity zone. To match trailing smaller amplitude phases, S-wave multiples are modeled in the core-rigidity zone. Cross correlation of ScP waveforms from 392 EarthScope seismometer sites surrounding the small array site shows discrete patches of similarity to the above model. These results point to a spatially large but patchy core-rigidity zone beneath Mexico and adds to the small set of previous core-rigidity zone observations globally.
Experimental study of the phase and amplitude observations of sub-ionospheric very low and low frequency signals is performed to analyze the response of the lower ionosphere during the August 21, 2017 total solar eclipse in the United States of America. Three subionospheric wave paths have been investigated. The length of the paths varied from 2200 to 6500 km, signal frequencies were 21.4 kHz, 25.2 kHz and 40.8 kHz. Two paths crossed the region of total eclipse and the third path was in the region of 40-60% of obscuration. None of the signals revealed any noticeable amplitude changes during the eclipse while negative phase anomalies (from -35° to -95°) were detected for all three paths. It was shown that the effective reflection height of the ionosphere in low and middle latitudes has been increased by 3.5-5 km during the eclipse.
Random and small-scale subsurface heterogeneities in velocity and/or density scatter the seismic wavefield when they have scale lengths on the order of the seismic wavelength. Seismic scattering is considered the origin of coda waves. Such inhomogeneities have an important effect on propagating waves, as they generate traveltime and amplitude fluctuations and may be the cause of attenuation or excitation of secondary waves. Understanding the effect of small-scale heterogeneities on the seismic wavefield is important for the characterization of the seismic source (e.g. source parameters of underground nuclear explosions) and to improve our knowledge of the Earth’s structure along the raypath. Several approaches and methods have been suggested to study the scattering of seismic waves and characterise subsurface heterogeneities. Here, we apply a combination of the analysis of the incoherent wavefield component and the coda decay with time to a dataset of over 350 teleseismic events (over 20000 traces) recorded at three seismic arrays (Warramunga, Alice Springs and Pilbara) in Australia. This combination allow us to obtain a series of parameters (correlation length, RMS velocity fluctuations of the heterogeneities and thickness of the scattering layer) that give us a measure of the spatial scale and the magnitude of the heterogeneities present in the lithosphere beneath the arrays. This is the first time such a large dataset is used for a study of these characteristics. Our new results show similar structures and scattering strength for Alice Springs and Warramunga, while revealing a different tectonic signature and stronger scattering in the case of Pilbara, possibly caused by the different thicknesses of crust and lithosphere between these regions and its different tectonic history. These stochastic models of the lithosphere are the first step in the development of a technique analogous to adaptive optics which, in this case, aims at removing the effect of the small-scale, near receiver structure from recorded wavefields, thus enabling us to improve our source characterization and to more clearly image the Earth’s interior.