ABSTRACT Fourier transfer function FTF(f), defined as a Fourier transform of the impulse response of the medium between a source and receiver, fully characterizes transfer properties in the frequency domain. Modulus of FTF(f) is one of the most used characteristics of seismic site response in both theoretical and numerical analyses of effects of local surface structures on seismic motion. It is especially useful for identifying resonance peaks and/or the amplified frequency bands. By definition, however, |FTF(f)| cannot provide information on temporal development of the site response. In this article, we introduce a novel characterization of seismic site response—time–frequency transfer function TFTF(t, f) that properly characterizes temporal development of frequency content of transfer function. This feature significantly helps in interpretating wavefield composition. We present theory that is readily applicable to any time-dependent output signal. We demonstrate superior properties of TFTF(t, f) in comparison with standard |FTF(f)| for selected quasi-1D canonical models of local surface sedimentary structures. In the model of a semi-infinite layer, TFTF(t, f) made it possible to distinguish modes of vertical resonance and modes of surface waves, which is impossible based on standard |FTF(f)|.

When developing a finite-difference (FD) scheme, one of the key aspects that must be addressed is the spatial discretization of material parameters and the implementation of material interfaces. Mittet (2017) and Moczo et al. (2022) recently suggested a novel approach based on the wavenumber limitation of the medium. They demonstrated that, due to spatial discretization, a model of the medium must be wavenumber-limited by a wavenumber km smaller than the Nyquist wavenumber. Mittet (2021a) and Valovcan et al. (2024) proved that the wavefield (numerically simulated or exact) in a km-limited medium can only be accurate up to km/2. Here, we numerically demonstrate a perfect subcell resolution (capability to sense the position of interface within a grid cell) of FD modeling based on the wavenumber-limited medium using a finite spatial low-pass filter. The finding that it is possible to use a finite-length filter for wavenumber limitation of the medium is of key importance for the next development of the concept in terms of computational efficiency. We demonstrate an unprecedented accuracy for a canonical model—a material interface between two homogeneous half-spaces. The time–frequency envelope and phase misfits between the FD solution and the exact analytical solution are surprisingly small and are the same for any position of the interface within a grid spacing. We show that the misfits for the reflected and transmitted waves are solely due to grid dispersion. We compare the accuracy of the FD solution for the wavenumber-limited medium with that for the harmonically averaged modulus (Moczo et al., 2002). The FD modeling based on the wavenumber-limited medium is considerably more accurate. The proof of concept of the wavenumber limitation of the medium should be followed by development of a practical way of wavenumber limitation using finite spatial filters in 2D and 3D elastic problems.

SUMMARY Seismic ambient-noise horizontal-to-vertical spectral ratios (H/V, HVSR) are widely used to characterize near-surface structure and to identify site resonant frequencies. In this article, we propose a non-parametric statistical descriptor of ambient-noise HVSRs. In single-station practice, HVSR curves are typically represented by the geometric mean of window-wise ratios at each frequency, that is, by the median of a lognormal model. However, distributions of window ensembles are frequently skewed or multimodal. The lognormal median and symmetric uncertainty bands can then depart from the most probable amplitudes and distort the spread. The lognormal mode can mitigate this bias in near-lognormal cases, but it remains tied to a unimodal parametric form. As a consequence, many workflows rely on strict window selection or rejection to justify lognormal assumption. To address this problem, we present a data-driven, non-parametric alternative that provides a marked improvement when ensembles deviate from lognormality. At each frequency, we estimate the probability density of the log-amplitude distribution using kernel density estimation (KDE), transform it to linear space and take the mode of the linear-domain density as the representative curve. Uncertainty is quantified by highest density intervals, which naturally accommodate asymmetry and multibranch distributions. Using an example of a near-lognormal microtremor, we demonstrate that the KDE mode closely tracks the lognormal mode while revealing a systematic upward shift of the commonly used lognormal median. Using a microtremor at a structurally complex site, we observe strongly multimodal amplitude and directional statistics. Both the lognormal-median and lognormal-mode curves fall between competing modes. However, the KDE-based modes and intervals follow the dominant branches and capture multimodal spread. Because the density is inferred directly from window-wise data, the method reduces the need for strong window rejection. We also compare the KDE-mode descriptor with the energy-ratio estimator which averages horizontal and vertical component energies before taking their ratio. Agreement between the energy-ratio estimator and the KDE mode is consistent with a stable single HVSR population, whereas discrepancies help identify frequency bands affected by multimodality, non-stationarity or directional effects. Finally, we extend the same framework to directional HVSR using the horizontal spectral matrix and quantify directional variability through a $\pi $-periodic circular KDE of the axial principal horizontal directions. We show that KDE-based modes, highest density intervals and directional spread provide robust distribution-aware observational diagnostics that can guide data selection, uncertainty assignment and frequency weighting in site-characterization and inversion workflows.

