Imaging Of Subsurface Structures Research Articles

What Extra Information Can Be Provided by Multi-Component Seismic Data: A Case Study of 2D3C Prospecting of a Copper–Molybdenum Mine in Inner Mongolia, China

With the decrease in shallow mineral reserves, deep mineral resources have become the focus of exploration. Seismic exploration, renowned for its deep penetration and high spatial resolution and precision, stands as a primary technique in geophysical exploration. In comparison to traditional P-wave seismic exploration, multi-component seismic techniques offer the advantage of simultaneously acquiring P-wave and S-wave data, overcoming the limitations of single P-wave impedance in predicting lithology and enabling high-precision imaging of subsurface structures. Constrained by field survey costs, the reflection seismic illumination is lower and results in a poor signal-to-noise ratio of multi-component seismic data in metallic ore exploration, which poses great challenges in imaging converted S-waves. Based on the seismic and geological characteristics of metallic ores, this study conducts imaging research on metallic ore models through synthetic data and field multi-component seismic data from a copper–molybdenum mine in Inner Mongolia, China. The emphasis is given to PS-wave pre-stack time migration based on precisely sorting the commonly converted point so as to explore the feasibility and technical advantages of multi-component seismic exploration in metal mines. Synthetic data and field data testing demonstrate that PS-wave imaging contains more abundant structural and lithological information compared to PP-waves, indicating promising prospects for the application of multi-component seismic data in metallic ore exploration.

Moho Imaging with Fiber Borehole Strainmeters Based on Ambient Noise Autocorrelation.

Moho tomography is important for studying the deep Earth structure and geodynamics, and fiber borehole strainmeters are broadband, low-noise, and attractive tools for seismic observation. Recently, many studies have shown that fiber optic seismic sensors can be used for subsurface structure imaging based on ambient noise cross-correlation, similar to conventional geophones. However, this array-dependent cross-correlation method is not suitable for fiber borehole strainmeters. Here, we developed a Moho imaging scheme for the characteristics of fiber borehole strainmeters based on ambient noise autocorrelation. S-wave reflection signals were extracted from the ambient noise through a series of processing steps, including phase autocorrelation (PAC), phase-weighted stacking (PWS), etc. Subsequently, the time-to-depth conversion crustal thickness beneath the station was calculated. We applied our scheme to continuous four-component recordings from four fiber borehole strainmeters in Lu'an, Anhui Province, China. The obtained Moho depth was consistent with the previous research results. Our work shows that this method is suitable for Moho imaging with fiber borehole strainmeters without relying on the number of stations.

Noise source localization using deep learning

Summary Ambient noise source localization is of great significance for estimating seismic noise source distribution, understanding source mechanisms and imaging subsurface structures. The commonly used methods for source localization, such as the matched field processing and the full-waveform inversion, are time-consuming and not applicable for time-lapse monitoring of the noise source distribution. We propose an efficient alternative of using deep learning for noise source localization. In the neural network, the input data are noise cross-correlation functions and the output are matrices containing the information of noise source distribution. It is assumed that the subsurface structure is a horizontally layered earth model and the model parameters are known. A wavefield superposition method is employed to efficiently simulate ambient noise data with quantities of local noise sources labelled as training datasets. We use a weighted binary cross-entropy loss function to address the prediction inaccuracy caused by a sparse label matrix during training. The proposed deep learning framework is validated by synthetic tests and two field data examples. The successful applications to locate an anthropogenic noise source and a carbon dioxide (CO2) degassing area demonstrate the accuracy and efficiency of the proposed deep learning method for noise source localization, which has great potential for monitoring the changes of the noise source distribution in a survey area.

Illuminating Urban Near-Surface with Distributed Acoustic Sensing Multimodal Noise Surface-Wave Imaging

Abstract Urban subsurface exploration requires high spatial and temporal resolution, cost-effective operation, and minimal interference with urban activities. Distributed acoustic sensing (DAS)—an innovative seismic observation tool—emerges as a promising solution for urban surveys. In this study, we repurposed a 7.9 km telecommunication cable traversing Hefei into a seismic observation array with 3850 channels spaced at 2 m intervals. Noise cross-correlation functions (NCFs) were constructed from recordings by iDAS2 and ZD-DAS interrogators along the entire cable. Spatial variation in the NCFs was observed and attributed to different traffic conditions. Employing the recently developed modified frequency–Bessel transform method to NCFs from the 2 km southern subsection of the optic cable, we extracted broadband, high-resolution multimodal dispersion curves. The inverted near-surface structure beneath the cable unveiled a sediment thinning trend from the center to the periphery of the Hefei basin, consistent with borehole inspections. The three-station interferometry (C3) method and beamforming with the Bessel kernel function are applied to mitigate challenges arising from the weak coupling between the cable and the Earth, as well as persistent localized noise sources. These techniques facilitated the acquisition of broadband surface waves. Distinct secondary scatters are observed in NCFs near channels 2090 and 2287, accompanied by a substantial velocity contrast of 30%–40%, suggesting the existence of a blind fault. The study reaffirms the significant potential of DAS arrays for high-resolution imaging of subsurface structures in challenging urban environments, emphasizing the importance of advanced processing techniques to enhance imaging accuracy and robustness.

