Composite memory variable velocity-stress viscoelastic modelling

Summary

The full 3-D viscoelastic wave equations have been reformulated,using the velocity-stress formulation for a curved grid, to include a set of memory variables for each particle velocity component. The earlier standard formulation includes a set of memory variables for each stress component. In 3-D this reduces the number of memory variables by three for each extra relaxation parameter above 1. This reformulation requires transforming the usual first order differential equation for velocity-stress modelling to a second order differential equation which again requires storage of two consecutive time step values of the particle velocities. No memory saving therefore is achieved using just one relaxation mechanism. However, the memory saved by the new formulation increases with increasing number of included relaxation mechanisms, and is around 30 per cent for 5–7 relaxation mechanisms, which is often required for accurate Q modelling. Incidentally, for media of heterogeneous Q, the new formulation is two to three times slower in simulations than the old formulation, so the new formulation should be used in memory critical applications and/or for heterogeneous media with layers/blocks of homogeneous Q, for which also CPU will be saved.

finite difference methods, Rayleigh waves, seismic wave propagation, topography, viscoelasticity, wave equation

1 Introduction

Full viscoelastic wave modelling, i.e. where physical attenuation and wave dispersion is included in the modelling, is important for applications as diverse as CTBT monitoring, fundamental earth structure mapping and exploration. Because of the relatively shorter wavelength/high frequency content in field exploration, viscoelasticity is more important here than for applications in earthquake modelling and teleseismic explosions (Hestholm & Ruud 2000; Hestholm 1999). However, accurate viscoelastic wave modelling was computationally unviable for a long time, until Carcione et al. (1988a, 1988b) and Carcione (1993) transformed the time convolution involved in its constitutive relation to a set of first order differential equations which could be solved computationally by numerical methods. Robertsson et al. (1994) extended this procedure from the displacement-stress formulation to the velocity-stress formulation, and Blanch et al. (1995) employed this formulation to develop procedures for modelling a desired Q. Full 3-D viscoelastic wave equations for a curved grid were then given in Hestholm (1999) in order to conform the velocity-stress formulation to free surface topography boundary conditions.

A successful method of solving the full viscoelastic wave equation with memory saving was developed by Xu & McMechan (1995, 1998). They reformulated the viscoelastic wave equations in the displacement-stress formulation to yield a set of memory variables for each displacement component as opposed to a set for each stress component in previous formulations. This exact procedure is applied in the present work to yield a set of memory variables for each velocity component in the velocity-stress formulation for a curved grid, to conform to implementation of our exact boundary conditions for free surface topography.

In the following, the standard form of the velocity-stress viscoelastic wave equations in a curved grid is given, containing an arbitrary number of relaxation mechanisms. I show their reformulation to include fewer memory variables, i.e. one set per velocity component. Then numerical aspects of the resulting new wave equations are treated in a separate section. A numerical example using the formulation is given next, where it is applied to a large domain of real surface topography in South-Western Norway. I give comparisons between the old and new formulations with respect to results, total memory requirement, and computational cost.

2 3-D Viscoelastic Wave Equations and their Reformulation

In the following, ρ is the density, π is the relaxation modulus for P waves, π = λ + 2μ where λ and μ are the Lamé parameters and μ is the relaxation modulus for S waves. τ^P and τ^S are the strain relaxation times for P and S waves respectively, and τ_σ is the stress relaxation time. The same τ_σ can be used both for P and S waves (Blanch et al. 1995). f_x, f_y and f_z are the components of the body forces, u, v and w are the particle velocity components and σ_xx, σ_yy, σ_zz, σ_xy, σ_xz and σ_yz are the stress components. r_xx, r_yy, r_zz, r_xy, r_xz and r_yz are the components of the memory variables. z₀(ξ, κ) is the function for an arbitrary free surface topography given in a rectangular, computational (ξ, κ, η)-grid of spatial coordinates, bounded by ξ = 0, ξ = ξ_max, κ = 0, κ = κ_max, η = 0 and η = η_max. For curved grid equations I define the following parameters, depending on the topography function and its spatial derivatives (Hestholm 1999),

 

(1)

 

(2)

 

(3)

A grid which is curved by vertical undulations in the two horizontal directions and has straight lines in the vertical direction (Hestholm 1999) is introduced. The extent of curvature is proportional to the distance from the bottom level of the grid, which is completely plane. The surface level of the grid coincides with the surface topography of the model. Then the velocity-stress formulation of the equations of motion, Hooke's law (the constitutive equations) and the equations for the memory variables together are the equations governing wave propagation in a linear isotropic viscoelastic medium. Hestholm (1999) gives them in a rectangular (ξ, κ, η)-grid, however in the following they are extended to an arbitrary number of L relaxation mechanisms (standard linear solids, SLSs) for more accurate modelling of a desired Q-behaviour,

 

