Dispersion error analysis can help to assess the performance of finite-element discretizations of the wave equation. Although less general than the convergence estimates offered by standard finite-element error analysis, it can provide more detailed insight as well as practical guidelines in terms of the number of elements per wavelength needed for acceptable results. We present eigenvalue and eigenvector error estimates for cubic Hermite elements on an equidistant 1-D mesh and on a regular structured 2-D triangular mesh consisting of squares cut in half. The results show that in 1D, the spectrum consists of 2 modes. If these are unwrapped, the spectrum is effectively doubled. The eigenvalue or dispersion error stays below 7% across the entire spectrum. The error in the corresponding eigenvectors, however, increases rapidly once the number of elements per wavelength decreases to one. In terms of element size, the dispersion error is of order 6 and the eigenvector error of order 4. The latter is consistent with the classic finite-element error estimate. In 2D, we provide eigenvalue and eigenvector errors as a series expansion in the element size and obtain the same orders. 2-D numerical tests in the timeand frequency-domain are included.@en

Non-periodic homogenization has proved to be an accurate asymptotic method for
computing long-wavelength equivalent media for the seismic wave equation, turning
small-scale heterogeneities and geometric complexity into smooth elastic properties.
Using homogenized media allows i) decreasing the computation cost of wave propagation
simulation and ii) studying the apparent, small-scale-induced anisotropy. After illustrating
these two aspects briefly, we propose to analyze in great detail the accuracy of body waves
simulated in homogenized 3D models of the subsurface. First, the behaviour of head-,
reflected and refracted waves with respect to source-receiver o\set, maximum frequency
and velocity contrast across a planar interface, is investigated. Then, we consider the SEGEAGE overthrust model to exemplify how the accuracy of simulated body waves anticorrelates with the distance to seismic source and the amount of apparent anisotropy. In
high apparent anisotropy regions, we show that the first-order correction provided by the
homogenization theory significantly improves the computed wavefield. The overall results
of this analysis better frame the use of homogenized media in seismic wave simulation.@en

To estimate the depth errors in a subsurface model obtained from the inversion of seismic data, the stationary-phase approximation in a two-dimensional constant-velocity model with a dipped reflector is applied to migration with a time-shift extension. This produces two asymptotic solutions: one is a straight line, and the other is a curve. If the velocity differs from the true one, a closed-form expression of the depth error follows from the depth and apparent dip of the reflector as well as the position of the amplitude peak at a non-zero time shift, where the two solutions meet and the extended migration image focuses. The results are compared to finite-frequency results from a finite-difference code. A two-dimensional synthetic example with a salt diapir illustrates how depth errors can be estimated in an inhomogeneous model after inverting the seismic data for the velocity model.

@en

The accuracy of a model obtained by full-waveform inversion can be estimated by analysing the sensitivity of the data to perturbations of the model parameters in selected subsurface points. Each perturbation requires the computation of the seismic response in the form of Born scattering data for a typically very large number of shots, making the method time consuming. The computational cost can be significantly reduced by considering the point where the subsurface parameters are perturbed as a Born scatterer. Instead of modelling each shot separately, reciprocity relations provide the Green functions from the sources to the scatterer in terms of Green’s functions from the scatterer to the sources. In this way, the Born scattering data from a single point in the isotropic elastic case for a marine acquisition with pressure sources and receivers can be expressed in terms of the Green functions for force and moment tensor sources located at the scatterer and only a small number of forward runs are required. A 2-D example illustrates how the result can be used to determine the hessian and local covariance matrix for the model parameters at the scatterer at the cost of 5 forward simulation.@en

The accuracy of a model obtained by multi-parameter full-waveform inversion can be estimated by analysing the sensitivity of the data to perturbations of the model parameters in selected subsurface points. Each perturbation requires the computation of the seismic response in the form of Born scattering data for a typically very large number of shots, making the method time consuming. The computational cost can be significantly reduced by placing sources of different types at the Born scatterer, the point where the subsurface parameters are perturbed. Instead of modelling each shot separately, reciprocity relations provide the wavefields from the shot positions to the scatter point in terms of wavefields from the scatterer to the shot positions. In this way, the Born scattering data from a single point in the isotropic elastic case for a marine acquisition with pressure sources and receivers can be expressed in terms of the wavefields for force and moment tensor sources located at the scatterer and only a small number of forward runs are required. A two-dimensional example illustrates how the result can be used to determine the Hessian and local relative covariance matrix for the model parameters at the scatterer at the cost of five forward simulations. In three dimensions, that would be nine.