Fine Shear-Wave Velocity Structures of Subsurface beneath the Guangdong–Hong Kong–Macao Greater Bay Area with Dense Seismic Array and SPAC Method

Abstract We apply the spatial autocorrelation (SPAC) method to construct the 3D subsurface shear-wave velocity structure model using the short-period dense seismic array (containing 725 nodal geophones) located at the Guangdong–Hong Kong–Macao Greater Bay area (GBA). We first divided the dense array into numerous subarrays, with each subarray consisting of nine nodal geophones, and obtained 562 subarrays that can provide 1D VS profiles of the same quantity. Then, the SPAC method and genetic algorithm are utilized to extract the dispersion curve of the Rayleigh wave from the raw microtremor data and invert VS structure, respectively. Finally, a 3D VS structure model from the surface to 3.3 km depth is derived by combining all 1D VS structures. Relatively low-velocity anomalies above 700 m are considered unconsolidated shallow sediments as well as relatively high-velocity anomalies beneath 1100 m are attributed to consolidated granite bedrock. Meanwhile, low-velocity anomalies that are identified through the vertical VS profile at a depth of about 900–3000 m can be contributed to the fractured zone, and striped low-velocity anomalies in the horizontal VS maps reveal the location of the deeply buried faults in the study area. The results also mean that the SPAC method combined with the records of short-period dense seismic array can be effectively applied to image subsurface structures in high-populated urban area. The development of this noise-resistance and environment-friendly geophysical technique provides a reliable and effective way to explore the complicated subsurface geological structures, which is of great significance to urban engineering construction and earthquake disaster reduction work in densely populated urban agglomerations.

Improving full-waveform inversion based on sparse regularization for geophysical data

Abstract Full-waveform inversion (FWI) is an advanced geophysical inversion technique. FWI provides images of subsurface structures with higher resolution in fields such as oil exploration and geology. The conventional algorithm minimizes the misfit error by calculating the least squares of the wavefield solutions between observed data and simulated data, followed by gradient direction and model update increment. Since the gradient is calculated by forward and backward wavefields, the high-accuracy model update relies on accurate forward and backward wavefield modelling. However, the quality of wavefield solutions obtained in practical situations could be poor and does not meet the requirements of high-resolution FWI. Specifically, the low-frequency wavefield is easily affected by noise and downsampling, which influences data quality, whereas the high-frequency wavefield is susceptible to spatial aliasing effects that produce imaging artefacts. Therefore, we propose using an algorithm called sparse relaxation regularized regression to optimize the wavefield solution in frequency-domain FWI, which is the forward and backward wavefield obtained from the Helmholtz equation, thus improving FWI's accuracy. The sparse relaxation regularized regression algorithm combines sparsity and regularization, allowing the broadband FWI to reduce the effects of noise and outliers, which can provide data supplementation in the low-frequency band and anti-aliasing in the high-frequency band. Our numerical examples demonstrate the wavefield optimization effect of the sparse relaxation regularized regression-based algorithm in various cases. The improved algorithm's accuracy and stability are verified compared to the Tikhonov regularization algorithm.

Broadband sub-wavelength terahertz subsurface imaging using a solid-immersion lens

