Geophysical Journal International

SummaryVelocity macro-model building is an essential step of the seismic imaging workflow. Indeed, obtaining acceptable results through migration or full waveform inversion is highly dependent on the kinematic accuracy of the background/initial velocity model. Two decades ago, stereotomography was proposed as an alternative to reflection traveltime tomography, the first relying on semi-automatic picking of locally coherent events associated with small reflection or diffraction segments tied to scatterers in depth by a pair of rays, while the latter on interpretive picking of laterally-continuous reflections. The flexibility of stereotomography paved the way for many developments that have shown the efficiency of the method whilst emphasizing on the complementary information carried out by traveltimes and slopes of locally-coherent events. A recent formulation recast stereotomography under a matrix-free formulation based on eikonal solvers and the adjoint-state method. In the latter, like in the previous works, the scatterer positions and the velocity field are updated jointly to tackle the ill-famed velocity-position coupling in reflection tomography. Following on from this adjoint-state formulation, we propose a new parsimonious formulation of slope tomography that offers the chance to restrain the problem to minimizing the residuals of a single data class being a slope, in search of a sole parameter class being the subsurface velocity field. This parsimonious formulation results from a variable projection, which is implemented by enforcing a consistency between the scatterer coordinates and the velocity macro-model through migration of kinematic attributes. We explain why the resulting reduced-parametrization inversion is more suitable for tomographic problems than the most common joint inversion strategy. We benchmark our method against the complex Marmousi model along with a validation through time domain full waveform inversion and then present the results of a field data case study.

SummaryA promising way to perform seismological studies in the Arctic region is deploying seismic stations on ice floes. The pioneering works by the Alfred Wegener Institute (AWI) Bremerhaven (Schlindwein et al. 2007; Laderach and Schlindwein, 2011) have demonstrated the efficiency of such floating networks to explore local and regional seismicity and to build 3D seismic models. However, problems remain, related to the identification of different types of seismic waves, particularly S-waves. Here, we perform 2D and 3D numerical simulations of seismic waves emitted by an earthquake to explore the possibility of recording different phases on the sea surface. We use different types of simple shear source models, namely, strike-slip, vertical displacement and normal faults. In the calculated wave field, we obtain three major types of seismic waves recorded on the sea surface: Pw, Sw and SPw (w denotes an acoustic wave in the water layer) and numerous multiple waves. The clarity of the recorded phases strongly depends on the type of wave, source mechanism, epicentral distance, thickness of the water layer and depth of the source. For example, the Pw phase is clearest for the strike-slip mechanism, less clear for the normal fault and almost invisible for the vertical displacement. The Sw phase is observable in all of these cases; however, it can be confused with the SPw phase that arrives earlier. In addition, at some distances, the Sw wave interferes with the multiple Pw2 wave and therefore is hardly detectable. In summary, the numerical simulations in a model with a water layer have demonstrated several non-obvious features of wave propagation that should be taken into account when analysing experimental data recorded on ice floes.

SummaryFrom a suite of 56 chemically-driven dynamo simulations with aspect ratio χ (inner to outer core radii) ranging from 0.10 to 0.44, we conduct the first systematic investigation of the impact of inner-core size on the reversing behaviour of dynamos. We show that the growth of the inner core leads to a transition between a “small inner-core” regime (χ ≤ 0.18), when the field produced is intermediately strong and dipolar, and a “large inner-core” regime (χ > 0.26), when the field is stronger and more dipolar. During that transition the field is weaker and slightly less dipolar. For aspect ratios 0.20 ≤ χ ≤ 0.22, reversal frequencies may be more sensitive to changes in the vigour of the convection, allowing high frequencies to be reached much more easily. Although other factors than the size of the inner core likely contribute to controlling the reversal frequency of the Earth’s dynamo, we hypothesise that the occurrence of such a transition for the Earth’s core between the end of the Precambrian and the end of the Devonian could possibly account for the manifestation of an unusual long-lasting episode of predominantly reversal hyperactivity and complex low intensity fields during that still poorly documented period of time.

