Geophysical Journal International

SummaryNumerical models of geothermal reservoirs typically depend on hundreds or thousands of unknown parameters, which must be estimated using sparse, noisy data. However, these models capture complex physical processes, which frequently results in long run-times and simulation failures, making the process of estimating the unknown parameters a challenging task. Conventional techniques for parameter estimation and uncertainty quantification, such as Markov chain Monte Carlo (MCMC), can require tens of thousands of simulations to provide accurate results and are therefore challenging to apply in this context. In this paper, we study the ensemble Kalman inversion (EKI) algorithm as an alternative technique for approximate parameter estimation and uncertainty quantification for geothermal reservoir models. EKI possesses several characteristics that make it well-suited to a geothermal setting; it is derivative-free, parallelisable, robust to simulation failures, and in many cases requires far fewer simulations to provide an accurate characterisation of the posterior than conventional uncertainty quantification techniques such as MCMC. We illustrate the use of EKI in a reservoir modelling context using a combination of synthetic and real-world case studies. Through these case studies, we also demonstrate how EKI can be paired with flexible parametrisation techniques capable of accurately representing prior knowledge of the characteristics of a reservoir and adhering to geological constraints, and how the algorithm can be made robust to simulation failures. Our results demonstrate that EKI provides a reliable and efficient means of obtaining accurate parameter estimates for large-scale, two-phase geothermal reservoir models, with appropriate characterisation of uncertainty.

SummaryWe present new numerical tools for geophysical inversion and uncertainty quantification (UQ), with an emphasis on blocky (piecewise-constant) layered models that can reproduce sharp contrasts in geophysical or geological properties. The new tools are inspired by an “old” and very successful inversion tool: regularized, nonlinear inversion. We combine Occam’s inversion with total variation (TV) regularization and a split Bregman method to obtain an inversion algorithm that we call blocky Occam, because it determines the blockiest model that fits the data adequately. To generate a UQ, we use a modified randomize-then-optimize approach (RTO) and call the resulting algorithm RamBO (randomized blocky Occam), because it essentially amounts to running blocky Occam in a randomized parallel for-loop. Blocky Occam and RamBO inherit computational advantages and stability from the combination of Occam’s inversion, split Bregman and RTO, and, therefore, can be expected to be robustly applicable across geophysics.

SummaryThe Guerrero seismic gap in the Mexican subduction zone exhibits a slip behaviour distinct from that of adjacent segments, which typically experience large earthquakes. With the acquisition of offshore seismic data in this region and the discovery of shallow tectonic tremors, the study of slow earthquakes has gradually increased. This study presents the detection of tectonic tremors and low frequency earthquakes (LFEs) in the Guerrero seismic gap using a combination of a modified envelope cross-correlation method and a matched filter applied to Ocean Bottom Seismometer (OBS) data for a continuous two-year observational period. The modified envelope cross-correlation method was used to detect and locate tremors, and the matched filter technique enabled the detection of LFEs. These methods allowed for better constraints on the depths of the detected events, offering new insights into tremors and LFE activity offshore the Guerrero seismic gap. Our results show that the spatial distribution of these phenomena, along with seismicity, residual gravity anomalies, and seafloor topography, suggests that a section of the shallow plate interface within the gap has experienced stable slip. This study builds on previous work by enhancing the detection and location accuracy of these slow earthquakes, contributing to a more comprehensive understanding of subduction dynamics in the region.

SummaryThe XSoDEx (eXperiment of Sodankylä Deep Exploration) project acquired in total 82 km of seismic reflection and refraction data to improve the understanding of the crustal composition and, consequently, the mineral systems of the Sodankylä region in Northern Finland. The Sodankylä region is part of the Central Lapland Greenstone belt, which is famous for its mineral resources. Here, we present the first subsurface images resulting from the seismic reflection and refraction data processing, and provide the first geological interpretation of the data. Our workflow comprised time domain signal processing, migration velocity model building and finally, focussing pre-stack depth migration. The results along the acquired seismic profile lines show a rich inventory of imaged reflectors throughout the upper crust, which in some parts can be correlated clearly with geological features at the surface and also show the complex structure of the lithological units of the Central Lapland Greenstone Belt in the investigation area. Moreover, the presumable top of the Archaean basement can be traced through all lines. The basement is partly bent up to the shallow subsurface. In places, the basement forms a dome-like outcrop. The derived results of the seismic data are in good accordance with earlier interpretations of adjacent seismic investigations. The XSoDEx seismic profiles connect the imaged reflective structures to these surveys, which were acquired over known mineral deposits.