(4)

 

(5)

 

(6)

 

(7)

 

(8)

 

(9)

 

(10)

 

(11)

 

(12)

 

(13)

 

(14)

 

(15)

 

(16)

 

(17)

 

(18)

Equations (4)–(18) are the momentum conservation equations, Hooke's law and the memory variable equations conformed to a medium bounded above by any free surface topography, given in the rectangular (ξ, κ, η)-grid.

I differentiate the momentum conservation equations (4)–(6) with respect to time and substitute for the time differentiated stresses from the stress-strain relations (7)–(12) with L relaxation mechanisms for each stress. Then we obtain

 

(19)

 

(20)

 

(21)

where f˙_x, f˙_y and f˙_z denote the time differentiated volume forces. The composite memory variable wave equations then become

 

(22)

 

(23)

 

(24)

with F_x = f˙_x, F_y = f˙_z and F_z = f˙_z being the new volume forces and

 

(25)

 

(26)

 

(27)

 

(28)

 

(29)

 

(30)

 

(31)

Each set of new memory variables m_ℓ, n_ℓ and o_ℓ is identified with each particle velocity component, u, v and w respectively. These new memory variables are given by

 

(32)

 

(33)

 

(34)

Differentiating the new memory variables m_ℓ, n_ℓ and o_ℓ with respect to t and substituting from the corresponding expressions for the old memory variables (13)–(18), gives the following first order equations for the new memory variables,

 

(35)

 

(36)

 

(37)

where the definitions (25)–(31) are used again. This set of equations assumes spatially (and temporally) constant values of all τ_σℓ, which is no practical restriction because spatial variability of τ^P_εℓ and τ^S_εℓ can account for spatially varying Q. I choose preselected, spatially constant values of 1/τ_σℓ = ω_ℓ, i.e. one specific frequency for each of the L mechanisms (Emmerich 1992; Xu & McMechan 1995).

Eqs (22)–(31) and (35)–(37) together constitute the new wave equations using composite memory variables for wave propagation in the interior of the medium. A set of memory variables is now associated with each particle velocity component instead of with each stress component as before. Exact boundary conditions for an arbitrary free surface topography are employed in this work and can be written as (Hestholm 1999; Hestholm & Ruud in preparation)

 

(38)

 

(39)

 

(40)

using definition (3) and ζ = λ/(λ + 2μ), d = [∂z₀(ξ, κ)/∂ξ], e = cos[arctan (d)], f = sin[arctan (d)], and p = [∂z₀(ξ, κ)/∂κ]e.

3 Numerical Aspects

For spatial differentiation high-order, cost-optimized finite-difference (F-D) operators were used (Kindelan et al. 1990; Holberg 1987). Staggered discretization stencils are used (Levander 1988; Virieux 1986), enabling (and requiring) us to define only the particle velocities and not the sums (25)–(31) at the surface topography. The 3-D boundary conditions (38)–(40) are discretized by second order, staggered F-D operators. Higher-order F-D methods will not decrease the numerical dispersion of Rg (fundamental mode Rayleigh) waves, only closer spatial sampling will (Xu et al. 1999), so second-order F-Ds may as well be used along the free surface topography. A direct solution for the particle velocities from second order F-D discretization of the system (38)–(40) is performed, and the resulting system of equations is unconditionally stable (Hestholm 1999; Hestholm & Ruud 2001). Below the free surface, the F-D order is gradually increased via 4 and 6 up to 8, which is the order used in the interior of the medium (Kindelan et al. 1990).

To time propagate the particle velocities from the second order differential equations (22)–(24) one may use an explicit central difference method, involving an addition of order Δt² in each time step (Δt being the time step). For small Δt though, this may lead to unfavourable rounding errors. Therefore, I apply the summed form of the method instead (Dahlquist & Bjørk 1974). With i the time step index, I define z(i + 1/2) = [u(i + 1)−u(i)]/Δt, where u is the discretized version of the full vector of particle velocities u˜(ξ, κ, η, t) = (u, v, w). z is the discretized version of a three component dummy vector z˜(ξ, κ, η, t) that has to be maintained over the entire computational domain at each time step. f˜(i) I define to be the ith time step of the discretized version of the three component vector f consisting of the right hand side of eqs (22)–(24), according to ∂²u˜/∂t² = f˜. Then we have from the central difference method

 