SummaryStoneley waves propagating in a borehole will produce reflected waves when they encounter a porous formation. In previous studies, the simplified Biot-Rosenbaum theory is used to calculate the Stoneley wave reflection coefficient for porous formations. Such simplified theory ignores the effect of formation frame elasticity, thus can't obtain accurately the Stoneley wave reflection for the porous formation with small stiffness. In this study, to take the effect of formation frame elasticity into account, we use the Biot's theory in the low-frequency limit to simulate the Stoneley wave reflection by employing the velocity-stress finite-difference time-domain (FDTD) method. In addition, as permeability of the formation varies in the axial direction, and the viscosity of the pore fluid also changes in a radial direction due to mud invasion, it is difficult to use the simplified theory to calculate the reflection coefficient of the Stoneley wave in such heterogeneous cases. Therefore, this study first investigates the effect of the formation permeability heterogeneity on Stoneley wave reflection coefficient by this FDTD method.The FDTD method is verified by a comparison with the real axis integration (RAI) method with respect to the Stoneley wave propagation in a borehole surrounded by a homogeneous porous formation. The reflection coefficient obtained by the FDTD method is smaller than that using the simplified theory, which shows that elasticity of the formation frame affects Stoneley wave reflection: the effect of elasticity on the reflection coefficient is greater when the formation frame is less rigid or when the porosity and permeability of the formation are lower. According to the simulation results of the FDTD method, a modified simplified theory which can improve the calculation accuracy of Stoneley wave reflection coefficient is proposed. Furthermore, the effects of the permeability heterogeneity on the Stoneley wave reflection are investigated: the reflection coefficient peaks change when permeability alters in an axial direction, and the peak interval increases. For the mud invasion model, the reflection coefficient is almost identical to that of the homogeneous model, which has the same permeability as the borehole wall of the mud invasion model.

SummarySlow slip events (SSEs) and other forms of slow and transient fault slip are becoming increasingly recognized as important, due to their influence on seismicity and potential to provide information on fault zone properties. The majority of our current knowledge of slow fault slip has been obtained from geodetic and seismologic field measurements, but laboratory data of slow slip are comparatively scarce. Here, I present the results of laboratory friction experiments conducted on 11 natural samples from major fault zones around the world, all obtained by scientific drilling. The experiments are conducted water saturated, at room temperature and 10 MPa effective normal stress, representative of in-situ conditions on the shallow portions of fault zones from which the samples were recovered and where slow slip is known to occur. A key component of these experiments is shearing at realistically slow driving rates of cm/yr, accurately simulating tectonic driving rates. In most samples, these cm/yr driving rates produce laboratory SSEs, which are instances of accelerating slip accompanied by a stress drop. The peak slip velocities and stress drops measured in these laboratory SSEs are comparable with those of natural SSEs measured or estimated from geodetic data. A strong correlation is observed between reduced pre-SSE velocity and higher peak slip velocity for the entire laboratory SSE dataset. In contrast to the velocity data, significant scatter is observed in the per cent stress drop measurements. The source of this scatter can be attributed to samples with a significant expandable clay component, which tend to exhibit larger stress drops. Results of velocity-stepping tests at cm/yr rates show a tendency for velocity-weakening friction not observed at higher sliding velocities, and that the materials with lower values of the rate-dependent friction parameter a-b tend to produce faster SSEs. Critical stiffness analyses within the framework of rate-and-state friction laws show that most of the SSEs observed in this study do not satisfy the condition for slip instability. The SSEs are more consistent with accelerating stable slip, although the stiffness condition allowing such behavior is not always satisfied. Considering the laboratory SSEs to be accelerating stable slip, I present a conceptual model for their nucleation. Key elements of the model are a healing-dominated departure from steady-state causing partial locking and velocity decrease, followed by a transition to a velocity-dominated phase representing the actual slip event. The model is consistent with observations from geodetic measurements and the experimental observations in this study. In general, characteristics of SSE-producing fault portions such as the ability to strengthen and store elastic strain energy released as stress drops may be expected to enhance coseismic slip from remotely nucleating earthquakes, an effect which may be quite limited but should be investigated further.

