MeV ultrafast electron diffraction (UED) is a widely used technique for ultrafast structural dynamics studies of matter in numerous areas. The development of the laser wakefield accelerator (LWFA) shows great potential for an advanced all-optical electron source based on LWFA in UED applications. Here we experimentally demonstrated that an LWFA-based device with a miniaturized permanent magnet beamline can generate and manipulate electron beams suitable for UED. During beam transport, the LWFA electron beams with intrinsically short duration undergo temporal stretching owing to the energy spread and are subsequently compressed by the following double-bend achromat. The optimized double-bend achromat can make the beamline isochronous such that the arrival time jitter induced by the shot-to-shot energy fluctuation can be eliminated, and allow the advantage of the natural laser-beam synchronization for LWFAs to emerge. With the energy filtering, the beam energy spread can be reduced to 3% (full-width at half-maximum), while a sufficient amount of charge (11.9 fC) per bunch for diffraction is retained. Using a laser-driven terahertz deflector, the beam length and arrival time jitter measured at the sample location are approximately 49.6 fs (root mean square (r.m.s.)) and 4.7 fs (r.m.s.), respectively, resulting in a temporal resolution of ~49.8 fs. Comprehensive start-to-end simulations indicate the potential of reducing the bunch length to ~10 fs (r.m.s.) with a lower energy spread of around 1.6%. Clear single-shot and multi-shot diffraction patterns of single-crystalline gold samples are obtained, and the derived lattice constant agrees well with the actual value. Our proof-of-principle experiments open the door to the detection of ultrafast structural dynamics using MeV LWFA beams, and pave the way for UED applications with sub-10 fs temporal resolution. Researchers demonstrate that a laser wakefield accelerator-based device with a miniaturized permanent magnet beamline can generate and manipulate electron beams suitable for ultrafast electron diffraction.

Frequency-domain viscoelastic anisotropic full waveform inversion provides a powerful framework for high-resolution subsurface characterization but remains computationally challenging at realistic scales. Although 2.5D formulations offer a favourable compromise between accuracy and computational cost, they introduce complex interactions among spectral decomposition, sparse linear system solutions, and sensitivity kernel evaluation, thereby complicating efficient implementation. This study presents the modernization of the SEIS25D_GQG FORTRAN90 code for 2.5D frequency-domain viscoelastic anisotropic full-waveform inversion through a scalable hybrid-parallel architecture. Although the underlying formulation has been established previously, the original implementation was limited to small-scale applications because of its sequential architecture and large memory footprint. The redesigned framework combines distributed-memory (MPI) and shared-memory (OpenMP) parallelism to exploit multiple levels of concurrency, including wavenumber-sample distribution, data-sample distribution, and model-block parallelism. The redesign also enables flexible solver selection, supporting both banded and sparse solvers, with the latter significantly reducing memory requirements. Performance evaluations conducted on a high-performance computing platform using synthetic models of increasing complexity show that forward modeling benefits from solver-level optimizations and wavenumber-level parallelism, while Fréchet derivative computation, identified as the dominant computational cost, achieves near-linear scaling. Together, these improvements reduce inversion runtimes by more than two orders of magnitude compared with the legacy implementation. Finally, the framework is validated using a modified Marmousi model, demonstrating the capability of the redesigned code to handle realistic-scale problems.

Traveltime computation is an effective way to simulate seismic wave propagation in isotropic and anisotropic media. It often requires the phase and group velocities along a ray direction. The phase and group velocities are not functions of the ray direction, but functions of the slowness direction, which motivates the need for efficient numerical approaches for computing the slowness direction. A generalized method was proposed to achieve this goal in 3D tilted transverse isotropic media, which involves solving a nonlinear equation of two unknowns, the inclination and azimuthal angles of the slowness direction, and is computationally demanding. To overcome the difficulty, the equation was recast by projecting it onto the local coordinates that align with the axis of symmetry of the medium. Under the local coordinate system, the nonlinear equation reduces to an equation of one unknown, the inclination angle, with the azimuthal angle obtained from the given ray direction. Hence, the nonlinear equation can be solved efficiently and its solutions can be used to calculate the phase and group velocities. As an application, the proposed method was used to compute traveltimes in 3D vertical transverse isotropic (VTI) media. The calculated group velocities are incorporated into the fast sweeping method to solve the original anisotropic Eikonal equation directly on uniform meshes. The feasibility of the proposed group velocity calculation method is evaluated in three 3D homogeneous anisotropic models, and the application on traveltime computation is tested on a 3D homogeneous VTI model and the British Petroleum anisotropic VTI model.

