Permeability quantifies the efficiency of rocks in transporting fluids; its accurate estimation is crucial for formation evaluation and fluid production forecasts. In shaly sandstones, permeability depends on pore structure and shale or clay content, which may be laminated and/or dispersed. Laminated shale introduces flow anisotropy, whereas dispersed clays constrict pore space, thus reducing porosity and increasing fluid friction. Existing models for permeability estimation typically focus on either laminated shale or dispersed clay, but they are limited in their reliability when both are present simultaneously. We introduce a new permeability estimation model designed specifically for shaly sandstones, accounting for the simultaneous presence of thin shale laminations and dispersed clay. It incorporates the volumetric concentrations of shale lamination and grain-coating clay and accounts for their respective permeability endpoints in conjunction with clean-sandstone permeability to improve the estimation of in situ rock permeability from well logs. Furthermore, we develop an advanced approach for implementing nuclear magnetic resonance (NMR) permeability models, particularly in the presence of thinly laminated shale and immiscible fluids. By integrating our NMR interpretation method with the new permeability model, we significantly improve the accuracy of the estimated formation permeability, hence obtaining a better match between rock core measurements and estimated permeability values. We validate and benchmark the permeability calculated with our model and the newly developed NMR interpretation method against core measurements from a published data set and two field cases of heterogeneous sandstones. In both field cases, we observe a 50% reduction in the deviation of permeability results from core measurements, indicating a significant improvement in the estimation of flow properties despite the associated rock complexity.

We develop a novel asymptotic local finite-difference (ALFD) method for solving the fractional Laplacian wave equation in the time domain. This numerical method uses the trigonometric function-based generating function to derive the quasidifference operator, facilitating the discretization of the fractional Laplacian without requiring fast Fourier transforms during simulations. Theoretical analysis and numerical experiments demonstrate that the ALFD method is accurate and efficient. Stability analysis using the eigenvalue method reveals that our method enjoys a relaxed stability condition (the minimum Courant number of 0.659492), providing greater flexibility in the selection of sampling intervals. Numerical dispersion anisotropy is significantly mitigated for various very small Q values ( Q < 10). The numerical solutions obtained via the ALFD method march closely with the analytical solutions derived from Green’s function. The computation speed of ALFD is significantly higher than the analytical method and the original global finite-difference (GFD) method. For a 2D model with 400 × 400 grid points, the computational time required by the ALFD method is approximately 9.96%, 5.40%, 4.94%, and 4.60% of that required by the analytical method, and only 2.54%, 1.37%, 1.25%, and 1.15% of that required by the GFD method, respectively. Numerical examples illustrate the effectiveness of the ALFD method in simulating heterogeneous attenuating media, highlighting its potential as a valuable tool for studying inverse problems in viscoelastic media.

Multicomponent seismic data can enhance seismic imaging of complex geologic structures, lithology classification, and hydrocarbon reservoir characterization. In the multicomponent data processing workflow, reconstructing velocities for P and S waves ([Formula: see text] and [Formula: see text], respectively) is a critical step. We previously develop an efficient eikonal equation-based adjoint state characteristic reflection traveltime tomography (ASCRT) method, which automatically picks reflection traveltimes and maps reflectors using fast kinematic imaging for acoustic media. In this study, this ASCRT method is further extended to handle elastic velocities. Here, [Formula: see text] and [Formula: see text] are hierarchically inverted using reflection traveltimes of primary PP and PS waves, which are more challenging to automatically pick. The independent picking of PP and PS waves is demonstrated through a combination of radiation pattern analysis and migration tests. In addition, we develop a novel [Formula: see text] inversion method for ASCRT, using traveltimes of PS converted waves. The traveltime calculation and the eikonal equation-based fast kinematic imaging technique can be used for PS converted waves. Furthermore, an efficient eikonal equation-based migration method for PP and PS waves is developed. Double-layer and multilayer velocity model inversions with a layer-stripping strategy demonstrate that the elastic ASCRT method effectively reconstructs reliable background elastic velocities. The migrated images indicate significant improvements after elastic ASCRT inversion. Our elastic ASCRT is applied to ocean-bottom data acquired in the East China Sea, illustrating that the reconstructed background elastic velocities match the well-log velocities and that the image quality of PP and PS waves is significantly improved.