(41)

with n being the current time step index. I apply the summed form as follows

 

(42)

to time propagate the particle velocities. The memory variable equations (35)–(37) may become close to instability in time for small τ_σℓ compared to Δt, therefore the unconditionally stable Crank-Nicholson method was used to propagate the memory variables in time (Robertsson et al. 1994; Hestholm 1999). Since the memory variable equations are first order ordinary differential equations, the usually implicit Crank-Nicholson scheme becomes explicit, and only marginally more expensive than conventional explicit schemes.

Instead of the six stress components σ_xx, σ_yy, σ_zz, σ_xy, σ_xz and σ_yz in the old formulation we have, in the new formulation, the seven sums of eqs (25)–(31) for which we have to maintain full arrays over the complete computational domain. Both the stresses and the sums have to be kept at one time step only. The extra array included by the seventh sum is outweighed by the gain achieved by the spatially homogeneous τ_σℓ. The introduction of the second order eqs (22)–(24) requires us to store two consecutive time steps of the particle velocities and hence three full arrays extra compared to before (the vector z in our case). For one relaxation parameter (L = 1) this will be outweighed by the gain of three fewer memory variables required. As a result there will be the same memory requirement when using only one relaxation parameter. However, for each relaxation mechanism more than 1, there will be a gain of three fewer memory variables needed.

For the displacement-stress formulation (Xu & McMechan 1995) there was already an immediate memory saving when using one relaxation mechanism. This standard formulation used seven memory variables when including one relaxation mechanism, and so the immediate gain was from these seven to three memory variables by using the new formulation. Also, the displacement-stress formulation already in its standard form included second order differential equations with respect to time. Therefore the new formulation did not involve three extra full-size variables (the vector z in our case). Hence the total memory saving for the displacement-stress formulation would be larger than for the velocity-stress formulation, also because we would have a gain of four rather than three fewer memory variables required per extra relaxation parameter included. It is important, however, to emphasize that this is not because of less total memory required by the new displacement-stress formulation compared with the new velocity-stress formulation, but rather because of the pre-eminence of the standard velocity-stress formulation in already requiring less memory than the standard displacement-stress formulation.

I used MPI- (Message Passing Interface) parallelized codes by domain decomposition of the standard and new forms of velocity-stress viscoelastic wave modelling on a curved grid. Depending on the number of processors I found a memory reduction between 12.5 and 32 per cent of the new formulation for a model of 500 × 500 × 100 grid points. For this model, I found the total memory saving to be about 12.6 per cent for two relaxation parameters, 19.7 per cent for three, 24.2 per cent for four, 27.5 per cent for five, 30 per cent for six and 31.5 per cent for seven relaxation parameters.

Unfortunately, I cannot report similarly positive results for CPU times. I found the new form of the velocity-stress wave equations to be between two and three times slower than the old form, on average. This is easy to explain, however. In contrast with the new form of the displacement-stress formulation, the new velocity-stress formulation does not contain the same number of computations as before. In the two new sets of equations, the momentum conservation equations (22)–(24), and the memory variable equations (35)–(37), Θ₁–Θ₇ have to be multiplied by different sets of coefficients before spatial partial differentiation. This means that the partial differentiations have to be done separately for each of these new sets of equations. In addition, after the first and before the second differentiation we have to redefine the variables to be differentiated by multiplying them by the new appropriate coefficients. All of this should justify an approximate factor of 2.5 for the extra CPU time compared to before. As in the case of memory requirement, this high factor of extra CPU time is not because the new velocity-stress formulation requires more CPU than the new displacement-stress formulation (they require about equal CPU), but rather because the standard velocity-stress formulation requires fewer computations than the standard displacement-stress formulation. The new and standard displacement-stress formulations require the same amount of computations. Because of the added CPU time of the new velocity-stress formulation, it should be used mainly for applications where the requirement on CPU is less restrictive than on memory.