SummaryNumerical modelling of seismic waves is a method for simulating the propagation of waves in the earth. The objective is to make a prediction of the seismogram when given an assumed structure of the subsurface. The diffusive-viscous theory can be used to describe the attenuation of seismic waves propagating in fluid-saturated rocks, and to study the relationship between the frequency dependence of reflections and fluid-saturation in a porous medium, since the generally used theories such as acoustic, elastic theory, among others, are unable to effectively characterise the subsurface of the earth. We derive the second-order Runge-Kutta finite-volume scheme for the diffusive-viscous wave equation and based on this scheme, we simulate the propagation of seismic waves in a fluid-saturated medium. The numerical results indicate that the propagating waves in the fluid-saturated media attenuate greatly when compared with those of the acoustic scenario.

SummaryWe present a new approach to estimate 3-D seismic velocities along a target interface. This approach uses an artificial neural network trained with user-supplied geological and geophysical input features derived from both a 3-D seismic reflection volume and a 2-D wide-angle seismic profile that were acquired from the Galicia margin, offshore Spain. The S-reflector detachment fault was selected as the interface of interest. The neural network in the form of a Multi-Layer Perceptron was employed with an autoencoder and a regression layer. The autoencoder was trained using a set of input features from the 3-D reflection volume. These set of features included the reflection amplitude and instantaneous frequency at the interface of interest, time-thicknesses of overlying major layers, and ratios of major layer time-thicknesses to the total time-depth of the interface. The regression model was trained to estimate the seismic velocities of the crystalline basement and mantle from these features. The ‘true’ velocities were obtained from an independent full-waveform inversion along a 2-D wide-angle seismic profile, contained within the 3-D dataset. The autoencoder compressed the vector of inputs into a lower dimensional space, then the regression layer was trained in the lower dimensional space to estimate velocities above and below the targeted interface. This model was trained on 50 networks with different initializations. Thirty-seven of these networks reached minimum achievable error of 2 per cent. The low standard deviation (< 300 m/s) between different networks and low errors on velocity estimations demonstrate that the input features were sufficient to capture variations in the velocity above and below the targeted S-reflector. This regression model was then applied to the 3-D reflection volume where velocities were predicted over an area of ∼400 km2. This approach provides an alternative way to obtain velocities across a 3-D seismic survey from a deep non-reflective lithology (e.g. upper mantle) where conventional reflection velocity estimations can be unreliable.

SummarySeismic inversion is essentially ill-posed because of inaccurate and insufficient data. Existing methods for improving stability and reducing the size of solution space are usually model-driven. These methods often rely on sophisticated mathematical models associated with prior knowledge or expectation of the structure, distribution, and correlation of subsurface parameters. In general, model-driven methods are established only for particular usage situations, which are known to have poor adaptivity. Also, because of their abstraction and simplification, these methods often have difficulty achieving satisfactory accuracy and resolution for complex geology. In this paper, we present a data-driven inversion framework using cooperative sparse representation (CSR) for multi-parameter simultaneous inversion. We first let the system learn a joint dictionary from well-log data using multi-component dictionary learning algorithm. Such a learned dictionary captures structural features and correlations of the elastic parameters. Then, in the inversion process, we replace the constraints of mathematical models with that of cooperative sparse representation over the learned dictionary on model parameters. This multi-parameter inversion framework incorporates the correlation among different parameters and is data-driven. It only relies upon the widely accepted assumption that subsurface parameters have a certain degree of consistency and lateral continuity within the same geological formation. Tests on a complex synthetic example and field data demonstrate that our algorithm is effective in improving the resolution and accuracy of solutions, especially for the density parameter that is insensitive to amplitude but extremely important for reservoir identification.