@en

To investigate the physical processes behind induced seismicities due to, for example, production of hydrocarbons from a reservoir, most of the earlier studies performed geomechanical simulations on a simple reservoir geometry. The effect of fluid depletion is, in general, simulated for such a simple geometry. Neglecting the contribution of realistic 3-D reservoir geometries can lead to a wrong estimation of the incremental stress field. A reliable estimate of the induced stress field is key to producing meaningful simulation results. We perform geomechanical simulations on a simple fault model as well as a more realistic model based on the known geological structures at the earthquake source-region in Zeerijp region, the Netherlands. Our results demonstrate that the angle of the fault intersection affects the incremental stress field, including the effective normal stress, the shear stress, and hence, the Coulomb stress and the SCU value. Our results also show a shift in the rupture pattern and the location of the maximum slip on the fault plane. We conclude that, to properly evaluate the effects of production activities and to simulate precisely the in-situ stress field and the induced seismicity, the incorporation of a realistic reservoir structure in modelling is essential.@en

The uncertainty of model parameters obtained by full-waveform inversion can be determined from the hessian of the least-squares error functional. Because the hessian is generally too costly to compute and too large to be stored, a segmented representation of perturbations of the reconstructed subsurface model in the form of geological units is proposed. This enables the computation of the hessian and the related covariance matrix on a larger length scale. A synthetic 2-D isotropic elastic example illustrates how conditional and marginal uncertainties can be estimated for the properties per geological unit by themselves and in relation to other units. A discussion on how the chosen length scale affects the result is included.@en

The stationary-phase method applied to migration with a time-shift extension in a 2-D constant-velocity model with a dipped reflector produces two solutions in the domain of the extended image: one a straight line and the other a curve. If the velocity differs from the true one, the depth error follows from the depth and apparent dip of the reflector as well as the depth of the amplitude peak at a non-zero time shift, where the two solutions meet and the extended image focuses. The results are compared to finite-frequency results from a finite-difference code. A 2-D synthetic example with a salt diapir illustrates how depth errors can be estimated in an inhomogeneous model after inverting the seismic data for the velocity model.@en

Finite elements with mass lumping allow for explicit time stepping when modelling wave propagation and can be more efficient than finite differences in complex geological settings. In two dimensions on quadrilaterals, spectral elements are the obvious choice. Triangles offer more flexibility for meshing, but the construction of polynomial elements is less straightforward. The elements have to be augmented with higher-degree polynomials in the interior to preserve accuracy after lumping of the mass matrix. With the classic accuracy criterion, triangular elements suitable for mass lumping up to a polynomial degree 9 were found. With a newer, less restrictive criterion, new elements were constructed of degree 5–7. Some of these are more efficient than the older ones. To assess which of all these elements performs best, the acoustic wave equation is solved for a homogeneous model on a square and on a domain with corners, as well as on a heterogeneous example with topography. The accuracy and runtimes are measured using either higher-order time stepping or second-order time stepping with dispersion correction. For elements of polynomial degree 2 and higher, the latter is more efficient. Among the various finite elements, the degree-4 element appears to be a good choice.

@en

Temporal dispersion correction of second-order finite-difference time stepping for numerical wave propagation modelling exploits the fact that the discrete operator is exact but for the wrong frequencies. Mapping recorded traces to the correct frequencies removes the numerical error. Most of the implementations employ forward and inverse Fourier transforms. Here, it is noted that these can be replaced by a series expansion involving higher time derivatives of the data. Its implementation by higher-order finite differencing can be sensitive to numerical noise, but this can be suppressed by enlarging the stencil. Tests with the finite-element method on a homogeneous acoustic problem with an exact solution show that the method can achieve the same accuracy as higher-order time stepping, similar to that obtained with Fourier transforms. The same holds for an inhomogeneous problem with topography where the solution on a very fine mesh is used as reference. The series approach costs less than dispersion correction with the Fourier method and can be used on the fly during the time stepping. It does, however, require a wavelet that is sufficiently many times differentiable in time.