This study develops a nonstationary blind deconvolution using spectral constraints (NBDSC), a novel method to address the ill-posed nature of wavelet inversion in nonstationary seismic data. NBDSC does not require prior knowledge about the [Formula: see text]-factor of the medium and source wavelet, which are typically unknown. Our power spectrum-guided waveform inversion method corrects a wavelet to match any desired spectral shape. This method can be introduced as a wavelet-smoothing constraint, yielding blind deconvolution with spectral constraints (BDSC). First, based on a segmented stationary convolution model (SSCM), the original nonstationary convolution problem is decomposed into multiple stationary subproblems. Next, an absorption-constrained wavelet power spectrum inversion method is used on the SSCM decomposition results to extract the time-varying wavelet power spectrum. Finally, the estimated time-varying wavelet power spectrum is used as a spectral constraint to perform BDSC on each subseismic trace. We use synthetic and field data examples to demonstrate how the NBDSC enhances the resolution of nonstationary seismic data. Compared with conventional blind deconvolution methods, NBDSC incorporates prior knowledge of the wavelet power spectrum, thereby reducing inversion ambiguity and improving noise robustness.

Despite their low porosity and permeability, tight gas-bearing sandstones contain substantial oil and gas trapped within their dense matrix. These unconventional reservoirs have been developed through advanced drilling and hydraulic fracturing technologies, unlocking new energy supplies and offering a sustainable solution to the growing global demand for energy. However, seismic exploration of these reservoirs presents unique challenges. The tight nature of the sandstones means that conventional seismic methods often struggle to provide the necessary details for successful exploration. Special interpretation strategies are required to accurately identify potential hydrocarbon-bearing zones due to the subtle contrasts in acoustic impedance between the tight sandstones and the surrounding beds. In our exploration of deep tight gas-bearing sandstone from northwestern China, we find that the matrix shear modulus effectively distinguishes between sandstone and mudstone lithologies, making it a useful parameter for describing tight gas-bearing sandstone. We develop a novel amplitude-variation-with-angle (AVA) inversion method to obtain matrix shear modulus. First, under deep tight geologic conditions with porosity less than 10%, the modulus of relatively cleaner sandstone reservoirs, saturated with a two-phase mixture of gas and water, is approximately equivalent to the dry modulus. Based on this, we approximate the saturated modulus in the Gray AVA approximation as the dry modulus. Next, we introduce the Nur critical porosity model to establish an AVA approximation using the matrix modulus, achieving the Aki-Richards approximation accuracy. Then, using this approximation, we establish a Bayesian inversion method for the matrix modulus. Theoretical model testing and practical data inversion demonstrate that this method is theoretically sound and practically feasible. After estimating Nur’s critical porosity parameter using logging data and core measurements, the inverted matrix shear modulus parameters successfully describe the vertical and lateral distribution of the deep tight sandstone.

Seismic petrophysical inversion aims to predict petrophysical properties, such as porosity, lithology, and fluid saturation, given a set of seismic measurements by combining inverse methods with amplitude-variation-with-offset (AVO) and rock-physics models. Stochastic optimization methods are often used; however, their applicability to real cases is limited by large computational costs due to the size of the data and the required number of iterations. We develop a deterministic approach based on the Gauss-Newton inversion using the Levenberg-Marquardt algorithm. This nonlinear optimization algorithm calculates the Hessian matrix using the sum of the outer products of the gradients. The novelty of the work is the analytical calculation of the Jacobian matrix of the seismic and rock-physics models. In our implementation, the partial derivatives of the seismic properties with respect to the petrophysical properties are expressed in a closed form. In our implementation, we test the inversion with a seismic AVO convolutional model combined with the soft sand model; however, our inversion method can be extended to any geophysical model. Our method has been tested and validated on a simulated synthetic data set based on field data measured offshore Norway. The results prove that our method is robust and efficient in predicting petrophysical properties from seismic data.