SummaryThe structure and processes of collisional orogens at the edges of cratons are complex and the extent to which the crust is affected by subduction of cratonic material, slab break-off and post-collisional magmatism is poorly constrained. The Carpathian arc is an Alpine orogenic system, involving the collision and subduction of the East European Craton and the closure of the Tethys ocean. Unusual Neogene back-arc magmatism developed after subduction ended ∼9 Ma ago and upper-mantle Mw>7 earthquakes are associated with a relic slab actively breaking-off in the Vrancea Zone. To investigate the crustal and uppermost mantle structure, we analyse teleseismic earthquakes recorded at 56 broadband seismic stations in Romania and Moldova. Using phase-weighted H − κ stacking of receiver functions, we estimate the bulk crustal thickness and Vp/Vs ratio, a good discriminant of crustal composition. This technique was used on all stations and assumes a single-layer crust with a discontinuity at its base. Additionally, by jointly inverting receiver functions and ambient noise we obtain shear wavespeed (Vs) models to ∼60 km depth beneath 34 stations and provide another independent measure of the Moho depth and its apparent seismic sharpness. Crustal thickness estimates from the two methods are inconsistent in regions where the Moho is gradational or an intra-crustal discontinuity is more dominant. We support this interpretation with a range of synthetic tests. Above the actively detaching Vrancea slab, multiple intracrustal discontinuities and a gradational Moho are identified on Vs profiles, and high Vp/Vs (>1.83) are consistent with magmatic underplating. On the back-arc side of the slab, magmatism-pervaded crust has low Vp/Vs signature (<1.73), consistent with felsic composition or the presence of water without extensive partial melt. In the back-arc Transylvanian Basin and Apuseni Mountains, the upper mantle has anomalously low Vs (<4 km/s to 60 km depth), the crust is dominantly felsic and thin (<35 km), with a Moho probably inherited from the obducted Neo-Tethys ophiolites. On the fore-arc side, sedimentary basins have low Vs (<2.5 km/s) and high Vp/Vs signatures (>1.8), indicating a high fluid-pressure regime. Finally, Precambrian units in the foreland are seismically heterogenous, with Moho depths of 30 km-50 km, and extend ∼100 km beneath the South and Southeast Carpathians, implying that the Vrancea slab underlies Precambrian-aged crust.

SummaryThe formation of fold-thrust belts at convergent margins is a dynamic process. Accretion of weak sediments to the front of the overriding plate results in crustal thickening and continued flexural subsidence of the underthrusting plate. Fold-thrust belts are often treated as a Coulomb wedge having self-similar geometries with a critical taper, and either a rigid or isostatically compensated base. In this paper we build upon this work by developing a new dynamic model to investigate both the role of the thickness and material properties of the incoming sediment, and the flexure in the underthrusting plate in controlling the behaviour and evolution of fold-thrust belts. Our analysis shows that the evolution of fold-thrust belts can be dominated by either gravitational spreading or vertical thickening, depending on the relative importance of sediment flux, material properties and flexure. We apply our model to the Makran accretionary prism and the Indo-Burman Ranges, and show that for the Makran flexure must be considered in order to explain the dip of the sediment-basement interface from seismic reflection profiles. In the Indo-Burman Ranges, we show that incoming sediment thickness has a first-order control on the variations in the characteristics of the topography from north to south of the Shillong Plateau.