The attenuation of seismic waves holds significant theoretical and practical importance for characterizing subsurface lithology and predicting oil and gas.Vertical Seismic Profile(VSP)data can integrate with well log and seismic data,facilitating the high-precision delineation of subsurface structures.However,conventional approaches using VSP data often focus solely on the downgoing wavefield,neglecting the upgoing wavefield due to its energy dominance,exceeding 90％,in the VSP spectrum.This paper proposed a joint viscoacoustic waveform inversion of upgoing and downgoing wavefields in zero-offset VSP.The separated upgoing and downgoing wavefields are obtained through decoupled forward modeling of VSP upgoing and downgoing waves.A joint inversion is then performed by establishing an objective function that effectively utilizes the transmission waves from the downgoing wavefield and the reflection waves from the upgoing wavefield,thereby enhancing the accuracy and convergence of viscoacoustic full-waveform inversion.Objective function analysis of different wavefields to velocity and quality-factor reveals varying contributions of these wavefields to the inversion of viscoelastic parameters.Finally,numerical experiments and field data inversion verify the feasibility of the joint inversion method,laying a useful tool for oil and gas analysis.

The Fr & eacute;chet derivatives or sensitivity kernels of the observed seismograms are fundamental to seismic full-waveform inversion (FWI). They quantitatively measure the seismogram variations caused by any physical parameter perturbation of the Earth's subsurface. The 3-D viscoelastic tilted transversely isotropic (TTI) media are often encountered in practices due to the presence of dip thin layers, joints, fractures or cracks, orientated grains or crystallization, and water or gas saturation. To image such subsurface, we have derived explicit 3-D frequency-domain Fr & eacute;chet derivatives of the seismogram spectrum with respect to 13 independent physical parameters of TTI rock, which include density, five elastic moduli, five Q-factors, and inclination and declination angles of the symmetric axis of rock structure. We have demonstrated a fully parallel implementation to compute the Fr & eacute;chet derivatives and conduct synthetic subsurface imaging experiments, in which the 13 independent parameters of the subsurface targets are successfully reconstructed. The experimental results have verified the correctness and validity of the derived 3-D Fr & eacute;chet derivatives for imaging viscoelastic TTI media.

Accurate and efficient modeling method is important for the study of seismic wave propagation and crucial for the full waveform inversion. High-order temporal numerical methods can improve the performance of viscoacoustic wave modeling at the aspects of computational accuracy and efficiency. Based on the staggered Adams–Bashforth time stepping scheme, we present the high-order recursive convolution method, which is obtained by the Taylor series expansion and offers arbitrary order by the number of terms retained, to calculate the temporal convolutions for viscoacoustic wave modeling. The theoretical analysis of the high-order recursive convolution method demonstrates the achievement of high accuracy of viscoacoustic wave modeling, and the proposed method beats the conventional high-order auxiliary-differential-equation method in terms of the computer memory cost and runtime. The 1D and 2D numerical examples that verify the superiority of the proposed method. In a 2D case with a 1000×000 grid, the fourth-order recursive convolution method reduces computation time to 73% and memory requirements to 88% compared to the fourth-order auxiliary differential equation method, while maintaining the same accuracy level. Numerical tests on the Marmousi model demonstrate the accuracy and feasibility of the high-order recursive convolution method for heterogeneous models.

This research demonstrates an innovative numerical technique to simulate seismic wave propagation of a practical point source in complex 2-D geological models, which encompass a free surface topography, an undulating seafloor, and acoustic, elastic isotropic, viscoacoustic and viscoelastic anisotropic rocks. This technique is particularly beneficial in scenarios where 3-D wave modeling is resource-intensive and may efficiently offer the 3-D wavefields from arbitrary 2-D geological models often encountered in practice. Based on the point-source viscoelastic wave equations in a 2-D heterogeneous tilted transversely isotropic (TTI) medium, representative of subsurface igneous and sedimentary rocks, we tailor the wave equations valid for different rocks and the boundary conditions of the free-surface topography and seafloor and adapt the conventional memory variable method and the newly developed Taylor-series recursive convolution method to solve such point-source comprehensive wave equations. To overcome the inherent computational intensity of the methods, we convert the complex domain into a real domain and implement a fully parallelized computing strategy to ensure that the runtime of the numerical simulation remains on par with that of common 2-D wave modeling. Our experimental validations confirm the accuracy of the Taylor-series recursive method to offer the 3-D wavefields in an arbitrary heterogeneous 2-D geological model having a free-surface topography or an undulating seafloor. Moreover, our applications of this technique to two benchmark practical 2-D geological models demonstrate its capability to replicate 3-D wavefields in arbitrary viscoelastic anisotropic media, and greatly help in interpreting offshore and onshore seismic data and generating an accurate image of the subsurface.