SummaryWe invert 122 147 P, S, and PmP phase arrival-times from 1549 local earthquakes for both isotropic and azimuthally anisotropic lithospheric P-wave velocity structure beneath the region of the Longmenshan Fault zone, China. The use of PmP data significantly improves the spatial resolution of the middle-lower crust tomography. Our results show that widespread low-Vp anomalies exist in the middle and lower crust of the Songpan-Ganzi block and the Chuandian block, which contribute most crustal anisotropy. Moderate and strong earthquakes mainly occurred in the high-Vp and low-Vp transition zone, and obvious low-Vp anomalies appear below the seismogenic zone, indicating that the occurrence of earthquakes is affected by crustal fluids. The upper-crust anisotropy is mainly controlled by the stress field and local faults. The fast Vp directions (FVDs) on the Longmenshan fault zone are NE-SW in the lower crust and uppermost mantle, suggesting that the material flow is blocked by the Sichuan basin, so the flow moves in the NE-SW direction. The FVDs in the Longmenshan fault zone are different from SKS splitting measurements, suggesting that the crust and lithospheric mantle are decoupled there. Our anisotropy results also suggest that the thickening deformation of the upper crust and the middle-lower crustal flow jointly control the uplift and deformation of the Longmenshan mountain.

SummaryThe sensitivity of Rayleigh wave amplitude to Earth structure has applications to seismic tomography, both in cases where amplitude information is used to supplement phase velocity data to improve images of elastic parameters, and to correct amplitudes for local Earth structure in attenuation tomography. We review the theoretical basis of the ray theoretical approximation, in which the wave amplitudes are controlled by a combination of geometrical spreading and local changes in energy density due to Earth structure. We focus mainly on the latter effect, which we term the constant energy flux approximation. We investigate the ray theoretical basis for this approximation, test it against a full waveform simulation that verifies its accuracy, and show how it can be used to compute the sensitivity of amplitude to elastic moduli and density. We investigate how perturbing these parameters in a set of simple Earth models affects Rayleigh wave amplitudes, and demonstrate that a slow velocity heterogeneity can cause either increased or reduced amplitudes, depending upon the depth of the heterogeneity and the observation frequency. Consequently, amplitude sensitivity can be either positive or negative, and its magnitude can vary significantly with frequency. Although an added complication, the very different behavior of phase velocity and amplitudes to changes in Earth structure implies that the two types of data are complementary and suggest the effectiveness of using both in Rayleigh wave tomography.

SummaryGeodetic velocity models, derived from Global Navigation Satellite System (GNSS) velocity solutions, interpolate between sparse GNSS measurements to provide a more comprehensive view of horizontal and vertical intraplate deformation. These models contribute to improved assessment of seismic and volcanic hazards, assist in validating geodynamic models, and enable the integration of diverse datasets for comprehensive Earth science studies. Most interpolation techniques are not adequate to model the velocity distribution, especially for the horizontal velocities where the correlation between the components has to be included. Here, we apply a recent extension of the least-squares collocation interpolation technique to the velocity field solution of the EUREF Permanent GNSS Network Densification (EPND) project (EPND_D2150). The effect of known plate boundaries is accounted for during the interpolation to avoid smoothing across European micro-plates, thereby preserving the high velocity gradients at plate boundaries. The velocity model EuVeM2022 covers Europe and Anatolia, and has a resolution of 0.1○. The model can be applied, amongst other things, to correct models used for ground motion services, support tectonic studies, or identify local deformation along coasts for use in sea-level research.