SummaryWe propose a new Bayesian method to reveal the Vs structure of the near surface of the earth using spatial autocorrelation (SPAC) functions and apply this new method to synthetic, broadband, and geophone datasets. The principle of SPAC is introduced, and an implementation of the Bayesian Monte Carlo inversion (BMCI) for modeling SPAC coherency functions is described. To demonstrate its effectiveness, BMCI is applied to synthetic tests, data from 14 SPAC array sites in the Salt Lake Valley (SLV), Utah, and two arrays (one broadband and one geophone) located in south central Utah. The Vs models derived from previous SPAC analysis of the 14 SLV sites differ by 10 per cent at most from those determined by BMCI and lie within uncertainties determined for the BMCI models. These agreements demonstrate the effectiveness of the BMCI method. The synthetic tests and applications to the SLV SPAC data show BMCI has great potential to resolve Vs structure down to at least 400 m. To achieve resolution for deeper Vs structure, longer duration deployments, wider array apertures, and additional seismometers or geophones can be employed. Additionally, when the target frequencies are greater than 0.1 Hz, there is no apparent disadvantage in using geophone data for BMCI compared to broadband data. Most significantly, BMCI places a quantifiable constraint on the uncertainties of the Vs models as well as Vs30.

SummaryKnowledge of the coastal ocean mass variations is important for understanding the ocean climate and sea level change. However, estimates of coastal ocean mass variations have been perplexed due to poor representativeness of previous Gravity Recovery and Climate Experiment (GRACE) data across the land/ocean boundary. We here use GRACE mascons solutions to investigate the coastal ocean mass variations (within 400 km band) at global and regional scales. GRACE mascons solutions are advanced by spatial constraints and leakage correction. For the global mean, it is found that mascons ocean mass variations are in rough agreement with those inferred from satellite altimeter and an ocean reanalysis; the agreement is better on seasonal scales (6.4 ± 0.5 mm annual amplitude for the mascons and 7.3 ± 0.7 mm annual amplitude for the inferred); on the other hand, large differences are shown for the linear trend (2.1 ± 0.1 mm/yr for the mascons and 3.2 ± 0.1 mm/yr for the inferred), indicating that it is far from closing the coastal sea level budget. At regional scales, high consistency between mascons and the inferred is observed over two shallow areas i.e. Europe coast and East China coast, whose annual amplitudes are 24 ± 5 and 42 ± 3 mm, respectively. Ocean mass variations are not well captured by mascons over other regions (e.g. Africa coast, Indian coast, and North and South America coast). Possible explanations include: (1) the steric component plays a more important role in sea level, consequently resulting in weak mass signal; and (2) the performance of mascons varies with locations. We find that mascons solutions show better variability (especially seasonal cycles) than do the traditional spherical harmonic coefficients, suggesting that mascons applied with spatial constraints improve the coastal ocean mass variability.

SummaryIn the absence of quantitative relationships, a cross-gradient function can be used to correlate unknown physical properties in a joint inversion of geophysical data sets. It introduces a structural correlation between properties. A commonly used, inexact approach adds a weighted cross-gradient term as a penalty to the cost function being minimized during inversion. This weighting factor needs to be tuned to balance the regularization and cross gradient terms. In this paper we propose non-linear weighting for the cross-gradient function which addresses the very different magnitudes of the cross gradient and regularization terms. This approach also couples the weighting factors for the regularization and correlation terms reducing the number of tuning parameters. The approach is investigated for a synthetic case. Results are also shown for the 3D joint inversion of high-resolution magnetic and gravity anomaly data from Southern Queensland in Australia with over 30 Million cells.

SummaryThe converted wave data (P-to-s/S-to-p), commonly termed as receiver functions, contain noises of various origins. Such noises may influence the modeling and may sometimes lead to over interpretations of the data. In order to suppress noise, we use a robust sparsity enhancing tool, i.e. the Seislet Transform (ST), to process receiver function data by applying regularization in the seislet domain. The transform utilizes the multiscale orthogonal basis and the basis-functions are oriented along the dominant seismic phases following local linearity. The inversion results of both the synthetic and field examples from the Hi-CLIMB network and station HYB from the Indian shield show an excellent performance over the original data sets.