This paper uses the plane waves in super-space and real space to give the mathematic definitions of the slowness vector, ray-velocity vector, traveltime and raypath of seismic body waves in a viscoelastic anisotropic media (VEAM). Then we explain their physical implications in terms of wave propagation. We show that the slowness vector may be decomposed into the homogeneous and inhomogeneous components, the ray-velocity vector must be a homogeneous complex vector, and the real and imaginary traveltimes represent the wave-phase propagation and wave-energy attenuation, respectively, and the raypaths of wave-phase propagation and wave-energy attenuation are different but both follow Fermat' principle.

Solving large sparse linear systems in 3D frequency-domain seismic wave modeling, especially in viscoelastic anisotropic media, poses significant challenges due to the increasing number of discrete moduli and nonzero elements in the linear system matrix. The computational load surpasses that of acoustic or viscoacoustic media, making it even more challenging when dealing with multisource problems. Popular scientific tools for solving a linear system, such as multifrontal massively parallel sparse direct solver (MUMPS), structured matrix package (STRUM- PACK), and portable extensible toolkit for scientific computation (PETSc), can be used, but their applicability to our specific problem has not been comprehensively evaluated. Our study aims to tackle the challenges in solving large, sparse, complex-valued symmetric linear systems with multiple right-hand side vectors for 3D frequency-domain seismic wave modeling. We have adopted the preconditioned conjugate gradient iterative algorithms as the foundation for our research, introducing two highly cost-effective parallel iterative solvers: the parallel symmetric successive overrelaxation conjugate gradient (P-SSORCG) and the parallel incomplete Cholesky conjugate gradient (P-ICCG). These novel solvers are subjected to a comprehensive comparative analysis against well-established scientific tools, such as MUMPS, STRUMPACK, and PETSc, in the context of 3D frequency-domain seismic wave modeling. We show their promising performances in a practical 3D SEG/EAGE overthrust model and demonstrate that the grouped P-SSORCG offers an efficient alternative to parallel direct solvers, particularly in situations wherein computational resources are limited.

Seismic wavefield forward modeling in anelastic (attenuating) media is a fundamental tool for both data processing and interpretation in modern seismic exploration. We propose a generalized recursive convolution (RC) formula to calculate the temporal convolutions directly, rather than solving many auxiliary partial differential equations of the memory variables when dealing with a viscoelastic medium. The new formula is obtained in terms of the Taylor series expansion and offers approximations of the convolutions to arbitrary order by the number of terms retained. We conduct theoretical and numerical comparisons of the new method with the commonly used memory variable method and other traditional RC methods. The comparisons show that the new method has the highest accuracy of all these RC methods and yields better performance with various stress relaxation times and time steps than the common leapfrog time-stepping scheme to solve the auxiliary partial differential equations of the memory variables. Our numerical examples verify the versatility and feasibility of the new method for viscoacoustic and viscoelastic wave modeling.

Under the high-frequency assumption, the slowness vector in a viscoelastic anisotropic medium is often defined by a complex-valued vector, whose direction is given by a complex unit vector that can hardly ever be explained by the intuitive physical reasoning related to the actual seismic wave propagation direction. We extend the existing conjugate real ray-tracing (C-RRT) method to compute the ray velocity vectors and then apply it to determine the slowness vectors for three body waves (qP, qSV, and qSH) in a viscoelastic anisotropic medium. Moreover, we dissect the slowness vector with two physical specifications — traveltime gradient and inhomogeneity components, to reveal the physical significance of the slowness vector and give illustrative examples of these components and the phase propagation (real traveltime) and wave-energy attenuation (imaginary traveltime) wavefronts of seismic waves. We also display the homogeneous and inhomogeneous components, as well as the inhomogeneity angles (or deviation angles between these two wavefronts) for the three body waves (qP, qSV, and qSH) in shale with different quality factors ( Q-factors). These results reveal the elusive link between the homogeneous slowness vector and the Q-factors of the medium. The theoretical and numerical dissections of the slowness vector and the extended C-RRT method are extremely helpful in understanding seismic wave propagation and offer a ray-tracing method for wave-energy compensation in seismic migration and reconstruction of subsurface images through seismic ray tomography in viscoelastic anisotropic media.