We demonstrate a broadband sub-wavelength terahertz (THz) imaging system by combining a solid-immersion lens (SIL) with a THz time-domain spectrometer (TDS). The imaging resolution of the system is enhanced in a wide spectral range (from 0.3 to 2.0 THz) with the aid of the SIL. The beam size is within 120–280 μm (0.28–0.8 λ) in the measured frequency range. Experimentally achieved imaging resolution enhancement results are consistent with the numerical simulation results. Unlike the previous reports on imaging systems using SIL, the beam size reduction by the SIL is more than a factor of nSIL (3.42) for all frequencies near the focal spot. An analytical calculation is employed to obtain insights into the beam size reduction mechanism and estimate the beam size in a given experimental condition. THz image of the WR1.9 waveguide is measured to highlight the broadband sub-wavelength imaging capability. The increased reflection from the waveguide near the cutoff frequency for TE10 and TE01 modes, inaccessible with other THz-TDS-based or continuous wave (CW) source-based THz imaging systems, is experimentally observed. In addition, we show the multispectral sub-wavelength images of subsurface structures by selecting the appropriate measurement time window of the temporal THz signal. THz images of a credit card show the sub-wavelength details of the card, including an IC chip, invisible due to the plastic cover. The experimental results highlight the potential applications of the imaging system in non-contact, non-destructive testing of the sub-wavelength scale structures in various fields.

Frequency-dependent Q simulation and viscoacoustic reverse time migration based on the fractional Zener model

Seismic attenuation is a basic physical property of the earth, which significantly affects the characteristics of seismic wavefields. Accurately simulating wave propagation in the earth is essential to image subsurface structures. Some prevailing methods (e.g., the standard linear solid and fractional Laplacian equation) to describe seismic wave propagation in attenuating media are mainly based on the constant- Q model (CQM), which is valid at room temperature and pressure. However, laboratory measurements suggest that the quality factor Q is a function of frequencies in some regions. To simulate the frequency-dependent Q effect, we derive a viscoacoustic wave equation from the stress-strain relationship of the fractional Zener model (FZM) with variable fractional orders. During the implementation, we separate the real and imaginary parts of the modulus and introduce a low-rank decomposition method to solve the FZM equation. Because the amplitude dissipation and phase dispersion are decoupled, we establish a compensated reverse time migration ( Q-RTM) algorithm to mitigate adverse effects caused by seismic attenuation and improve the quality of seismic migration in frequency-dependent attenuating media. A two-layer model and the BP gas chimney model are used to perform Q-RTM tests. A low-pass filter with a Tukey window function is applied to suppress numerical instability during the compensation. Numerical results demonstrate that our FZM Q-RTM approach can produce high-resolution images with corrected reflector positions and amplitudes. Because the CQM equation ignores the frequency dependence of Q, it may lead to overcompensation in Q-RTM.

Thermal and imaging subsurface structure of Farasan and Dahlak Islands, southern Red Sea, derived from the interpretation of marine geophysical data

Thermal and imaging subsurface structure of Farasan and Dahlak Islands, southern Red Sea, derived from the interpretation of marine geophysical data

The utilization of ghost reflections retrieved by seismic interferometry for layer‐specific characterization of the shallow subsurface

AbstractSeismic interferometry (SI) retrieves the Green function between two receiver locations using their recordings from a boundary of sources. When using sources and receivers only at the surface, the virtual‐source gathers retrieved by SI contain pseudo‐physical reflections as well as ghost (non‐physical) reflections. These ghost reflections are the results of the cross‐correlation or auto‐correlation (AC) of primary reflections from two different depth levels, and they contain information about the seismic properties of specific layers in the subsurface. We investigated the application of ghost reflections for layer‐specific characterization of the shallow subsurface using SI by AC. First, we showed the technique's potential using synthetic data for a subsurface model with a lateral change in velocity, a gradient in depth for velocity, a thickness change and a velocity change of the target layer. Then, we applied the technique to shallow subsurface field data. We also focused on improving the retrieval of ghost reflections by removing the free‐surface multiples and muting undesired events in active‐source gathers before applying SI. Our results demonstrate that the ghost reflections can be used advantageously to characterize the layer that causes them to appear in the results of SI. Consequently, they can also provide valuable information for imaging and monitoring shallow subsurface structures.

Receiver Function Adjoint Tomography for Three‐Dimensional High‐Resolution Seismic Array Imaging: Methodology and Applications in Southeastern Tibet

AbstractA new technique for P‐wave receiver function (PRF) inversion, within the framework of wave equation‐based adjoint tomography and referred to as receiver function adjoint tomography (RFAT), has been developed to obtain models of Vp, Vs, and density. This innovative technique fits the synthetic PRFs with observed PRFs and can better image the lateral variations of Vs from the crust to the uppermost mantle than traditional 1‐D PRF inversion. We utilized RFAT to perform high‐resolution imaging beneath a dense seismic array in Southeastern Tibet, revealing low‐velocity zones extending from the uppermost mantle to the crust, as well as an eastward dipping Moho under the Red River Fault (RRF). Our inversion results provide direct evidence for the existence of a distinct asthenospheric upwelling channel beneath the RRF, and further highlight the effectiveness of RFAT for accurately imaging subsurface structures.