However, there are applications where the new formulation could be used to advantage. These include cases of viscoelastic homogeneous media or viscoelastic homogenous blocky media. In each such block, the coefficients in front of Θ₁–Θ₇ (the ones that depend only on the viscoelastic relaxation times) may be moved outside of the differentiation operator preceding them. Then the sums to be differentiated will be the same for both the momentum conservation equations (22)–(24) and the memory variable equations (35)–(37), and hence also no intermediate redefinition of the sums have to be performed. This leads to the same number of spatial partial differentiations and internal medium computations as for the standard velocity-stress formulation, and therefore the total CPU time will also be slightly reduced as a result of the new formulation's smaller memory requirement. I found this to be the case in numerical tests assuming viscoelastic homogeneity, for a grid of 500 × 500 × 120 points over 20 parallel processors and using one to five relaxation parameters. It is important to note that this advantageous form of the new equations may be used also in heterogenous media, as long as the medium consists of blocks or domains each of which has a homogeneous viscoelastic definition. Jumps of Q between blocks are then smoothed implicitly. In these cases we may save both a significant amount of memory and a slight amount of CPU by the new formulation.

4 Accurate Q-modelling near Fjord Topography in South-Western Norway

In this section I present an application of the new memory saving version of velocity-stress viscoelastic wave modelling, using a homogeneous medium of v_P = 6000 m/s, v_S = 3464 m/s, and density = 2000 kg/m³. The total domain is 100 × 100 × 24 km³, gridded by 500 × 500 × 120 points, so the spatial grid distance is 200 m. It is covered by an area of 100 × 100 km² of digitized elevation data from fjord-topography in South-Western Norway. Fig. (1) shows elevation (a) and gradients (b) from this area. In the simulation an incident dipping plane wave is approaching this topography from below. The plane wave consists of a plane of Ricker wavelets, each with a central frequency of 2.5 Hz, and it has the same dip of 3.6° with both the negative x and the positive y-axes, i.e. with both the westward and northward directions. It simulates a teleseismic explosion approaching from the North-West. I applied values of relaxation times from Table 5 of Blanch et al. (1995), according to an accurate constant Q modelling of Q≈20 by 5 relaxation mechanisms. Two snapshots of this simulation are shown in Fig. (2). The corresponding snapshots of the same simulation using the standard velocity-stress formulation are shown in Fig. (3).

In both simulations we confirm the strongest scattered waves to be generated from the strongest varying topography and/or steepest gradients. This can be confirmed by comparing the snapshots with Fig. (1). Since the simulations are performed for a homogeneous medium, all the wave scattering shown is due to the surface topography. The results from using the new formulation (Fig. 2) seem to contain a spectrum of far lower frequencies than those from using the standard formulation (Fig. 3). The reason is that the equations of the new formulation (22)–(24) and (35)–(37) are time-differentiated compared to the equations of the standard formulation (4)–(18). They will therefore contain relatively more high-frequency energy in their wave propagation than will the standard formulation (because their amplitudes are multiplied by the factor iω in the frequency domain). This relative high-frequency wave propagation of the new formulation will then become significantly more attenuated and dispersed with time than the corresponding propagation of the old formulation, by the homogeneous viscoelastic medium of the low Q of 20. When this attenuated field hits the surface topography, it has far less scattering efficiency than when using the standard formulation. This is seen when comparing Fig. 2 with Fig. 3. Hence we have confirmed different attributes of the new and old formulations in exhibiting different wave energy content of the same source signal.

Figure 1

a) Topography of a 100 × 100 km fjord topography area used in 3-D viscoelastic simulations. Labels are in kilometres. Blue is sea and fjords. (b) Gradients of the area.

Figure 2

Snapshots of the vertical particle velocity component w at two times after a plane wave is released near the surface topography of Fig. (1). The snapshots are taken along the surface topography and the new viscoelastic velocity-stress formulation for a curved grid is used. A constant Q≈20 is accurately modelled using five relaxation mechanisms.

Figure 3

Snapshots of the vertical particle velocity component w at two times after a plane wave is released near the surface topography of Fig. (1). The snapshots are taken along the surface topography and the standard viscoelastic velocity-stress formulation for a curved grid is used. A constant Q≈20 is accurately modelled using five relaxation mechanisms.

The left part of Fig. (4) is the west-east oriented topography profile taken from the middle position (y = 50 km) of the map in Fig. (1), whereas the right part of Fig. (4) is the south-north oriented topography profile taken from the middle position (x = 50 km) of the area. Seismograms for the profiles are shown in Figs (5)–(8). In each topography profile 50 receivers are covering the distance of 26–75 km and are interspaced by 1 km. Fig. (5) shows the west-east profile from the simulation of the new viscoelastic formulation (corresponding to Fig. 2), and confirms the relative low-frequent appearance and attenuated waves when comparing with Fig. (6). This is the corresponding west-east profile from the simulation of Fig. 3, i.e. using the standard formulation. In both cases we see that the prominent peaks at 52, 58 and 65 km of the left profile of Fig. (4) (approximately receivers 27, 33 and 40), are giving rise to prominently scattered waves.