SummaryThe leaky-mode dispersion extracted from seismograms and noise cross-correlation functions has gained lots of attention in recent years. It has been reported that leaky modes can provide constraints for subsurface structures, especially for P-wave velocities, which may compensate for the limits of traditional surface wave methods. For stable and reliable dispersion-curve inversion, the quantitative analysis of leaky-mode sensitivity is of great importance, which, however, has rarely been studied systematically. Limited by the forward modeling methods, the previous methods for calculating leaky-mode sensitivity are usually hindered by issues of mode skipping, low efficiency, etc. To this end, we propose an effective method that can calculate the leaky-mode sensitivity for various types of models based on the previously proposed forward modeling method named the semi-analytical spectral element method (SASEM). Using the intermediate results of SASEM, we derive analytical expressions for the sensitivity kernels with only matrix operations, which endows the sensitivity calculation procedure with high accuracy and reliability identical to the SASEM. In addition, we suggest a novel modal classification scheme to distinguish different kinds of leaky modes based on the sensitivity features. This scheme facilitates the stable identification of the most attractive guided-P modes from numerous normal and leaky modes, which removes obstacles in the dispersion-curve inversion using guided-P modes to constrain P-wave velocities. Several numerical tests are performed to demonstrate the high accuracy of the sensitivity calculation method and the effectiveness of the modal classification method. To assess the roles of leaky modes in the retrieval of underground structures, we perform comprehensive sensitivity analyses of leaky modes using both crust-scale and near-surface models. Besides the general conclusion that the joint inversion using normal and leaky modes can effectively retrieve P- and S-wave velocities, the feasibility of constraining models with the apparent Σ modes and the split guided-P mode dispersion curves has been demonstrated.

SummaryThe uniaxial compressive strength ${\sigma }_{c}$ of rocks is a key material property in a wide range of applications. Models for ${\sigma }_{c}$ typically either require numerical solutions, restricting their wide utility, or are empirical and therefore confined to a specific case. Here, we study the theoretical pore-emanated crack model and provide an analytical emulator function that matches the 2D and 3D solutions to a high degree of accuracy over all porosities, $\phi $. A key input to both the full solution and to our emulator functions is the pore radius, assumed in the model to be circular or spherical, in a porous rock. In most porous lithologies, including sandstone, the notion of a pore radius is poorly defined since they are built from compacted or lithified grains. And so here we explore statistical methods to find a characteristic pore length scale, ${l}_2$, from an initial particle radius; this method is provided as an easy-to-use supplementary tool. We advocate for the use of our 3D function ${\sigma }_{c} \approx 1.57{K}_{Ic}/( {\phi }^{0.43}\sqrt {\pi }{l}_{2} )$ where ${K}_{Ic}$ is the fracture toughness of the solid matrix. A compilation of ${K}_{Ic}$ values for minerals and rocks allows us to explore the effect of this parameter and to make recommendations for appropriate values in the model. We compare our simple emulator function for ${\sigma }_{c}$ with existing datasets across a wide range of sandstones to demonstrate the utility of this law for applied cases. We find that our function performs particularly well for relatively low porosity sandstones ($\phi \mathbin{\lower.3ex\hbox{$\buildrel<\over {\smash{\scriptstyle \sim}\vphantom{_x}}$}} 0.15$) representative of mature basin systems from a diagenetic point of view; we discuss alternative models that are more appropriate for higher porosity sandstones.

SummaryThe Makran Subduction Zone is a distinctive segment within the Alpine-Himalayan system, where one of the final remnants of the once-expansive Neo-Tethys Ocean is being subducted beneath the Eurasian Plate. Limited seismic data has left several questions unanswered about the structure of the subducting oceanic lithosphere, the transition from the wide and thick Makran accretionary prism to the Zagros Collision Zone, variations in sedimentary cover thickness along and perpendicular to the accretionary prism, and fluctuations in the thickness of sedimentary cover within the fore-arc Jaz Murian Depression (JMD). In this study, we utilize ambient-noise and earthquake surface wave tomography within a period range of 5–50 s to construct a high-resolution 3D shear-wave velocity model down to a depth of 60 km for the Iranian Makran and northern Oman. Using a new dataset from 65 seismic stations located in southeastern Iran and northern Oman, our analysis reveals a sharp velocity contrast within the oceanic lithosphere of the Gulf of Oman, just north of Muscat, with abnormally low-velocity oceanic lithosphere extending westward from this contrast, revealing subduction of a segmented oceanic lithosphere beneath the Makran. Our study finds no lithospheric-scale seismic velocity contrast along the ZMP fault, as usually thought as a transition boundary between the Zagros and Makran. Our velocity model shows that the wide accretionary prism of western Makran consists of two zones: a southern low-velocity zone associated with younger sediments and a northern high-velocity zone corresponding to older sediments. A considerable thinning of the sedimentary cover is observed east of longitude 59° E within the coastal Makran tectono-stratigraphic unit, aligning with the structural trend of the Pan-African Semail Gap Fault observed both onshore and offshore Oman. Additionally, a thick sedimentary basin is located beneath the eastern section of the JMD, with the thickness decreasing towards the west.