@en

Finite elements with polynomial basis functions on the simplex with a symmetric distribution of nodes should have a unique polynomial representation. Unisolvence not only requires that the number of nodes equals the number of independent polynomials spanning a polynomial space of a given degree, but also that the Vandermonde matrix controlling their mapping to the Lagrange interpolating polynomials can be inverted. Here, a necessary condition for unisolvence is presented for polynomial spaces that have non-decreasing degrees when going from the edges and the various faces to the interior of the simplex. It leads to a proof of a conjecture on a necessary condition for unisolvence, requiring the node pattern to be the same as that of the regular simplex.

@en

Intersecting faults are often ignored in the geomechanical simulation of induced seismicity. To investigate the effects of fault intersection and the resulting reservoir geometry on induced seismicity, caused, for instance, by gas extraction, we have developed 3D geomechanical models considering two intersecting normal faults and the surrounding horst structure. We simulate the stress field and the dynamic fault reactivation in a uniformly depleted reservoir. We observe that a smaller intersection angle increases the incremental Coulomb stress at the lower reservoir juxtaposition, thus changing the temporal rupture pattern of the seismic event. In our dynamic simulation, the rupture propagates from the main fault to the secondary fault. We conclude that the fault intersection has important effects on the induced seismicity and should be taken into account when evaluating the seismicity risk in a specific region.

@en

When solving the wave equation with finite elements, mass lumping allows for explicit time stepping, avoiding the cost of a lower-upper decomposition of the large sparse mass matrix. Mass lumping on the reference element amounts to numerical quadrature. The weights should be positive for stable time stepping and preserve numerical accuracy. The standard triangular polynomial elements, except for the linear element, do not have these properties. Accuracy can be preserved by augmenting them with higher-degree polynomials in the interior. This leaves the search for elements with positive weights, which were found up to degree 9 by various authors. The classic accuracy condition, however, is too restrictive. A sharper, less restrictive condition recently led to new mass-lumped tetrahedral elements up to degree 4. Compared to the known ones up to degree 3, they have less nodes and are computationally more efficient. The same criterion is applied here to the construction of triangular elements. For degrees 2 to 4, these turn out to be identical to the known ones. For degree 5, the number of nodes is the same as for the known element, but now there are infinitely many solutions. Some of these have a considerably larger stability limit for time stepping. For degree 6, two elements are found with less nodes than the known ones. For degree 7, one element with less nodes was found but with a negative weight, making it useless for time stepping with the wave equation. If the number of nodes is the same as for the classic element, there are now infinitely many solutions. Numerical tests for a homogeneous wave-propagation problem with a point source confirm the expected accuracy of the new elements. Some of them require less compute time than those obtained with the more restrictive accuracy criterion.

@en

Dynamic geomechanical modeling can generate the seismic wavefield caused by a fault rupture. In dynamic fault-rupture modeling, the source is considered to be finite, with a limited extent both in space and in time. This contrasts with the definition of a point source, which is generally assumed to explain the seismic wavefield caused by an earthquake. Most earlier seismic inversion studies, including those of the induced earthquakes caused by depletion of the Groningen gas field, were performed assuming a point source. Still, finding a point-source reference from the seismic wavefield, even when generated by finite faulting, is important in order to calibrate the geomechanical simulation with field-seismic observations. To this end, we have developed a workflow that links geomechanical forward modeling to seismic moment-tensor inversion. We have tested this workflow for the dynamic rupture considering a realistic 3D layered earth model. At first, we simulate the triggering of dynamic fault slip at the center of a fault plane. Next, we invert the seismograms recorded by receivers located on or near the surface to obtain the full moment-tensor point-source representation and the location of the earthquake. The results of inversion show similar waveforms for both the point source and the finite source. The location of the inverted point source is within 400 m from the center of the slip patch. The double-couple components of the inverted moment tensor also match with the strike and the dip of the fault plane.@en