SummaryWe present a method for measuring the dispersion of higher mode surface wave phase velocities from a single seismogram using a hierarchical transdimensional Bayesian approach. The 1-D shear velocity profiles down to 800km depth between sources and seismic stations are regarded as the controlling parameters to tune the phase velocities of fundamental and higher modes. The misfits between synthetics and real waveforms indicate whether the phase velocities are recovered well from the data. We use Monte Carlo Markov chains (MCMC) to approximate the posterior distribution of each model parameters, and assess the uncertainties from these probability density functions. These techniques can test models of varying dimensions while being parsimonious, thereby letting the data themselves control the complexity of the solution. Another advantage is that the algorithm can decide how much data noise is needed in order to explain the data without overfitting them. The data noise can be treated as an unknown and different noise levels can be applied to the different time windows considered. The posterior noise distributions can then be used as an indicator of the quality of the waveform fit within each frequency-time window. We considered phase velocities between 50s and 200s for each mode, and performed a reliability analysis to determine which modes and periods are reliably constrained. In this paper, we first present the method and demonstrate its feasibility with synthetic tests, which show that the technique is robust. We then illustrate it with applications to real data. We applied the method to two paths sampling Australia using earthquakes at regional distances, and obtained results that agree well with previous studies. The new method can be used in regional and global tomographic studies to obtain phase velocity maps and 3-D models of seismic velocities and anisotropy at depths that are not well resolved by fundamental mode surface waves or body waves.

SummaryThis paper develops an analytical solution for the transient response of a semi-infinite solid with a fluid surface subjected to an arbitrary line source in the solid. The line source is first decomposed into P- and S-pulses, and the wave fields for transient responses of the fluid and semi-infinite solid are presented. Applying the Fourier and Laplace transform methods, the analytical solution in the transform domain, which is an integral solution, is obtained. Using de Hoop's method, a suitable distortion of the contour is provided to change the path of integration and the positions of the singularities and branch points, and a new form of integration is obtained. As the new integration has a form similar to that of the Laplace transform pair, the inverse Laplace transform is directly obtained, and the analytical solution in the time domain is developed. In practice, the compressional wave velocity in the solid is usually larger than that in the fluid. For the cylindrical S-pulse line source, the head wave can always be observed in the fluid-solid system. For the cylindrical P-pulse line source, the head wave can be observed in the fluid-solid system only when the wave is refracted from the solid to the fluid. Numerical examples are provided to discuss the behaviour of the fluid-solid system.

SummaryContinuous long-term sound sources are recorded at hydroacoustic station H03S, a three-element hydrophone array south of Robinson Crusoe Island, between 23 April 2014 and 20 August 2017. The origin of the signal between 3-20Hz is investigated by using cross-correlation, array processing using plane wave beamforming, and spectral analysis. One-bit normalization is successfully applied as a cross-correlation preprocessing step in order to suppress undesired earthquake events in the data. Travel times retrieved from averaged cross-correlations do not yield a coherent array direction of arrival. Averaged envelopes of the cross-correlations, however, indicate a coherent signal approaching H03S from a south-southwest direction. Beamforming indicates two dominant back azimuth directions: 172°-224° (Antarctica) and 242° (Monowai Volcanic Seamount). This continuous source field creates possibilities to investigate the applicability of acoustic thermometry at hydrophones H03 S1-S2. Cross-correlation and array processing indicate significant directional variation in local modal propagation, most likely related to the steep slope in the bathymetry near H03S. In addition, it is demonstrated that the ambient noise field is not sufficiently equipartitioned. It is shown that this causes a large error in the estimated temperature, primarily due to the short receiver spacing. These large errors have not been addressed in previous studies on deep-ocean acoustic thermometry. Hence, it is shown that acoustic thermometry does not perform well on small arrays such as H03S. The power spectral density yields a strong broadband signal in January - March, most likely related to iceberg noise. A narrow banded signal around 17Hz between April and September corresponds to whale calls. The best-beam sound pressure levels towards Antarctica are compared to ERA5 climatologies for sea ice cover and normalized stress into the ocean, supporting the hypothesis of iceberg noise.