AbstractIn this paper we examine the dynamic pressure torque acting on a bumpy core-mantle boundary (CMB) at diurnal timescale in a frame tied to the planet. This torque possibly contributes to the CMB coupling constants determined from nutation observations and could affect the interpretation of these constants in terms of different CMB coupling mechanisms. We revisit the work of Wu and Wahr (1997) who have used seismic estimates for the topography at the CMB and computed the associated pressure torque effect on nutations. These authors showed that some topography wavelengths can lead to amplifications in nutations. For example, they found that the effects on the retrograde annual nutation can be at the milliarcsecond level for a degree-5 spherical harmonics of the topography. While Wu and Wahr (1997) only go up to degree 6 in their development in spherical harmonics and use a numerical technique, we go up to degree 20 and employ an analytical approach to solve the equations and to further study the Earth’s nutations. The approach is similar to the one we used for the effects of the pressure torque on the tidal variations of the length of day (LOD) (a companion paper, Puica et al., 2023). Unlike the numerical approach, this has the advantage of highlighting the mathematical dependencies between the different spherical harmonics involved in the development of the topographic torque and to highlight the frequency dependence of the results and thereby the possible resonances with inertial waves. By doing so, we can isolate and estimate the magnitude of the influence of each topographic coefficient on nutation. We show that only the core flattening may have an important role on nutation and that the other large wavelengths of the topography have a very small contribution, less than that obtained by Wu and Wahr (1997).

SummaryWe conducted comparative measurements of thermal properties of samples from nine cores of the ICDP COSC-1 borehole and four widely used rock references, using a steady-state and a transient divided-bar device, a transient plane source device, a modified Ångström device, as well as two optical thermal conductivity scanners. In addition, a caloric method provided benchmark values for specific heat capacity. A complementary thin-section analysis of the COSC-1 samples allowed us to calculate specific heat capacity according to Kopp’s law and thermal conductivity according to commonly used mixing models. Our results demonstrate agreement between the various test methods within ±10% for about one half of the investigated samples. Furthermore, almost all results for specific heat capacity agree with the predictions of Kopp’s law, though the significance of this correspondence is limited owing to large uncertainties in the experimental and theoretical values. The results for thermal conductivity fall within the most extreme theoretical bounds that account for anisotropy but for an amphibolite. Thermal anisotropy seems to contribute significantly to the deviations between results of the different transient methods that, however, cannot be reconciled by the available theoretical relations for apparent thermal conductivity of transversely isotropic materials. The combination of characteristic investigation volume of the individual methods and sample heterogeneity has to be considered responsible for variability of results, too, an issue whose clarification is calling for dedicated numerical modelling in the future, with the prospect to characterise thermal heterogeneity from observed differences.

SummaryElastodynamic Green’s functions are an essential ingredient in seismology as they form the connection between direct observations of seismic waves and the earthquake source. They are also fundamental to various seismological techniques including physics-based ground motion prediction and kinematic or dynamic source inversions. In regions with established 3D models of the Earth’s elastic structure, such as southern California, 3D Green’s functions can be computed using numerical simulations of seismic wave propagation. However, such simulations are computationally expensive which poses challenges for real-time ground motion prediction and uncertainty quantification in source inversions. In this study, we address these challenges by using a reduced-order model (ROM) approach that enables the rapid evaluation of approximate Green’s functions. The ROM technique developed approximates three-component time-dependent surface velocity wavefields obtained from numerical simulations of seismic wave propagation. We apply our ROM approach to a 50 km × 40 km area in greater Los Angeles accounting for topography, site effects, 3D subsurface velocity structure, and viscoelastic attenuation. The ROM constructed for this region enables rapid computation (≈0.0001 CPU hours) of complete, high-resolution (500 m spacing), 0.5 Hz surface velocity wavefields that are accurate for a shortest wavelength of 1.0 km for a single elementary moment tensor source. Using leave-one-out cross validation, we measure the accuracy of our Green’s functions for the CVM-S velocity model in both the time domain and frequency domain. Averaged across all sources, receivers, and time steps, the error in the rapid seismograms is less than 0.01 cm/s. We demonstrate that the ROM can accurately and rapidly reproduce simulated seismograms for generalized moment tensor sources in our region, as well as kinematic sources by using a finite fault model of the 1987 MW 5.9 Whittier Narrows earthquake as an example. We envision that rapid, accurate Green’s functions from reduced-order modeling for complex 3D seismic wave propagation simulations will be useful for constructing real-time ground motion synthetics and source inversions with high spatial resolution.