Geomechanical modelling is generally used to simulate the nucleation of induce d earthquakes in, for instance the Groningen gas field. We apply quasi static simulation to investigate the stress changes from gas production. When a fault reaches a critical state, dynamic simulation provides information on the dynamic rupture during ea rthqu ake nucleation and the resulting wavefield . With the use of geomechanical modelling, it is possible to investigate the effects of the model parameters, e.g., depletion pattern and friction parameters. I n the modelling, the dynamic rupture at a finite fault is simulated both in space and time. The generated seismic wavefield from such a finite source is considered to be more realistic than the resulting wavefield from a point source. T he latter is often assumed in previous studies on the inversion of in duced earthquake data in the Groningen area. To link the wavefield generated by a geomechanically simulated finite source to the field seismic data for an earlier earthquake, we apply the same full moment tensor inversion to the waveform of a finite and of a point source . The inverted moment tensor from the field seismic observation provides a constraint to our geomechanical simulation. This allows us to perform a more realistic simulation of an induced earthquake.@en

Finite elements with mass lumping allow for explicit time stepping when modelling wave propagation and can be more efficient than finite differences in complex geological settings. In 2D on quadrilaterals, spectral elements are the obvious choice. Triangles are more flexible for meshing, but the construction of polynomial elements is less straightforward. So far, elements up to degree 9 have been found. Some years ago, an accuracy criterion that is sharper and less restrictive than the customary one led to new tetrahedral elements that are considerably more efficient than those previously known. Applying the same criterion to triangular elements provides infinitely many new elements of degree 5, with the same number of nodes as the old one, and two elements of degree 6 with less nodes than the known ones. Their efficiency, measured in terms of the compute time needed to obtain a solution with a given accuracy, is determined for a homogeneous problem and compared to that of the old elements of degree 1 to 8. For moderate accuracy, elements of degree 3 are the most efficient. For high accuracy, one of the new degree-6 elements performs best.@en

ABSTRACTWhen modeling wave propagation, truncation of the computational domain to a manageable size requires nonreflecting boundaries. To construct such a boundary condition on one side of a rectangular domain for a finite-difference discretization of the acoustic wave equation in the frequency domain, the domain is extended on that one side to infinity. Constant extrapolation in the direction perpendicular to the boundary provides the material properties in the exterior. Domain decomposition can split the enlarged domain into the original one and its exterior. Because the boundary-value problem for the latter is translation-invariant, the boundary Green functions obey a quadratic matrix equation. Selection of the solvent that corresponds to the outgoing waves provides the input for the remaining problem in the interior. The result is a numerically exact nonreflecting boundary condition on one side of the domain. When two nonreflecting sides have a common corner, the translation invariance is lost. Treating each side independently in combination with a classic absorbing condition in the other direction restores the translation invariance and enables application of the method at the expense of numerical exactness. Solving the quadratic matrix equation with Newton?s method turns out to be more costly than solving the Helmholtz equation and may select unwanted incoming waves. A proposed direct method has a much lower cost and selects the correct branch. A test on a 2D acoustic marine seismic problem with a free surface at the top, a classic second-order Higdon condition at the bottom, and numerically exact boundaries at the two lateral sides demonstrates the capability of the method. Numerically exact boundaries on each side, each computed independently with a free-surface or Higdon condition, provide even better results.@en

Induced seismicity from a gas-producing region such as Groningen is believed to be caused by reservoir depletion due to long-term gas production. However, because of the complexity and uncertainty regarding the underground structure and composition, it is difficult to quantify the effect on induced seismicity due to gas production. Here we use finite-element modelling to investigate the seismogenic potential of a pre-existing fault reactivated due to fluid depletion, considering different model settings. By applying quasi-static poroelastic loading representing reservoir depletion, the stress and strain fields are derived from the resulting displacement field. The equilibrium of the fault is then evaluated using either rate-and-state or slip-weakening behaviour for friction. When the critical state is reached on the fault, where the shear stress is greater than the friction, the reactivation of the fault takes place. This reactivation is simulated by using a dynamic solver to observe the propagation and the arrest of the dynamic faulting, as well as the resultant wavefield due to seismic slip. By comparing the depletion value at both aseismic and seismic ruptures, and looking at the stress distribution on the fault, the pattern of rupture nucleation, and the resulting seismic wavefield, we are able to evaluate separately the effect of different model settings, including the geometry and material property of both caprock and reservoir, reservoir depletion pattern, and the friction law. Furthermore, by combining our study with the observed seismic wavefield, it is possible to obtain useful insights on the spatial variation in the source region.@en