SummaryIt is difficult to separate additive random noise from spatially coherent signal using a rank-reduction method that is based on the truncated singular value decomposition (TSVD) operation. This problem is due to the mixture of the signal and the noise subspaces after the TSVD operation. This drawback can be partially conquered using a damped rank reduction method, where the singular values corresponding to effective signals are adjusted via a carefully designed damping operator. The damping operator works most powerfully in the case of a small rank and a small damping factor. However, for complicated seismic data, e.g., multi-channel reflection seismic data containing highly curved events, the rank should be large enough to preserve the details in the data, which makes the damped rank reduction method less effective. In this paper, we develop an optimal damping strategy for adjusting the singular values when a large rank parameter is selected so that the estimated signal can best approximate the exact signal. We first weight the singular values using optimally calculated weights. The weights are theoretically derived by solving an optimization problem that minimizes the Frobenius-norm difference between the approximated signal components and the exact signal components. The damping operator is then derived based on the initial weighting operator to further reduce the residual noise after the optimal weighting. The resulted optimally damped rank reduction method is nearly an adaptive method, i.e., insensitive to the rank parameter. We demonstrate the performance of the proposed method on a group of synthetic and real five-dimensional seismic data.

SummaryThe accuracy and efficiency of numerical simulations of seismic wave propagation and earthquake ground motion in realistic models strongly depend on discrete grid representation of the material heterogeneity and attenuation. We present a generalization of the orthorhombic representation of the elastic medium (Kristek et al. 2017) to the viscoelastic medium to make it possible to account for a realistic attenuation in a heterogeneous viscoelastic medium with material interfaces. An interface is represented by an averaged orthorhombic medium with rheology of the Generalized Maxwell body (GMB-EK, equivalent to the Generalized Zener body). The representation is important for the possibility of applying one explicit finite-difference scheme to all interior grid points (points not lying on a grid border) no matter what their positions are with respect to the material interface. This is one of the key factors of the computational efficiency of the finite-difference modelling. Smooth or discontinuous heterogeneity of the medium is accounted for only by values of the effective (i.e. representing reasonably averaged medium) grid moduli and densities. Accuracy of modelling thus very much depends on how the medium heterogeneity is represented/averaged. We numerically demonstrate accuracy of the developed orthorhombic representation. The orthorhombic representation neither changes the structure of calculating stress-tensor components nor increases the number of arithmetic operations compared to a smooth weakly heterogeneous viscoelastic medium. It is applicable to the velocity-stress, displacement-stress and displacement FD schemes on staggered, partly-staggered, Lebedev and collocated grids. We also present an optimal procedure for a joint determination of the relaxation frequencies and anelastic coefficients.

SummaryWe introduce a novel block rational Krylov method to accelerate three-dimensional time-domain marine controlled-source electromagnetic modelling with multiple sources. This method approximates the time-varying electric solutions explicitly in terms of matrix exponential functions. A main attraction is that no time stepping is required, while most of the computational costs are concentrated in constructing a rational Krylov basis. We optimize the shift parameters defining the rational Krylov space with a positive exponential weight function, thereby producing smaller approximation errors at later times and reducing iteration numbers. Furthermore, for multi-source modelling problems, we adopt block Krylov techniques to incorporate all source vectors in a single approximation space. The method is tested on two examples: a layered seafloor model and a 3D hydrocarbon reservoir model with seafloor bathymetry. The modelling results are found to agree very well with those from 1D semi-analytic solutions and finite-element time-domain solutions using a backward Euler scheme, respectively. Numerical benchmarks show that the block method benefits from better memory and cache efficiency, resulting in about 1.26 to 1.48-fold speedup compared to non-block methods. Further efficiency gains are achieved through optimized rational Krylov techniques, allowing our approach to outperform classical time stepping schemes.