Seismic Full Waveform Inversion Accelerated by Overlapping Data Input and Computation

Seismic full waveform inversion (FWI) is a powerful technology to obtain high-precision and high-resolution images of subsurface structures. However, FWI is a data-intensive algorithm that needs to read extensive seismic data from disks, which significantly affects its performance. We proposed a portable parallel framework to improve FWI by overlapping data input and computation (ODIC). The framework is based on POSIX threads (Pthreads), which is a standard thread API library and can create a parent thread and a child thread in the FWI process. The former is used to perform computation and the latter to read data from disks, both running simultaneously. This framework has two attractive features. First, it is broadly applicable; it can run on almost any computer from a laptop to a supercomputer. Second, it is easy to implement; it can be readily applied to existing FWI programs. A 3D FWI example shows that the framework speeds up FWI considerably.

Retrieval of surface waves from high-speed-train-induced vibrations using seismic interferometry

Train traffic has recently been realized as a potential seismic source for imaging and monitoring shallow subsurface structures. High-speed trains usually move on viaducts and generate seismic wavefields through bridge piers. This kind of source is not randomly distributed in space and time and thus indicates the characteristic of self-correlation. We find that correlated sources of one high-speed train running over piers make time-delayed spurious signals in the crosscorrelation functions. The direct-wave interferometry with correlated sources is proposed and demonstrates that stacking the interferometric responses from high-speed trains with different speeds may reduce the spurious imprint of source self-correlations, but it cannot remediate the problem. The influence of the source wavelet can be eliminated by using cross coherence. Stationary-phase segment selection and stacking can reduce the bias generated by trains moving out of the stationary-phase zone of station pairs. The validity is demonstrated by using synthetic data with inline and out-of-line linear arrays. We apply the proposed scheme to field data with nine high-speed train events. Multimodal Rayleigh waves were retrieved and shallow shear-wave velocities were inverted from multimodal dispersion data in a range of 4–25 Hz. The results bring a novel approach for the imaging and monitoring of the shallow subsurface using high-speed-train-induced vibrations.

Rift Basins and Intraplate Earthquakes: New High‐Resolution Aeromagnetic Data Provide Insights Into Buried Structures of the Charleston, South Carolina Seismic Zone

AbstractThe delineation of faults that pose seismic risk in intraplate seismic zones and the mapping of features associated with failed rift basins can help our understanding of links between the two. We use new high‐resolution aeromagnetic data, previous borehole sample information, and reprocessed seismic reflection profiles to image subsurface structures and evaluate recent fault activity within the Charleston seismic zone, the associated Mesozoic South Georgia rift basin, and surrounds. The new aeromagnetic data provide an unprecedented view of buried basement structures. NE‐ and NW‐trending lineaments of various lengths throughout the survey area are interpreted as Paleozoic orogenic structures and Mesozoic dikes, respectively. Within the rift basin, 15‐ to 20‐km long ESE‐trending lineaments are associated with faults in pre‐Cretaceous strata of the reflection data and are interpreted as Mesozoic rift structures. Various intersections and terminations of interpreted faults suggest rift‐related reactivation of Paleozoic faults and corresponding inheritance for Mesozoic structures. The reflection data show that several Paleozoic and Mesozoic faults are associated with deformation in Cretaceous and younger sediments, suggesting reactivation in the more recent passive margin setting. Two of these faults, one NE‐striking and one ESE‐striking, are coincident with surficial landforms, suggesting Quaternary slip; the ESE‐striking fault is also well‐aligned with a plan‐view offset in modern seismicity. A favorable orientation for reverse motion on ESE‐striking Mesozoic faults, a possible sub‐basin, and potentially weakened lithosphere are failed rift basin features that may influence intraplate seismicity within the Charleston seismic zone.

Characteristic reflector-based wave equation reflection traveltime inversion