The classic paper by Brown and Korringa, giving the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid, has been understood for almost 50 years to be a generalization (to nonhomogeneous solid) of earlier work on a similar topic by Gassmann. However, this understanding is based on a single passage of their paper, which contains a logical error, and so is not valid. The primary result of Brown and Korringa, which contains an additional parameter (compared with Gassmann), is correct as written in 1975 for all porous isotropic rocks (with one or many minerals) with uniform pore pressure, and stationary statistics of microgeometry. Their parameter [Formula: see text] is understood to be the Mean compressibility, not the matrix compressibility. Their parameter [Formula: see text] is not required, being eliminated in favor of the overall solid compressibility [Formula: see text]. The extension of the isotropic argument to bulk anisotropic rocks is straightforward (but with a complication due to thin layering). Their results should be used (in place of Gassmann) for analysis of all 4D seismic differences, if the differences are shown, by future experimentation, to be significant in practice.

Traditional deterministic inversion algorithms provide a single inverse model with limited appraisal of the uncertainty in the model parameters. Transdimensional Bayesian inversion provides a solution to address this issue. We develop a hierarchical transdimensional Bayesian inversion scheme for interpreting 2D magnetotelluric data. The subsurface resistivity model is represented using Voronoi cells, and the noise level in the data is treated as an unknown parameter and is estimated during the inversion process. The flexible parameterization and noise estimation make the scheme completely data driven. To increase the efficiency of forward computation, a multiresolution (MR) scheme is implemented. In this scheme, the grid is formed by stacking grids with different horizontal resolutions. The use of the MR approach enhances the efficiency of forward computations without any significant loss in the accuracy of the solutions. The transdimensional inversion is applied to two synthetic data sets to examine the resolution of complex model structures. In addition, the algorithm is applied to the COPROD2 data set. The results demonstrate that the transdimensional inversion effectively images complex subsurface structures and accounts for the unknown noise level in the data. The generated model ensemble provides valuable insights into the subsurface resistivity model through various statistical measures, with direct uncertainty quantification via standard deviation.

In this paper, we develop an average-derivative based second-order staggered-grid finite-difference (ADSSFD) scheme for performing 3D frequency-domain elastic wave modeling. The ADSSFD scheme applies the idea of average derivative to discretize the first-order derivative of the stress tensors and combines the anti-lumped mass strategy to obtain a set of optimal finite-difference (FD) coefficients. We combine the Levenberg-Marquardt method with simulated annealing to minimize the misfit between the normalized phase velocity and one, ensuring the solution avoids local minima. Dispersion analysis indicates that when the upper limit of the phase velocity error is set to 1%, the ADSSFD scheme can reduce the minimum number of grid points per shear wavelength to 3.6. Numerical experiments indicate that under the same numerical setting, the ADSSFD scheme exhibits higher accuracy than the classic fourth-order staggered-grid FD (FSFD) scheme and the compact second-order staggered-grid FD scheme. In addition, we find the feasibility of using the ADSSFD scheme to perform simulations in the 3D overthrust model. We use a parallel iteration solver to solve the large sparse linear system obtained after the discretization. Because the condition number of the matrix in the ADSSFD scheme is smaller than that of the FSFD scheme, the ADSSFD scheme is cheaper to run.

The magnetotelluric (MT) method is highly advantageous in geothermal resource exploration because of its deep probing capability, low cost, and sensitivity to key geothermal system elements, such as heat sources, reservoirs, and caprocks. To address the challenge of balancing the MT inversion resolution and computational efficiency, we develop a Gauss-Newton optimization inversion algorithm based on the solution space dimensionality reduction (SSDR) technique. This algorithm uses SSDR to establish algebraic mapping relationships between the edges of grid cells, effectively reducing the degrees of freedom in the linear equation system required for forward and pseudoforward modeling, thereby saving computational costs. In addition, during the inversion process, a multiple right-hand direct-iterative hybrid solver is used to solve the dimensionally reduced linear equations, further enhancing the computational efficiency of the inversion. We then validate and analyze the accuracy and efficiency of the forward and inversion algorithms using various synthetic models. Finally, we apply this inversion strategy to field MT data from the Dabie Mountain geothermal area in eastern China. On the basis of the inversion results, we then analyze the thermal genesis model of the region, providing technical support for the exploration and development of deep geothermal resources.