SummaryWe aim to improve our comprehension of the seismic process and to identify possible long-term predictability tools of strong earthquakes through the simulation performed by a new-generation simulator code based on a well-elaborated model of the earthquake sources. We applied our previously tested physics-based earthquake simulator to the Nankai megathrust fault system, characterised by a 13 centuries historical record of strong earthquakes. Our results show these significant seismicity patterns characterizing the seismic cycles: the average stress increases almost linearly, while its standard deviation decreases more and more rapidly as the next major earthquake approaches; the co-seismic stress drop and the simultaneous increase of the standard deviation mark the beginning of the new seismic cycle; and the b-value tends to increase some decades before major earthquakes and exhibits correlation with the occurrence rate. Our results encourage further investigations about the application of simulators in support of other methodologies of earthquake forecasting.

SummaryA grain-based stress corrosion model is built from 3DEC-GBM (i.e., a three-dimensional discrete element grain-based model). The model employs the effective stress law and stress corrosion theory to study the time-dependent and time-independent deformation at the mesoscale of the sandstone with varying confining and pore pressures. The simulations adequately explain complex macroscopic time-independent behavior in terms of the mesoscale interaction of grains, and tension cracks were the dominant crack propagation pattern in the simulation for different confining and pore pressures. The traditional creep behavior observed in laboratory brittle creep experiments could also be accurately reproduced by the proposed model. The simulations show that the percentage of tension cracks in rock fractures decreases with increasing confining pressure and pore pressures. Increasing the applied differential stress and reducing the effective pressure can shorten time-to-failure and increase the creep strain rate, respectively. We conclude that the proposed model is an appropriate tool to analyze the deformation behavior of sandstone under coupled hydro-mechanical loading in both the short and long term.

SummaryThis paper introduces a comprehensive framework for modelling both instantaneous and time-dependent elastic softening in anisotropic materials at high pressure and temperature. This framework employs Landau Theory, minimizing the Helmholtz energy by varying isochemical parameters (q) that capture structural changes, atomic ordering, and/or electronic spin states. This allows for internally consistent predictions of volume, unit cell parameters, the elastic tensor, and other thermodynamic properties, while allowing large symmetry-breaking strains. The formulation is validated using the stishovite-to-post-stishovite transition. It is demonstrated that, near this transition, both stishovite and post-stishovite exhibit auxetic behaviour in several directions, with post-stishovite also displaying negative linear compressibility along the long axis of its unit cell (either the a or b axis). The new formulation is implemented in the open-source BurnMan software package.

SummarySeismic and electrical surveys are the most employed geophysical exploration applications for understanding the subsurface earth. Differential effective medium (DEM) models are the models to interpret the seismic and electrical survey data with the greatest success. However, cementation exponent and pore aspect ratio as the indispensable geometric parameters in the electrical and elastic DEM models are independent, making the models not suitable for the joint elastic-electrical modelling, a key requirement for the joint interpretation of seismic and electrical exploration data to better understand the increasingly complex hydrocarbon reservoirs. We show how cementation exponent and pore aspect ratio are correlated in three Berea sandstone samples with changing porosity resulting from varying pore pressure. We find that cementation exponent inverted from the electrical DEM model shows a strong positive linear correlation with pore aspect ratio obtained from the elastic DEM model as an implicit function of porosity induced by increasing pore pressure. We also find that the established linear correlation can enable the DEM models to calculate one physical property (e.g., elastic or electrical) from the geometric parameter describing the other property (e.g., electrical or elastic). The results reveal how the elastic and electrical geometric parameters are linked, and provide a consistent microstructure that enables the existing elastic and electrical DEM models to be suitable for the joint elastic-electrical modelling of rocks undergoing varying pore pressure.