Numerical simulations of earthquakes and seismic wave propagation require accurate material models of the solid Earth. In contrast to purely elastic rheology, poroelasticity accounts for pore fluid pressure and fluid flow in porous media. Poroelastic effects can alter both the seismic wave field and the dynamic rupture characteristics of earthquakes. For example, the presence of fluids may affect cascading multifault ruptures, potentially leading to larger-than-expected earthquakes. However, incorporating poroelastic coupling into the elastodynamic wave equations increases the computational complexity of numerical simulations compared to elastic or viscoelastic material models, as the underlying partial differential equations become stiff. In this study, we use a Discontinuous Galerkin solver with Arbitrary High-Order DERivative time stepping of the poroelastic wave equations implemented in the open-source software SeisSol to simulate 3-D complex seismic wave propagation and 3-D dynamic rupture in poroelastic media. We verify our approach for double-couple point sources using independent methods including a semi-analytical solution and a finite-difference scheme and a homogeneous full-space and a poroelastic layer-over-half-space model, respectively. In a realistic carbon capture and storage reservoir scenario at the Sleipner site in the Utsira Formation, Norway, we model 3-D wave propagation through poroelastic sandstone layers separated by impermeable shale. Our results show a sudden change in the pressure field across material interfaces, which manifests as a discontinuity when viewed at the length scale of the dominant wavelengths of S or fast P waves. Accurately resolving the resulting steep pressure gradient dramatically increases the computational demands, requiring high-resolution modelling. We show that the Gassmann elastic equivalent model yields almost identical results to the fully poroelastic model when focusing solely on solid particle velocities. We extend this approach using suitable numerical fluxes to 3-D dynamic rupture simulations in complex fault systems, presenting the first 3-D scenarios that combine poroelastic media with geometrically complex, multifault rupture dynamics and tetrahedral meshes. Our findings reveal that, in contrast to modelling wave propagation only, poroelastic materials significantly alter rupture characteristics compared to using elastic equivalent media since the elastic equivalent fails to capture the evolution of pore pressure. Particularly in fault branching scenarios, the Biot coefficient plays a key role in either promoting or inhibiting fault activation. In some cases, ruptures are diverted to secondary faults, while in others, poroelastic effects induce rupture arrest. In a fault zone dynamic rupture model, we find poroelasticity aiding pulse-like rupture. A healing front is induced by the reduced pore pressure due to reflected waves from the boundaries of the poroelastic damage zone. Our results highlight that poroelastic effects are important for realistic simulations of seismic waves and earthquake rupture dynamics. In particular, our poroelastic simulations may offer new insights on the complexity of multifault rupture dynamics, fault-to-fault interaction and seismic wave propagation in realistic models of the Earth's subsurface.

SUMMARY We present a time-domain distributional finite-difference scheme based on the Lebedev staggered grid for the numerical simulation of wave propagation in acoustic and elastic media. The central aspect of the proposed method is the representation of the stresses and displacements with different sets of B-splines functions organized according to the staggered grid. The distributional finite-difference approach allows domain-decomposition, heterogeneity of the medium, curvilinear mesh, anisotropy, non-conformal interfaces, discontinuous grid and fluid–solid interfaces. Numerical examples show that the proposed scheme is suitable to model wave propagation through the Earth, where sharp interfaces separate large, relatively homogeneous layers. A few domains or elements are sufficient to represent the Earth’s internal structure without relying on advanced meshing techniques. We compare seismograms obtained with the proposed scheme and the spectral element method, and we show that our approach offers superior accuracy, reduced memory usage, and comparable efficiency.

ABSTRACT Analysis of equations of motion by Moczo et al. (2022) led to the conclusion that the discrete (grid) representation of the heterogeneous medium must be wavenumber bandlimited up to the Nyquist frequency. This is a consequence of the spatial discretization. Mittet (2021a) reported that if the discrete grid model of medium coincides with the true medium up to some wavenumber, the simulated wavefield is accurate only up to a half of this wavenumber. Here, we present results of the systematic and comprehensive analysis focused on the principal limits of accuracy of numerically simulated wavefields. First, we analyze wavenumber spectra of (1) exact wavefields in a heterogeneous elastic medium, (2) wavenumber bandlimited wavefields, and (3) spatially discretized wavefields. Then, we derive spatial dependence of the frequency spectrum of waves generated by a finite source, and perturbing wavefields due to a small perturbation of the medium and due to a small wavenumber bandlimited perturbation of the medium. We analyze an interaction of an incoming wave with the medium perturbation through a change of phase difference and through wavenumber spectra. We draw conclusions on the wavenumber limitation of wavefields in the wavenumber bandlimited heterogeneous medium. We numerically verify the fundamental finding using exact solutions. The main consequence for the finite-difference (FD) modeling based on spatial discretization of the computational domain is: Due to spatial sampling, the medium must be wavenumber limited up to the Nyquist frequency. Then, the wavefield should not be sampled by less than four spatial grid spacings per shortest wavelength to obtain sufficiently accurate results. This applies to any heterogeneous FD scheme.