The broad peak around 52 km from the right profile of Fig. (4) is seen to cause the most scattering in the seismograms covering this profile in Figs (7) (the new formulation) and (8) (the old formulation). The previously confirmed feature of the new formulation of a more attenuated and dispersed wavefield is similarly confirmed in these seismograms.

Figure 4

Vertically exaggerated (10 times) topographic west-east profile (left) and south-north profile (right) each at the middle of the domain covered by the surface topography of Fig. (1). Horizontal axes are in kilometres and vertical axes are in 100m.

Figure 5

Seismogram of the first horizontal particle velocity u for the simulation of Fig. (2). The receivers are located along the surface of the west-east oriented topographic profile of Fig. (4).

Figure 6

Seismogram of the first horizontal particle velocity u for the simulation of Fig. (3). The receivers are located along the surface of the west-east oriented topographic profile of Fig. (4).

Figure 7

Seismogram of the first horizontal particle velocity u for the simulation of Fig. (2). The receivers are located along the surface of the south-north oriented topographic profile of Fig. (4).

Figure 8

Seismogram of the first horizontal particle velocity u for the simulation of Fig. (3). The receivers are located along the surface of the south-north oriented topographic profile of Fig. (4).

The grid boundary reflections of the simulation of the new formulation (Fig. 2) are stronger and faster than for the standard formulation (Fig. 3). The reason is that I have only used the exponential damping procedure of Cerjan et al. (1985) in the present implementation of the new formulation. The best procedure I found for viscoelastic absorbing boundaries (Hestholm 1999) depends on explicitly setting the stress relaxation time to a high value ∼100 along the grid boundaries and changing each relaxation time linearly or by a cosine taper from its value at the boundary to its value in the interior. This procedure, however, is not feasible by the spatially constant stress relaxation times (one per modelled frequency) of the new memory saving formulation. This outlined algorithm for viscoelastic boundary absorption is however applied in the implementation of the standard formulation, using a linear change of relaxation times over 20 grid points and 10 grid points closest to the grid boundaries having constant values corresponding to a low Q. On top of this exponential damping (Cerjan et al. 1985) is applied at the 20 grid points closest to the boundaries. These 20 grid points of the method of Cerjan et al. (1985) were the only absorbing boundaries applied in the present implementation of the new formulation.

The simulation of the standard formulation used 8.6 GB of run time memory over 20 processors, whereas the simulation of the new formulation used 6.25 GB over equally many processors, representing a memory saving of 27.4 per cent. The total run time for the new formulation was about 8.5 hours, however, whereas for the standard formulation it was about 3 hr, i.e. viscoelastic homogeneity was not assumed. The Cray (SGI) Origin 2000 parallel machine at the University of Bergen, Dept. of Informatics, was used, which has a total capacity of 24 GB over 128 processors.

5 Conclusions

The full 3-D viscoelastic wave equations in the velocity-stress formulation is reformulated to include one set of memory variables per particle velocity instead of one set per stress. Almost 30 per cent memory is saved in the new formulation for accurate modelling of Q typically used in applications, although it is 2–3 times slower than the old formulation for fully heterogeneous media. The new formulation is therefore recommended for memory critical applications and/or for modelling of heterogeneous blocky media, where each block has a homogeneous Q. For the latter type of media there is a slight saving in CPU time together with the significant general saving in memory. Hence the new formulation should be used for such media. Results from applying the new formulation in combination with our exact boundary conditions for free surface topography are shown for a homogeneous medium of constant Q = 20 covered by a large area of fjord topography from South-Western Norway. I compare the results with a corresponding simulation using the old formulation. Areas of the strongest topographic variations and steepest gradients are clearly causing the strongest scattering in simulations, which are parallelized to accommodate the large run time memory even in the new memory saving version.

6 Acknowledgments

I appreciated discussions with and support from Prof. George McMechan (University of Texas at Dallas) and Prof. Eystein Husebye (University of Bergen, Norway). Discussions with Dr. Bent Ruud (University of Bergen, Norway) are greatly appreciated. I would like to thank two anonymous reviewers for constructive reviews. This research was supported by the Norwegian Research Council, the ERL/MIT Industry Consortium, the Cold Regions Research and Engineering Laboratory, US Army Corps of Engineers, under contract DACA89-99-C-0002, the Norwegian Defense Research Establishment, under contract FFI-0070, and by the Norwegian Supercomputer Committee through a grant of computing time.