Computation of the ray-velocity vector is crucial in seismic ray tracing for the three body waves (qP, qSV, qSH) in viscoelastic anisotropic media. The primary challenge is dealing with the likely cusps or triplications of the qSV wavefronts, which makes it theoretically difficult to track qSV raypaths and the reflection and transmission of these body waves in such media. We review three traditional methods, namely, g-Hamiltonian, p-Hamiltonian, and explicit c-derivative, and then present two new approaches called implicit c-derivative and g*-Hamiltonian to tackle the challenges of seismic ray tracing. We theoretically prove the equivalence of these five methods to calculate the group-velocity vector in a viscoelastic anisotropic medium, apply these methods to some rock samples, and investigate the applicability of each method. Our results indicate that if the body wave is homogeneous (i.e., its propagation and attenuation wavefronts are parallel to each other) or if the body wavefront has no cusps or triplications, then all the methods offer a consistent solution of the ray-velocity vector. If the body wave is inhomogeneous (its propagation and attenuation wavefronts are at different angles or cross each other) and cusps and triplications occur in the wavefronts, then all the methods but the g*-Hamiltonian one fail to give the proper solution of the ray-velocity vector. This work demonstrates that the innovative g*-Hamiltonian method is the only approach to overcome the theoretical difficulty of seismic ray tracing in viscoelastic anisotropic media.

Accurate seismic wave modeling of viscoelastic anisotropic medium is a fundamental tool for seismic data processing, interpretation and full waveform inversion. Also, free water surface, topographic relief and irregular seabed are often encountered in practical seismic surveys. Thus, basing on the General Maxwell Body, we proposed a generalized matrix form of the velocity-stress seismic wave equation, which becomes valid for composite viscoelastic anisotropic media and satisfies the boundary conditions in presence of topographic free surfaces and irregular fluid–solid interfaces. We theoretically show that the viscoelastic effect of a medium may be considered as the intrinsic body sources accumulated in wavefield history and computed by a recursive convolution formula accurately and efficiently. We also demonstrated that such a generalized viscoelastic wave equation may be solved with the curvilinear MacCormack finite difference method and validated the accuracy and feasibility of the proposed method. The modeling results in homogeneous and heterogeneous media match well with the analytical solutions and the references yielded by the spectral element solutions.

Many rocks exhibit electrical anisotropic characteristics, leading to artifacts in isotropic inversion of dc resistivity data. To mitigate this issue, we employ the flexible unstructured finite element (FE) method for 3-D anisotropic forward modeling and anisotropic inversion for dc resistivity data. Synthetic inversions further validate algorithm feasibility. The code not only replicates the inversion of principal axis resistivities accomplished by previous researchers but also reconstructs challenging all angle parameters. Comparative analysis of inversion using isotropic and arbitrary anisotropic assumptions reveals that anisotropic inversion accurately recovers the parameters. Models comprising conductive and resistive targets, respectively, are utilized to stimulate two realistic scenarios. The angle parameters of conductive target are challenging to recover compared to a conductive one with the same anisotropy. We conduct comprehensive experiments to assess the results across various data types. Borehole-generated data enhances the principal resistivities resolution, aligning with the direction of the lines connecting to the boreholes, rather than all ones, whereas conventional surface configurations are unable to accurately capture discrepancies between principal resistivities and Euler angles. Moreover, the investigations confirm the resolvability of each anisotropy component, especially angle parameters, which depend on borehole distribution. These experiments on data-acquisition configurations greatly enhance the practical possibility of successfully collecting data containing anisotropy information, achieving reconstruction of angular parameters, and obtaining better images of principal axis resistivity in practical applications. Even with the correct anisotropic assumption, unreasonable data-acquisition configurations produce inversion results with distortion. Testing with the model containing multiple anisotropic targets further confirms the algorithm's effectiveness.