By analyzing the equations of motion and constitutive relations in the wavenumber domain, we gain important insight into attributes determining the accuracy of finite-difference (FD) schemes. We present heterogeneous formulations of the equations of motion and constit-utive relations for four configurations of a wavefield in an elastic isotropic medium. We Fourier-transform the entire equations to the wavenumber domain. Subsequently, we apply the band-limited inverse Fourier transform back to the space domain. We analyze conse-quences of spatial discretization and wavenumber band limitation. The heterogeneity of the medium and the Nyquist-wavenumber band limitation of the entire equations has important implications for an FD modeling: The grid representation of the heterogeneous medium must be limited by the Nyquist wavenumber. The wavenumber band limitation replaces spatial derivatives both in the homogeneous medium and across a material inter-face by continuous spatial convolutions. The latter means that the wavenumber band limi-tation removes discontinuities of the spatial derivatives of the particle velocity and stress at the material interface. This allows to apply proper FD operators across material interfaces. A wavenumber band-limited heterogeneous formulation of the equations of motion and constitutive relations is the general condition for a heterogeneous FD scheme.

ABSTRACT It is well known that higher-order and thus longer-stencil finite-difference operators (FDOs) can be advantageously used for evaluating spatial derivatives in the finite-difference schemes applied to smoothly heterogeneous media. This is because they reduce spatial grid dispersion. However, realistic models often include sharp material interfaces. Can high-order long-stencil FDOs be applied across such material interface? We address this question by comparing exact spatial derivatives against derivatives approximated by FDOs with respect to the interface representation, velocity contrast, and order of the FDO. The interface is considered in an arbitrary position with respect to the spatial grid. The material interface exactly represented by the Heaviside step function causes a large error of the FDO spatial derivative near the interface. The maximum error near the interface practically does not depend on the order of the FDO. There are only small differences in errors among FDOs of different orders elsewhere. The larger the velocity contrast, the larger the error. If the material interface is represented using a wavenumber band-limited Heaviside function, the error is smoothed and several times smaller. The error in the wavenumber band-limited model decreases with an increasing order of the FDO. Our findings combined with those by Moczo et al. (2022) lead to the important conclusion: The wavenumber band-limited representation of the material interface is not only a necessary consequence of discretization of the original physical model but also significantly reduces the error in evaluating a spatial derivative using the FDO.

We present a new methodology of the finite-difference (FD) modelling of seismic wave propagation in a strongly heterogeneous medium composed of poroelastic (P) and (strictly) elastic (E) parts. The medium can include P/P, P/E and E/E material interfaces of arbitrary shapes. The poroelastic part can be with (i) zero resistive friction, (ii) non-zero constant resistive friction or (iii) JKD model of the frequency-dependent permeability and resistive friction. Our FD scheme is capable of subcell resolution: a material interface can have an arbitrary position in the spatial grid. The scheme keeps computational efficiency of the scheme for a smoothly and weakly heterogeneous medium (medium without material interfaces). Numerical tests against independent analytical, semi-analytical and spectral-element methods prove the efficiency and accuracy of our FD modelling. In numerical examples, we indicate effect of the P/E interfaces for the poroelastic medium with a constant resistive friction and medium with the JKD model of the frequency-dependent permeability and resistive friction. We address the 2-D P-SV problem. The approach can be readily extended to the 3-D problem.

Many applications from the fields of seismology and geoengineering require simulations of seismic waves in porous media. Biot's theory of poroelasticity describes the coupling between solid and fluid phases and introduces a stiff reactive source term (Darcy's Law) into the elastodynamic wave equations, thereby increasing computational cost of respective numerical solvers and motivating efficient methods utilising High-Performance Computing. We present a novel realisation of the discontinuous Galerkin scheme with Arbitrary High-Order DERivative time stepping (ADER-DG) that copes with stiff source terms. To integrate this source term with a reasonable time step size, we utilise an element-local space-time predictor, which needs to solve medium-sized linear systems – each with 1,000 to 10,000 unknowns – in each element update (i.e., billions of times). We present a novel block-wise back-substitution algorithm for solving these systems efficiently, thus enabling large-scale 3D simulations. In comparison to LU decomposition, we reduce the number of floating-point operations by a factor of up to 25, when using polynomials of degree 6. The block-wise back-substitution is mapped to a sequence of small matrix-matrix multiplications, for which code generators are available to generate highly optimised code. We verify the new solver thoroughly against analytical and semi-analytical reference solutions in problems of increasing complexity. We demonstrate high-order convergence of the scheme for 3D problems. We verify the correct treatment of point sources and boundary conditions, including homogeneous and heterogeneous full space problems as well as problems with traction-free boundary conditions. In addition, we compare against a finite difference solution for a newly defined 3D layer over half-space problem containing an internal material interface and free surface. We find that extremely high accuracy is required to accurately resolve the slow, diffusive P-wave at a or near a free surface, while we also demonstrate that solid particle velocities are not affected by coarser resolutions. We demonstrate that by using a clustered local time stepping scheme, time to solution is reduced by a factor of 6 to 10 compared to global time stepping. We conclude our study with a scaling and performance analysis on the SuperMUC-NG supercomputer, demonstrating our implementation's high computational efficiency and its potential for extreme-scale simulations.