Reconstructing the deep background velocity model is crucial for imaging subsurface structures. The reflection traveltime inversion (RTI) is capable of recovering the kinematic information in the velocity model under a limited seismic observation aperture. The estimation of the reflection traveltime residuals between the synthetic and observed seismic waves is the key to the RTI. However, it is difficult to accurately obtain the rapidly changing reflection traveltime shifts in the prestack data domain. Estimating the traveltime difference is especially difficult where the discrepancy between the noisy observed data and demigrated synthetic data includes a waveform distortion, inconsistent reflected events, and so on. To overcome this problem, we have developed a new reflection traveltime residual estimation method. In this method, after identifying a characteristic reflector in the prestack depth migration section with an initial velocity, the demigrated and observed seismic reflections in the data domain corresponding to the characteristic reflector share the same traveltime at the zero offset. Thus, their traveltime residual can be automatically tracked by their crosscorrelation with a modified dynamic programming algorithm based on a zero-offset control point, which can avoid traveltime picking or windowing of relevant reflection arrivals in the prestack seismic data domain. Then, the estimated reflection traveltime residual is applied to the RTI. To improve the inversion accuracy, we have further developed a layer-stripping velocity model-building workflow. Synthetic and field data tests demonstrate that the proposed method can obtain a robust reflection traveltime residual and further achieve a correct low-to-intermediate wavenumber velocity model update.

Weighted stacking of ambient noise correlations based on the coherency of 3C data

Seismic interferometry (SI) is widely used in imaging subsurface structures at different scales by extracting empirical Green's functions (EGFs) of surface wave propagation from ambient noise. The EGFs may be misestimated, however, because the anisotropic distribution of noise sources violates the SI's basic assumption of energy equipartitioning. A common solution is to stack all short-duration noise cross-correlation functions (NCFs) over a long period of time to attenuate the effect of anisotropic noise sources. The effectiveness of the common method is limited by the time of observation, and it may even be counterproductive. The influence of anisotropic noise sources may be reinforced if these sources possess strong energy. To obtain reliable EGFs, we develop a new weighted stack method named the components-coherency-weighted stack (CCWS). We calculate the coherency of different components of noise recorded by a single station to scan the azimuth of the main source in one segment. We select the segments according to the azimuth for NCF computation and then stack the NCFs with a weight relevant to the coherency. We first test the method on simulated data containing strong azimuthal noise sources and show that our method yields high quality EGFs. The method is then used for two real data examples at a regional scale of hundreds of kilometers and an engineering scale of tens of meters. Two examples confirm the ability of our approach to improve the EGFs and its applicability at different scales.

Imaging the source region of recent megathrust earthquakes along the Chile subduction zone: A summary of results from recent experiments

In addition to the largest documented earthquake, the Mw 9.5 Valdivia earthquake of 1960, the Chilean subduction zone has experienced three great (>Mw 8) earthquakes since 2010 for which detailed models of deformation prior to, during and after the earthquake are available. The 1960 and 2010 Mw 8.8 Maule earthquake affected the south-central Chilean margin, where the trench is sediment filled, whereas the 2014 Mw 8.2 Iquique and 2015 Mw 8.3 Illapel earthquakes occurred farther north where there is little sediment in the trench. We summarize recent results from three marine seismological projects designed to elucidate the impact that geologic structure may have had on slip during these events. Because analysis of these datasets is ongoing, this report should be considered as an interim summary rather than a comprehensive synthesis. None-the-less, some patterns can be discerned. That we can resolve 5 m of margin uplift by differencing swath bathymetric data highlights the value of obtaining ground truth marine geophysical data from potential seismic gaps. Controlled source seismic imaging indicates that subduction of nearly all incoming sediment in the sediment-filled section of the trench is a widespread characteristic of the subduction zone segment that ruptured in the two largest events and is likely a generalizable characteristic that can be used to anticipate whether a subduction zone is likely to rupture in great tsunamigenic earthquakes. While this is not a new idea, our data provide additional insights into plate boundary evolution at depth in this scenario. We also find that subduction of anomalously rough basement topography ∼2 million years ago may have led to indentation, large scale internal deformation, and erosion of the overlying wedge north of the 2010 rupture zone and that the process of wedge recovery may have resulted in an abrupt offset along strike of the plate boundary that served as a barrier to slip propagation in 2010. In general, structures generating wedge disruption may no longer be detectable following their subduction deep beneath the forearc wedge, but their passage can be inferred from residual structures mapped within the wedge, and their impact on slip during large earthquakes may be long lasting. We anticipate that further analysis of the data from these experiments, and from future experiments that build on these results, will further elucidate the relationship between geologic structure and slip such that images of subsurface structure become an additional tool in the toolbox of observations used to evaluate future hazard.

Implicit Seismic Full Waveform Inversion With Deep Neural Representation