Point-source to line-source transformation for seismic data prior to 2-D full waveform inversion is recommended for efficiently producing high-resolution images of the subsurface. The traditional filter of the transformation is derived by comparison of 3-D and 2-D Green’s function inhomogeneous acoustic media, so it doesn’t guarantee accurate transformation in complex heterogeneous viscoelastic anisotropic media. We adhere to a straightforward filter for compensating amplitude, incorporating offsets and phase shifts, and adjusting the stretching factors to rescale amplitude differences across various components. Additionally, we devise new stretching factors to accommodate time-domain transformations. Testing the proposed method with multi-layer and Marmousi models validates better amplitude compensation in both near and far offsets than the traditional hybrid method.

Crossing conjugate normal faults（CCNFs） are extensively developed in many hydrocarbon-producing basins, generally existing in the form of incomplete CCNFs. Nevertheless, the effect of the non-conjugate zone of the CCNFs on the conjugate relay zone post late tectonic action has not been previously studied. We use 3D elastic-plastic modeling to investigate the influence of incomplete（i. e.,partially intersecting） CCNFs on the pattern of deformation of strata in the intersection region. A series of model simulations were performed to examine the effects of horizontal tectonic extension, fault size, and fault depth on the deformation of conjugate relay zones of incomplete CCNFs. Our analyses yielded the following results.（1） The model of incomplete conjugation predicts a convex-up style of deformation in the conjugate graben region superimposed on overall subsidence under applied horizontal tectonic extension.（2） The degree of convex-up deformation of the conjugate graben depends on the influence of the non-conjugate zone on the conjugate relay zone, which varies with the amount of horizontal tectonic extension, fault size, and fault burial depth.（3） Our results indicate that incomplete CCNFs can form convex-up deformation, similar to that in the Nanpu Sag area and provide a sound understanding of hydrocarbon migration and accumulation.

A g*-Hamiltonian method for tracing real rays was developed that can handle cusps and triplication of quasi-shear wave in a general viscoelastic anisotropic medium. We demonstrate that the g*-Hamiltonian method can produce homogeneous ray-velocity vectors (with parallel real and imaginary parts) and the slowness vectors of reflected and transmitted waves at the interface based on the real Snell's law (RSL), which constrains only the continuity of the real parts of the slowness vectors, or the real slowness direction (RSD) method, which ignores the inhomogeneous component of the slowness vector. These methods are based on the characteristic lines with different Hamiltonians. Our research indicates that these methods are limited to pre-critical incidence ranges. Moreover, we derived a complex energy velocity vector (energy flux velocity) expression, which is always homogenous. We found that directions of corresponding energy velocity calculated with complex Snell's law (CSL) at a contact of two viscoelastic anisotropic materials well match the solutions of the RSL and RSD methods for all R / T waves except post-critical incidence in which the RSL and RSD methods fail to obtain homogenous ray velocities. The RSL and RSD methods result in discrepancies in the ray quality factor, R / T coefficients, and energy ratios, especially for post-critical incidence. However, the exact critical angle determined by the RSD method approximates the 'critical' angle for anelastic/inhomogeneous waves, which was a previous challenge. Our calculations suggest that the energy velocity and energy quality factor obtained with the CSL method can be used for real ray tracing at the interface of viscoelastic anisotropic media, and the complex energy flux velocity vector is always exactly homogeneous. For the post-critical incidence, the RSL and RSD methods fail because the ray quality factor drastically changes from the infinite down to near 2, which contradicts homogeneous ray velocity even in elastic anisotropic materials for RSD method. However, the energy flux quality factor for the elastic-anisotropic material is all infinite, even for post-critical incidence, which is correct. We also show that the CSL method has the same efficiency as the RSD method.

Full-waveform inversion (FWI) is an appealing data-fitting approach to produce high-resolution seismic imaging of subsurface geological structures. Implementation of FWI in frequency domain usually requires spectral data of seismograms only at some frequencies, making it more convenient. Compared to the common approximations such as acoustic or viscoacoustic, 3D frequency-domain viscoelastic FWI is quite challenging due to the fact that more additional parameters are needed to be inverted resulting in highly increased computational cost and unwanted crosstalk between quality factors and wave velocities. We derived explicit Fréchet derivatives for sensitivity kernels for 3D frequency-domain viscoelastic FWI and validated it numerically. Then, we performed 3D multiparameter viscoelastic FWI in frequency domain in parallel environments. Our results demonstrate successful inversion of five key parameters—density, p-wave velocity, s-wave velocity, and associated quality factors—thus providing a preliminary validation of the feasibility and efficacy of our approach for multiparameter FWI.