SummaryThe Tibetan Plateau, a critical region influencing both local and global atmospheric circulation, climate dynamics, hydrology, and terrestrial ecosystems, is undergoing climate-driven changes, including glacial retreat, permafrost thaw, and groundwater changes. Despite its importance, implementing continuous and systematic observations have been challenging due to the area’s high altitude and extreme climate conditions. In this context, seismic interferometry emerges as a cost-effective method for the continuous monitoring of subsurface structural changes driven by environmental factors and internal geophysical processes. We investigate subsurface evolution using four years of seismic data from nine stations on the northeastern Tibetan Plateau, by applying coda wave interferometry across multiple frequency bands. Our findings highlight seismic velocity changes within the frequency bands 5–10 Hz, 0.77–1.54 Hz, and 0.25–0.51 Hz, revealing depth-dependent seasonal and long-term changes. Near-surface and deeper strata exhibit similar seasonal patterns, with velocities increasing in winter and decreasing in summer driven by changes in hydrological processes, while intermediate ice-water phase strata show contrasting behavior due to thermal elastic strain. Long-term trends suggest that the upper subsurface layer is affected by melting water and precipitation originating from Kunlun Mountains, whereas deeper layer reflect groundwater level variations influenced by climate change and human activities. This study provides insights into the environmental evolution of the Tibetan Plateau and its impact on managing local groundwater resources.

SummarySeismic tomography is a principal method for studying mantle structure, but imaging of Earth’s wavespeed anomalies is conditioned by seismic wave sampling. Global models use misfit criteria that may strive for balance between portions of the data set but can leave important regional domains underserved. We evaluate two full-waveform global tomography wavespeed models, GLAD-M25 and SEMUCB-WM1, in the mantle below the Pacific Ocean. The region of the South Pacific Superswell contains multiple hotspots which may be fed by plumes anchored in the Large Low Shear-Velocity Province at the base of the mantle. The uneven distribution of seismic receivers worldwide leaves several candidate plumes beneath various hotspots poorly resolved. We assess the regional quality of GLAD-M25 relative to its global performance using a partition of the seismic waveform data used in its construction. We evaluate synthetic waveforms computed using the spectral-element method to determine how well they fit the data according to a variety of criteria measured across multiple seismic phases and frequency bands. The distributions of travel-time anomalies that remain in GLAD-M25 are wider for trans-Pacific paths than globally, suggesting comparatively insufficiently resolved seismic velocity structure in the region of interest. Hence, Pacific-centered regional inversions, based on (augmented) subsets of the global data set have the potential to enhance the resolution of velocity structure. We compare GLAD-M25 and SEMUCB-WM1 by cross-validation with a new, independent, data set. Our results reveal that short- and long-wavelength structure is captured differently by the two models. Our findings lead us to recommend focusing future model iteration on and around the Pacific Superswell and adding data that sample new corridors, especially using ocean sensors, to better constrain seismic velocity structure in this area of significant geodynamic complexity.

SummaryWe present a new 3D crustal P-wave velocity (VP) model for the greater Alpine region (GAR). We use and merge three different high-quality datasets for local earthquake tomography covering 24 years, starting from January 1st, 1996, up to December 31st, 2019. We processed and repicked the waveforms from the events reported by the European-Mediterranean Seismological Centre with M > 3.0 inside the greater Alpine region for the period between May 2007 and December 2015 using a recently developed automated arrival time-picking procedure (ADAPT framework). This allows bridging the data gap between previously published (pre-2007) datasets and the recently published AlpArray research seismicity catalogue and thus provides a high-quality, highly consistent set of P-wave arrival times covering 24 years. With this data set we derived a new minimum 1D VP model and associated station delays covering the entire GAR. Subsequently, we performed a series of local-earthquake-tomography (LET) inversions obtaining a 3D VP model with a horizontal node spacing of 20×20 km and between 7 to 15 km variable vertical spacing in the well-resolved area of investigation, thus improving the spatial and uniformly high-resolution coverage compared to previous LET studies in the area. For well-known major crustal structures, such as, e.g. the geophysical Ivrea body, deep foreland basins and main orogenic crustal roots, our tomographic results correlate well with features documented by various previous seismic studies in the region. This correlation increases our confidence in the model's accuracy throughout the well-resolved area. Additionally, our model reveals previously poorly known, or unknown crustal features and it documents details in the Moho topography throughout the region. Eventually, we present a LET-Moho map (VP isoline of 7.25 km/s) for the GAR with spatially nearly uniform resolution and document its comparison with previously published Moho maps. The new regional 3D VP crustal model also correlates well with a previously published VS crustal model obtained by ambient noise tomography. These comparisons document the new LET results of combined 3D VP crustal velocities and Moho topography being intrinsically consistent and reliable within the region of high resolution. Hence, in addition to further improving our understanding of crustal structure geometries in the GAR, our results also provide pivotal information for a future reference seismic 3D crustal model of the region.