AbstractFull waveform inversion (FWI) is arguably the current state‐of‐the‐art amongst methodologies for imaging subsurface structures and physical parameters with seismic data; however, important challenges are faced in its implementation and use. Keys amongst these are (a) building a suitable initial model, from which a local minimum is unlikely to be reached, and (b) availability of tools for evaluation of uncertainty. An algorithm we refer to as implicit full waveform inversion (IFWI), designed using continuously and implicitly defined deep neural representations, appears in principle to address both of these issues. We observe in IFWI, with its random initialization and deep learning optimization, improved convergence relative to standard FWI model initialization and optimization. Models close to the global minimum, capturing relatively high‐resolution subsurface structures, are obtained. In addition, uncertainty analysis, though not solved in IFWI, is meaningfully addressed by approximating Bayesian inference with the addition of dropout neurons. Numerical experimentation with a range of 2D geological models is suggestive that IFWI exhibits a strong capacity for generalization, and is likely well‐suited for multi‐scale joint geophysical inversion.

Deformation of the Qinling belt revealed byP-wave velocity and azimuthal anisotropy tomography

SUMMARYThe Qinling belt is a transitional zone lying among three units: the North China block (NCB), the South China block (SCB) and the northeastern Tibetan Plateau (NETP). Owing to the interaction of these units, complex deformation has occurred in the Qinling belt. Although many studies have been conducted to understand the deformation mechanism in the Qinling belt, some key issues are still under debate, such as whether middle-lower crustal flow exists beneath the western Qinling belt (WQB). High-resolution images of subsurface structures are essential to shed light on the deformation mechanism. In this paper, high-resolution images of the velocity structure and azimuthal anisotropy beneath the Qinling belt are obtained by using an eikonal equation-based traveltime tomography method. Our seismic tomography inverts 38 719 high-quality P-wave first arrivals from 1697 regional earthquakes recorded by 387 broad-band seismic stations. In the WQB, our tomography results show low-velocity anomalies but relatively weak anisotropy in the middle-lower crust. These features suggest that middle-lower crustal flow may not exist in this area. In the central Qinling belt (CQB), we find low-velocity anomalies in the middle-lower crust; however, the fast velocity directions no longer trend E–W but vary from NNE–SSW to N–S. These characteristics can be ascribed to the convergence and collision between the NCB and the SCB. In addition, we find strong low-velocity anomalies in the uppermost mantle beneath the CQB, which may indicate delamination of the lower crust. In the southern Qinling belt, we observe significant high-velocity anomalies in the upper crust beneath the Hannan–Micang and Shennong–Huangling domes. These high-velocity anomalies indicate a mechanically strong upper crust, which is responsible for the arc-shaped deformation process of the Dabashan fold. Based on the P-wave velocity and azimuthal anisotropic structures revealed by the inversion of high-quality seismic data, the deformation of the Qinling belt is affected mainly by the convergence between the NCB and the SCB rather than by the middle-lower crustal flow from the Tibetan Plateau.

Geophysical Imaging of the Shallow Geyser and Hydrothermal Reservoir Structures of Spouter Geyser, Yellowstone National Park: Geyser Dynamics I

AbstractYellowstone National Park (YNP) is home to roughly 500 geysers, making it the most concentrated geyser field in the world. Recent studies exploring the mechanics of geyser eruptions utilizing laboratory models were limited by a lack of knowledge of the subsurface geometry of the geyser features. This study presents results from active hydrogeophysical surveys, including ground penetrating radar, nuclear magnetic resonance (NMR), seismic refraction, electrical resistivity tomography, and transient electromagnetics to image subsurface geyser structure and constrain the geophysical response of the reservoir structure of Spouter Geyser, Yellowstone National Park. Previous geophysical studies on similar geyser systems in Yellowstone characterize the hydrothermal reservoir that supplies the geyser with hydrothermal fluids and vapors as high porosity, low density structures. Whereas our imaging highlights the hydrothermal conduits and reservoir as high resistivity and high velocity structures, interpreted as silica precipitate decreasing the porosity and increasing the bulk modulus (i.e., from unconsolidated media to semi/fully consolidated). Ground penetrating radar identifies sinter thickness around the geyser, NMR provides porosity measurements of the geyser reservoir, and transient electromagnetics maps the depth to the Biscuit Basin Rhyolite bedrock throughout the field site. Electrical resistivity tomography and seismic refraction image the shallow geyser conduit structures to ∼15 m depth and establish the extent of the geyser hydrothermal reservoir structure. Furthermore, a spatial correlation of the electrical resistivity and seismic velocity results establishes the interpreted hydrothermal geyser reservoir as a high resistivity, high velocity structure located to the northeast of Spouter Geyser at depths greater than 15 m.