# Non-periodic homogenization of 3-D elastic media for the seismic wave equation

Non-periodic homogenization of 3-D elastic media for the seismic wave equation Summary Because seismic waves have a limited frequency spectrum, the velocity structure of the Earth that can be extracted from seismic records has a limited resolution. As a consequence, one obtains smooth images from waveform inversion, although the Earth holds discontinuities and small scales of various natures. Within the last decade, the non-periodic homogenization method shed light on how seismic waves interact with small geological heterogeneities and ‘see’ upscaled properties. This theory enables us to compute long-wave equivalent density and elastic coefficients of any media, with no constraint on the size, the shape and the contrast of the heterogeneities. In particular, the homogenization leads to the apparent, structure-induced anisotropy. In this paper, we implement this method in 3-D and show 3-D tests for the very first time. The non-periodic homogenization relies on an asymptotic expansion of the displacement and the stress involved in the elastic wave equation. Limiting ourselves to the order 0, we show that the practical computation of an upscaled elastic tensor basically requires (i) to solve an elastostatic problem and (ii) to low-pass filter the strain and the stress associated with the obtained solution. The elastostatic problem consists in finding the displacements due to local unit strains acting in all directions within the medium to upscale. This is solved using a parallel, highly optimized finite-element code. As for the filtering, we rely on the finite-element quadrature to perform the convolution in the space domain. We end up with an efficient numerical tool that we apply on various 3-D models to test the accuracy and the benefit of the homogenization. In the case of a finely layered model, our method agrees with results derived from Backus. In a more challenging model composed by a million of small cubes, waveforms computed in the homogenized medium fit reference waveforms very well. Both direct phases and complex diffracted waves are accurately retrieved in the upscaled model, although it is smooth. Finally, our upscaling method is applied to a realistic geological model. The obtained homogenized medium holds structure-induced anisotropy. Moreover, full seismic wavefields in this medium can be simulated with a coarse mesh (no matter what the numerical solver is), which significantly reduces computation costs usually associated with discontinuities and small heterogeneities. These three tests show that the non-periodic homogenization is both accurate and tractable in large 3-D cases, which opens the path to the correct account of the effect of small-scale features on seismic wave propagation for various applications and to a deeper understanding of the apparent anisotropy. Numerical solutions, Computational seismology, Seismic anisotropy, Theoretical seismology, Wave propagation, Wave scattering and diffraction 1 INTRODUCTION Seismic waves are one of the most powerful tools to image the Earth’s interior. Giving access to the geometry of geological structures and to the distribution of mechanical properties within our planet, they lead to a better understanding of geodynamic processes and resource potentials. In the last decades, the seismic tomography and imaging community took advantage of the increasing computational power and the development of efficient numerical methods to improve its techniques and results. Accurate solutions to the forward modelling (see Virieux et al.2011, for a review of the various numerical methods available to model the seismic wave propagation) and the inverse problem (e.g. Pratt et al.1998; Tromp et al.2005; Plessix 2006; Fichtner et al.2008; Métivier et al.2013; Brossier et al.2015; Warner & Guasch 2016) now allow accounting for full seismic waveforms and getting robust and well-resolved models of the Earth at different scales (from the subsurface in exploration geophysics to the entire globe in seismology). One of the important remaining challenges in seismic wave simulation and inversion is the understanding of and the correct account for the effect of small heterogeneities on wave propagation. By small heterogeneities we here mean structures which are smaller than the minimum seismic wavelength propagating in the medium. Seismic waves indeed have a finite frequency band which implies the existence of a minimum wavelength, whereas heterogeneities within the Earth can occur at all scales. When propagating through small heterogeneities, seismic waves naturally average the elastic properties of the medium. A deeper understanding of this averaging process could significantly improve the interpretation of seismic inversion results, pointing out what small-scale features can be ‘hidden’ behind the smooth images one usually obtains. A proper knowledge of the averaging is of major interest for the forward modelling as well, because replacing a given discontinuous and highly heterogeneous model by an effective (or, equivalently, ‘upscaled’ or ‘long-wave equivalent’) medium greatly eases the numerical simulation of wave propagation. When present, small-scale features indeed control the spatial sampling of the model and, consequently, the time sampling too. This is because of a stability condition that all the wave simulators have to ensure (Courant et al.1928):   $$\Delta t \le C \left[ \frac{\Delta x}{V_P} \right]_{\text{min}},$$ (1)Δt denoting the time-step of the simulation, C a constant and $$\left[ \frac{\Delta x}{V_P} \right]_{\text{min}}$$ the minimum ratio of grid spacing Δx and P-wave speed VP. Small heterogeneities therefore induce massive, possibly prohibitive, computation costs. In addition, an accurate simulation usually requires all the physical discontinuities to be honoured by the mesh, which can lead to enormous meshing efforts, especially when dealing with hexahedra (e.g. Casarotti et al.2008; Peter et al.2011). For these reasons, working with effective media is much preferable. Theoretical studies for understanding the effective properties of heterogeneous elastic media was initiated in the sixties by Hashin & Shtrikman (1963) and Hill (1965). These studies both use the pioneering ideas of Voigt (1889) and Reuss (1929) to derive upper and lower-bound effective constants of periodic materials. After these works, many analytical derivations followed (e.g. Kutsenko et al.2013, and the references within), handling more and more complex periodic units or improving accuracy and efficiency of the effective properties calculation. In the case of waves propagating in finely layered media, Backus (1962) derived formula for long-wave equivalent elastic coefficients which are still widely used in seismic exploration. Fruitfully, these formulae quantify the seismic anisotropy produced by fine layering, showing that a stack of isotropic layers can explain the anisotropy observed on wave measurements. Schoenberg & Muir (1989) extended the Backus theory by including any kind of anisotropy within each layer, which enables the account for several sets of fractures in the layers. Further studies on the upscaling of fractured media then came out, such as Sayers & Kachanov (1991); Mauge & Kachanov (1994); Sayers & Kachanov (1995); Schoenberg & Sayers (1995); Schoenberg & Helbig (1997); Sayers (1998); Grechka & Kachanov (2006); Grechka (2007); Carcione et al. (2012). At larger scale, particular efforts have been focused on smoothing the Earth’s crust to ease the simulation of long-period surface waves (Capdeville & Marigo 2008; Fichtner & Igel 2008; Lekić et al.2010). More recently, the understanding of the seismic anisotropy in the upper-mantle drew the attention (Fichtner et al.2013a; Wang et al.2013; Bodin et al.2015). In this context, Jordan (2015) developed an elegant effective medium theory for random media in which statistics of the local anisotropy (produced by lattice-preferred orientation of mineral grains for instance) are separated from those of ellipsoidal geometric heterogeneities (i.e. shape-preferred orientation of simple geological structures). In this paper, the upscaling method in consideration is the two-scale homogenization technique. This method emerged in the seventies from research in micromechanics for predicting the macroscopic response of composite and random materials to either static or dynamic excitations (Bensoussan et al.1978; Sanchez-Palencia 1980; Papanicolaou & Varadhan 1981). Since then, the technique has been successfully applied to many physical processes, such as heat transfer (e.g. Allaire & Habibi 2013), Stokes flow (e.g. Hornung 1997), neutronic diffusion (e.g. Allaire & Capdeboscq 2000), magnetization (e.g. Santugini-Repiquet 2007) and elastic wave propagation (e.g. Boutin & Auriault 1993; Fish & Chen 2001, 2004; Parnell & Abrahams 2008; Bacigalupo & Gambarotta 2014). In this last field, the two-scale homogenization was adapted to non-periodic media within the last decade (Capdeville & Marigo 2007; Capdeville et al.2010a,b; Guillot et al.2010; Cance & Capdeville 2015; Capdeville et al.2015), which opened the path to the upscaling of general elastic media, with no constraint on the shape and size of the heterogeneities. The method has been tested in 1-D (Capdeville & Marigo 2007; Capdeville et al.2010a) and 2-D (Capdeville et al.2010b; Guillot et al.2010; Capdeville et al.2015); the goal of the present paper is to show its efficiency (in terms of accuracy and computation speed) on 3-D cases. To solve the homogenization problem, we will rely on a finite-element method, but it is worth noting here that the algorithm proposed by Capdeville et al. (2015) could be used as well. This latter is up and running on continuous or pixel-based 3-D media. The non-periodic two-scale homogenization actually is different to what some authors called numerical homogenization (e.g. Zijl et al.2002; Grechka 2003; Gao et al.2015) or heterogeneous multi-scale methods (e.g. Engquist et al.2009; Abdulle & Grote 2011). In these classes of techniques, the effective elastic tensor is computed by solving the wave equation in the static regime (i.e. at zero frequency so that a simple elastostatic problem arises) with a set of unit stresses applied on a representative volume. Such a volume implicitly delineates the small and the large scales. It fully depends on the mesh and the numerical technique that will then be used for simulating seismic waves. Even though our non-periodic homogenization method also includes the numerical resolution of a static problem, this latter does not involve a representative volume. Our static problem is actually defined on the whole medium at once with volume forces equal to the divergence of the elastic tensor (eq. 33). Low-pass filtering the strain and the stress which result from these forces then leads to the effective properties (eq. 30). Our filter only depends on the minimum wavelength to be propagated in the medium and a precision factor, so that the scales are separated by the waves rather than imposed by a numerical method. This makes non-periodic homogenization results physically meaningful and independent from wave equation solvers. Homogenized properties actually reveal what the waves ‘see’. In particular, Capdeville et al. (2013) showed that models obtained by full-waveform inversion are homogenized models. Applying the non-periodic homogenization method in 3-D does not require further theoretical derivations than those developed in 2-D by Capdeville et al. (2010b) and Guillot et al. (2010). In a first part, the present paper recalls the ideas and concepts of the method. From this theoretical part emerge the practical issues which have to be tackled to get homogenized properties, namely (i) the resolution of an elastostatic problem and (ii) a filtering process. The second part of the article describes the implementation of a parallel finite-element method for solving the 3-D elastostatic problem and provides some details on the 3-D filtering. In a third part, the accuracy and performance of the resulting homogenization code are challenged on various elastic models: (i) a finely layered medium for which analytical expressions of the upscaled properties are available, (ii) a large and highly heterogeneous medium in which reference seismograms can be computed and (iii) a realistic geological model made of multiple folded and faulted horizons. Finally, we discuss the possible improvements of our code and the numerous perspectives opened by the 3-D homogenization method in terms of solving forward and inverse problems. 2 THE NON-PERIODIC HOMOGENIZATION THEORY IN A NUTSHELL We here attempt to give a synoptic overview of the non-periodic homogenization theory, starting from basics introduced in micromechanics for 1-D periodic materials (e.g. Bensoussan et al.1978), then extending these basics to non-periodic media (Capdeville et al.2010a), and finally moving to the general 2-D/3-D non-periodic case (Capdeville et al.2010b; Guillot et al.2010). Our goal is to provide a digest of the theory to make the non-periodic homogenization understandable by a large number of geophysicists, so we intentionally skip some technical details of the whole Capdeville et al. (2010a,b) derivation to focus on the main ideas and concepts of the method. 2.1 1-D periodic homogenization The homogenization theory relies on an ansatz for the solution of the physical problem in consideration. In elastodynamics, the displacement u(x; t) and the stress σ(x; t) involved in the 1-D case (x being the space variable and t being the time) are postulated to be two-scale asymptotic expansions:   $$u(x;t) = \sum _{i=0}^{+\infty } \varepsilon ^i u_i\left(x,\frac{x}{\varepsilon };t\right),$$ (2)  $$\sigma (x;t) = \sum _{i=0}^{+\infty } \varepsilon ^i \sigma _i\left(x,\frac{x}{\varepsilon };t\right),$$ (3)where $$\varepsilon =\frac{l}{\lambda _m}$$ is the ratio of the size of the periodic cell which constitutes the 1-D medium to the minimum wavelength propagating in this medium (Fig. 1). By definition, l is microscopic and λm is macroscopic. ε, which is called scaling parameter, is therefore much smaller than 1. It enables to explicitly separate the scales within coefficients ui and σi of series (2) and (3), x capturing the large-scale variations and $$y=\frac{x}{\varepsilon }$$ handling the small-scale variations. Because any change in y induces a very slight change in x, the two variables can be treated independently and the spatial derivative operator can be written   $$\nabla = \varepsilon ^{-1}\nabla _y + \nabla _x.$$ (4) Figure 1. View largeDownload slide 1-D periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in an infinite 1-D medium made of periodic cells whose size l is much smaller than λm. Figure 1. View largeDownload slide 1-D periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in an infinite 1-D medium made of periodic cells whose size l is much smaller than λm. The introduction of the small-scale variable $$y=\frac{x}{\varepsilon }$$ also allows us to rewrite the physical parameters E(x) (Young modulus) and ρ(x) (density) of the medium as λm-periodic quantities E(y) and ρ(y). In other words, the l-periodic bar in the initial one-variable problem is now seen as a medium containing only small scales which are repeated in space with a λm-periodicity. Furthermore, coefficients ui and σi are assumed to be λm-periodic in y as are E and ρ. This assumption imposes that small-scale variations of the displacement and stress fields are due to local small-scale structures of the medium. In seismology, such a phenomenon is commonly called site effect. Plugging eqs (2)–(4) into the elastodynamic problem (i.e. the wave equation and Hooke’s law) yields a cascade of coupled equations which can be solved for each i using the average over the periodic cell   $$\langle f \rangle (x) = \frac{1}{\lambda _m}\int _{-\frac{\lambda _m}{2}}^{\frac{\lambda _m}{2}} f(x,y) \, \text{d}y, \, \forall f \! : \! \mathbb {R}^2 \! \rightarrow \! \mathbb {R},$$ (5)and the periodicity in y of the problem. Doing so, it turns out that the zeroth-order terms u0 and σ0 do not depend on the microscopic variable y. This result formalizes the poor sensitivity of the wavefield to small heterogeneities. Going further, one shows that u0 and σ0 are the solution of the so-called homogenized problem, which is a classical elastodynamic problem with homogeneous effective properties E⋆ and ρ⋆ such that   $$\rho ^{\star } = \langle{\rho (y)}\rangle$$ (6)and   $$E^{\star } = \left\langle {E(y)[1+\nabla _y\chi (y)]}\right\rangle .$$ (7)In this last equation, χ is the so-called first-order corrector. It is the solution of the so-called cell problem, which is   $$\nabla _y\left[E(1+\nabla _y\chi )\right] = 0$$ (8)defined on the cell with periodic boundary conditions. Because an analytical solution exists for this problem, eq. (7) finally reduces to   $$E^{\star } = \left\langle {\,\frac{1}{E(y)}\,}\right\rangle ^{-1}.$$ (9)In conclusion, the quantity to average for obtaining the zeroth-order long-wave equivalent Young modulus E⋆ is the inverse of the initial Young modulus E (eq. 9), whereas the zeroth-order long-wave equivalent density ρ⋆ simply is the average of the initial density ρ (eq. 6). Contrary to u0 and σ0, the higher-order terms of series (2) and (3) contain small-scale variations. For instance, the first-order displacement term is y-dependent through the first-order corrector:   $$u_1(x,y) = \chi (y) \nabla _xu_0(x) + \langle {u_1}\rangle (x).$$ (10)Examples of the contribution of this last term to the whole wavefield can be found in Capdeville et al. (2010a). In the present work, we focus on the order 0 so we do not provide any further details on the higher-order terms. A crucial point in the homogenization theory is that the asymptotic convergence can be proved mathematically. One actually does not directly show the convergence of series (2) and (3) towards the exact wavefield; one uses the so-called Γ-convergence instead (Dal Maso 1993). Rather than studying a single problem for the physically relevant value of ε, one considers a sequence of problems indexed by ε which is now regarded as a small parameter going to zero. In other words, one builds a fictitious sequence of problems in which the periodic cell becomes smaller and smaller. In this context, one demonstrates that the exact solution converges to the solution of the homogenized problem u0 and σ0 (Nguetseng 1989; Allaire 1992). Such a demonstration provides rigorous mathematical foundations for the homogenization theory. 2.2 1-D non-periodic homogenization Let us now consider any given distribution of the properties E and ρ within a finite 1-D medium. The length of the medium is denoted by L. Again, our goal here is to find long-wave equivalent properties. We denote by λm the minimum wavelength of the wavefield propagating in the bar (Fig. 2). Inspired by the periodic homogenization, we will represent the displacement u(x; t) and the stress σ(x; t) as power series in ε (to be redefined) with coefficients ui and σi as functions of x and $$y=\frac{x}{\varepsilon }$$. Plugging these series into the wave equation and Hooke’s law will yield a cascade of equations. To get a solution for these equations, let us bring a periodicity into the problem by imposing periodic boundary conditions at the border of the domain. We end up with an infinite medium made of a periodic cell, so we can certainly benefit from the mathematical results of the previous section (eqs 6 and 9). Nevertheless, two big issues remain: Because the boundary conditions are forced to be periodic, the effective properties computed near the border of the domain are not meaningful for any kind of physically interesting conditions such as Dirichlet or Neumann conditions. In the present paper, we do not investigate this issue. Disregarding boundary effects, we focus on the computation of accurate effective properties in the interior of the domain. Precisions on what we exactly mean by the interior of the domain are given in Part 3. Contrary to the periodic case, the size of the periodic cell is not microscopic. The cell can contain various sizes of heterogeneity, including macroscopic scales (i.e. sizes equal to or larger than λm), so ε defined as the ratio of the size L of the periodic cell to the minimum wavelength λm no longer is a relevant scaling parameter to separate the scales through variables x and $$y=\frac{x}{\varepsilon }$$. The goal of the present section is (i) to redefine ε and (ii) to separate the scales within the mechanical properties of the bar in a way that allows for a solution of our homogenization problem. Figure 2. View largeDownload slide 1-D non-periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in a finite 1-D medium of size L. Contrary to the periodic case, no assumption is made on the distribution of the mechanical properties E(x) and ρ(x) within the medium. Figure 2. View largeDownload slide 1-D non-periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in a finite 1-D medium of size L. Contrary to the periodic case, no assumption is made on the distribution of the mechanical properties E(x) and ρ(x) within the medium. 2.2.1 Redefining ε Because there is no quantity for defining the small scales yet, we introduce λ0 < λm: all the heterogeneities whose size is smaller than λ0 are considered as small. λ0 therefore is analogous to l in the periodic case. Using it, we can bring in a new scaling parameter $$\varepsilon _0=\frac{\lambda _0}{\lambda _m} < 1$$ and we can choose ε ≤ ε0. Similarly to the periodic case, ε is a formal quantity which can go to zero to prove the Γ-convergence of the asymptotic solution. It can be seen as the ratio of a small length λ ≤ λ0 to the minimum wavelength of the wavefield to be propagated λm, λ becoming smaller and smaller when studying the convergence. In practice, the only physically relevant value of ε is ε0. With this definition of ε, variables x and $$y=\frac{x}{\varepsilon }$$ can be treated independently. We denote by Ψ the set of functions $$f(x,y) \! : \mathbb {R}^2 \! \rightarrow \! \mathbb {R}$$ such that the x part of f carries the large-scale variations while the y part of f handles the small-scale variations. As in the periodic case, we impose coefficients ui and σi of the asymptotic expansions (2) and (3) to belong to Ψ. Moreover, we can write the spatial derivative operator ∇ as in eq. (4). 2.2.2 Separating the scales within the 1-D medium In the non-periodic case, the main issue we face is that the medium contains both microscopic and macroscopic scales, so we cannot write E and ρ as a function of y only. Using a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$ (Appendix A) to separate large-scale and small-scale variations, we will have to find a proper way to build $$E^{\varepsilon _0}(x,y)$$ and $$\rho ^{\varepsilon _0}(x,y)$$ from E(x) and ρ(x). By proper we here mean which allows for a solution to the cascade of equations that arise when plugging eqs (2)–(4) into the elastodynamic problem. Let us introduce $$\eta =\frac{L}{\lambda }$$ (which implies that $$\varepsilon =\frac{L}{\eta \lambda _m}$$) and assume that we properly built $$E^{\varepsilon _0}(x,y)$$ and $$\rho ^{\varepsilon _0}(x,y)$$ from E(x) and ρ(x). Because E and ρ are L-periodic in x, $$E^{\varepsilon _0}$$ and $$\rho ^{\varepsilon _0}$$ are ηλm-periodic in y. Assuming that coefficients ui(x, y) and σi(x, y) are also ηλm-periodic in y and using the average over the periodic cell   $$\langle{f}\rangle (x) = \frac{1}{\eta \lambda _m}\int _{-\frac{\eta \lambda _m}{2}}^{\frac{\eta \lambda _m}{2}} \! f(x,y) \, \text{d}y , \forall f \! : \! \mathbb {R}^2 \! \rightarrow \! \mathbb {R}\, ,$$ (11)we can easily derive the cascade of equations. Similarly to the periodic case, it turns out that the zeroth-order displacement u0 and stress σ0 do not depend on y. Again, these fields are the solution of the so-called homogenized problem, which is a classical elastodynamic problem involving effective properties $$\rho ^{\varepsilon _0\star }$$ and $$E^{\varepsilon _0\star }$$ such that   $$\rho ^{\varepsilon _0\star }(x) = \langle {\rho ^{\varepsilon _0}(x,y)}\rangle$$ (12)and   $$E^{\varepsilon _0\star }(x) = \left\langle {E^{\varepsilon _0}(x,y)[1+\nabla _y\chi ^{\varepsilon _0}(x,y)]}\right\rangle .$$ (13)Note that the effective properties no longer are homogeneous in the non-periodic case. Moreover, they depend on the choice of ε0. The first-order corrector $$\chi ^{\varepsilon _0}(x,y)$$ also depends on ε0. It necessarily belongs to Ψ and it is ηλm-periodic in y. Similarly to the periodic case, it is the solution of the cell problem, that is,   $$\nabla _y\left[E^{\varepsilon _0}(1+\nabla _y\chi ^{\varepsilon _0})\right] = 0$$ (14)with periodic boundary conditions. As shown by Capdeville et al. (2010a), function $$\nabla _y\chi ^{\varepsilon _0}(x,y)$$ therefore satisfies   $$\nabla _y\chi ^{\varepsilon _0} = -1 + \left\langle {\,\frac{1}{E^{\varepsilon _0}}\,}\right\rangle ^{-1} \frac{1}{E^{\varepsilon _0}}.$$ (15)This last equation implies that Eq. (13) reduces to   $$E^{\varepsilon _0\star }(x) = \left\langle {\,\frac{1}{E^{\varepsilon _0}(x,y)}\,}\right\rangle ^{-1},$$ (16)which along with eq. (12) tells the quantity to average to get the zeroth-order effective properties. $$\frac{1}{E^{\varepsilon _0}(x,y)}$$ belongs to Ψ, which is the condition to meet for building $$E^{\varepsilon _0}(x,y)$$ properly. Another equation from our cascade tells that $$\rho ^{\varepsilon _0}(x,y)$$ has to lie in Ψ (Appendix B), so we also have a condition for building $$\rho ^{\varepsilon _0}(x,y)$$ properly. Following these conditions, we easily form   $$E^{\varepsilon _0}(x,y) \, = \left[ \mathcal {F}^{\varepsilon _0}\left\lbrace \frac{1}{E}\right\rbrace \!(x) \, + \left[\frac{1}{E}-\mathcal {F}^{\varepsilon _0}\left\lbrace \frac{1}{E}\right\rbrace \right]\!(y) \right]^{-1}$$ (17)and   $$\rho ^{\varepsilon _0}(x,y) \, = \, \mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x) \, + \, [\rho -\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace ](y).$$ (18)Fig. 3 illustrates the mathematical construction. In this figure, symbol g represents either ρ or $$\frac{1}{E}$$, and FT stands for Fourier transform. From any distribution of g over a length L extended to $$\mathbb {R}$$ by periodicity (Fig. 3a), large scales (Fig. 3b) and small scales (Fig. 3c) are extracted using $$\mathcal {F}^{\varepsilon _0}$$. Within the small scales, space variable x is replaced by $$y=\frac{x}{\varepsilon }$$ (Fig. 3d), thus changing the L-periodicity in x to a ηλm-periodicity in y. Then, the two-variable quantity $$g^{\varepsilon _0}(x,y)$$ is formed by simply adding the large scales (which are a function of x) to the small scales (which are a function of y). Fig. 3(e) represents $$g^{\varepsilon _0}$$ for a particular value of $$x= \bar{x}$$; it shows the small scales oscillating around $$\mathcal {F}^{\varepsilon _0}\lbrace g\rbrace ( \bar{x})$$. Figure 3. View largeDownload slide Separation of the scales within a 1-D medium of length L. g here represents either ρ or $$\frac{1}{E}$$. From any distribution of g extended to $$\mathbb {R}$$ by periodicity (top row), we build $$g^{\varepsilon _0}(x,y)$$ (bottom row) using a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$. In this figure, $$g^{\varepsilon _0}(x,y)$$ is represented for a given $$x= \bar{x}$$. See the text for more details. Figure 3. View largeDownload slide Separation of the scales within a 1-D medium of length L. g here represents either ρ or $$\frac{1}{E}$$. From any distribution of g extended to $$\mathbb {R}$$ by periodicity (top row), we build $$g^{\varepsilon _0}(x,y)$$ (bottom row) using a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$. In this figure, $$g^{\varepsilon _0}(x,y)$$ is represented for a given $$x= \bar{x}$$. See the text for more details. 2.2.3 Final results Including eq. (18) into (12) yields   \begin{eqnarray} \rho ^{\varepsilon _0\star }(x) & =& \left\langle {\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x) + [\rho -\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace ](y)}\right\rangle \nonumber \\ &=& \left\langle {\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x)}\right\rangle + \langle {[\rho -\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace ](y)}\rangle \nonumber \\ & =& \mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x). \end{eqnarray} (19) Similarly, including eq. (17) into (16) leads to   $$E^{\varepsilon _0\star }(x) = \left[ \mathcal {F}^{\varepsilon _0}\left\lbrace \frac{1}{E}\right\rbrace (x) \right]^{-1}.$$ (20)Eqs (19) and (20) tell that the zeroth-order effective properties of any 1-D medium can be computed by filtering the density and the inverse of the Young modulus. This result can be easily intuited from the solution inferred in the periodic case (eqs 6 and 9). We have here given a rigorous demonstration for it. To derive a zeroth-order solution for our non-periodic homogenization problem, we introduced a new scaling parameter ε0. Because the obtained effective medium depends on this parameter (eqs 19 and 20), u0 and σ0 also depend on it. For sake of simplicity, we did not index these two fields by ε0 as Capdeville et al. (2010a) did. Nevertheless, we will study the ε0-convergence of these fields in a specific case (Section 4.2). 2.3 3-D non-periodic homogenization In the 3-D case, Hooke’s law involves a fourth-order tensor  C to relate the stress and the strain, so the linear elastic behaviour of a given medium no longer is fully described by E only. Nevertheless, the homogenization theory developed for 1-D non-periodic media is still valid up to the zeroth-order effective properties (12) and (13) which are now written   $$\rho ^{\varepsilon _0\star }(\boldsymbol {x}) = \left\langle {\rho ^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})}\right\rangle _3$$ (21)and   \begin{eqnarray} {\bf C}^{\varepsilon _0\star }(\boldsymbol {x}) = \left\langle {\, {\bf C}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) : \left[ {\bf I}+ \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) \right) \right] \,}\right\rangle _3\!,\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (22)where I is the fourth-order identity tensor, t is the transpose operator, : is the tensor contraction [A:  B]ijkl = AijmnBmnkl and 〈〉3 is the average on y of any function $$f(\boldsymbol {x},\boldsymbol {y}) \! : \! \mathbb {R}^{3\times 2} \! \rightarrow \! \mathbb {R}$$ over the periodic cell. The first-order corrector $$\boldsymbol {\chi }^{\varepsilon _0}$$ now is a third-order tensor. It is periodic in y and it necessarily belongs to Ψ3 (the extension of Ψ to 3-D). It is the solution of the following cell problem:   $$\boldsymbol {\nabla }_y \cdot \left\lbrace {\bf C}^{\varepsilon _0} : \left[ {\bf I}+ \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} \right) \right] \right\rbrace = \boldsymbol {0}$$ (23)with periodic boundary conditions. Contrary to the 1-D case, there is no analytical solution for this problem here (unless the medium is layered transverse isotropic, as demonstrated by Guillot et al. (2010) and Lin et al. (2017)), so we are not able (i) to write $${\bf C}^{\varepsilon _0\star }$$ as a function of $${\bf C}^{\varepsilon _0}$$ only and (ii) to get a mathematical condition on $${\bf C}^{\varepsilon _0}$$ for properly separating the scales within the elastic properties of the medium. To overcome this issue, a procedure inspired by Papanicolaou & Varadhan (1981) is proposed. Introducing $${\bf G}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})$$ and $${\bf H}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})$$ such that   \begin{eqnarray} {\bf H}^{\varepsilon _0} = {\bf C}^{\varepsilon _0}:{\bf G}^{\varepsilon _0} \end{eqnarray} (24)  \begin{eqnarray} \phantom{{\bf H}^{\varepsilon _0}}= {\bf C}^{\varepsilon _0} : \left[{\bf I}+\frac{1}{2} (\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}+{^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0})\right], \end{eqnarray} (25)this procedure enables an implicit construction of $${\bf C}^{\varepsilon _0\star }$$ such that both $${\bf H}^{\varepsilon _0}$$ and $$\boldsymbol {\chi }^{\varepsilon _0}$$ are periodic in y and belong to Ψ3, which are the conditions to meet to get a proper solution to our homogenization problem. Here are the steps of the procedure: Eq. (23) with periodic boundary conditions is solved using C(y) instead of $${\bf C}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})$$. C(y) is just the original elastic tensor in which the space variable x has been changed by y = ε0x, meaning that all the scales within the medium are considered as small in this first step. This yields a cell problem that we can solve numerically (Part 3). The solution of such a cell problem is called the starting corrector χs(y). It can be seen as the static response of the medium to local unit strains expressed by the identity tensor I. From χs, two tensors Gs and Hs are built:   $${\bf G}_s(\boldsymbol {y}) = {\bf I}+ \frac{1}{2} \left[ \boldsymbol {\nabla }_y\boldsymbol {\chi }_s(\boldsymbol {y}) + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }_s(\boldsymbol {y}) \right],$$ (26)  $${\bf H}_s(\boldsymbol {y}) = {\bf C}(\boldsymbol {y}) : {\bf G}_s(\boldsymbol {y}).$$ (27)Gs can be seen as the unit strains plus the strains associated with the starting corrector. It is called strain concentrator. In the same way, Hs is called stress concentrator. Using $$\mathcal {F}^{\varepsilon _0}_3$$ (the extension of $$\mathcal {F}^{\varepsilon _0}$$ to 3-D, see Appendix C), the large-scale and small-scale variations are separated in Gs and Hs to form $${\bf G}^{\varepsilon _0}$$ and $${\bf H}^{\varepsilon _0}$$:   $${\bf G}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) = {\bf I}+ \left[{\bf G}_s - \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace \right ](\boldsymbol {y}) : \left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace (\boldsymbol {x}) \right]^{-1} ,$$ (28)  \begin{eqnarray} {{\bf H}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})}\nonumber\\ {\quad = \left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf H}_s\rbrace (\boldsymbol {x}) + [{\bf H}_s - \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf H}_s\rbrace ](\boldsymbol {y}) \right] :\left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace (\boldsymbol {x}) \right]^{-1}.} \end{eqnarray} (29)Through definitions (24) and (25), $${\bf G}^{\varepsilon _0}$$ and $${\bf H}^{\varepsilon _0}$$ are the strain and the stress concentrators associated with the first-order corrector $$\boldsymbol {\chi }^{\varepsilon _0}$$. They are here constructed by separating the scales in Gs and Hs (eqs 28 and 29). From such a construction, one can demonstrate that $${\bf H}^{\varepsilon _0}$$ and $$\boldsymbol {\chi }^{\varepsilon _0}$$ are indeed periodic in y and belong to Ψ3 (Appendix D). Using (25) in (22), we note that $${\bf C}^{\varepsilon _0\star } = \left\langle {{\bf H}^{\varepsilon _0}}\right\rangle_3$$. Introducing (29) into this latter equation, we end up with   $${\bf C}^{\varepsilon _0\star }(\boldsymbol {x}) = \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf H}_s\rbrace (\boldsymbol {x}) : \left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace (\boldsymbol {x}) \right]^{-1}.$$ (30) The zeroth-order effective elastic tensor given by the procedure therefore is obtained by filtering the stress and the strain associated with the starting corrector. This means that one just needs to go through steps (i), (ii) and (iv) to get this tensor in practice. Step (iii) would be necessary to obtain $$\boldsymbol {\chi }^{\varepsilon _0}$$ (through the construction of $${\bf C}^{\varepsilon _0} = {\bf H}^{\varepsilon _0} : {\bf G}^{\varepsilon _0^{\ -1}}$$ and the resolution of eq. 23 for instance) in order to access the first-order displacement. Because we focus on the zeroth-order solution in the present paper, we will not deal with this latter aspect. Furthermore, there is no demonstration for the minor and major symmetries of $${\bf C}^{\varepsilon _0\star }$$ yet. We do not address this point here, but we have checked the symmetries (up to a certain precision) in all the applications that we present in Part 4. As regards the density, an equation similar to (B1) also emerges in 3-D, implying that $$\rho ^{\varepsilon _0}$$ must lie in Ψ3, so the scales can be separated as they were in 1-D:   $$\rho ^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) \, = \, \mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace (\boldsymbol {x}) \, + \, \left[\rho -\mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace \right ](\boldsymbol {y}).$$ (31)Using (31) in (21), the effective density comes out:   $$\rho ^{\varepsilon _0\star } = \mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace.$$ (32)As in the 1-D case, $$\rho ^{\varepsilon _0\star }$$ is obtained by simply filtering the initial density ρ. 3 IMPLEMENTATION OF THE 3-D NON-PERIODIC HOMOGENIZATION While the non-periodic homogenization theory involves many specific concepts and quantities, its zeroth-order result can be stated quite shortly. From a given 3-D medium described by its density ρ and its elastic tensor C, a long-wave equivalent density $$\rho ^{\varepsilon _0\star }$$ can be computed by just filtering ρ (eq. 32), and a long-wave elastic tensor $${\bf C}^{\varepsilon _0\star }$$ can be obtained by Calculate the starting corrector χs from the cell problem (23) using C with periodic boundary conditions instead of $${\bf C}^{\varepsilon _0}$$, that is,   $$\boldsymbol {\nabla }_y \cdot \left\lbrace {\bf C}: \left[ {\bf I}+ \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }_s + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }_s \right) \right] \right\rbrace = \boldsymbol {0}.$$ (33) Build Gs and Hs from χs using eqs (26) and (27), and then filter these two quantities to get $${\bf C}^{\varepsilon _0\star }$$ (eq. 30). Differential eq. (33) can be seen as a classic elastostatic problem with a specific load consisting in the divergence of the elastic tensor C. Such a divergence yields a third-order tensor ∂iCijkl. Thanks to the symmetry of C (Cijkl = Cijlk), this tensor reduces to six force vectors. To determine the six corresponding displacements, eq. (33) can be solved in two different ways: using a strong-form iterative scheme based on Fast Fourier Transform (Moulinec & Suquet 1998; Capdeville et al.2015) or using a more classic weak-form finite-element approach (e.g. Hughes 2012). The latter is adapted to strongly discontinuous media. It is the method we implement here, relying on tetrahedral meshes. Because we want to investigate the behaviour of the solution with respect to the discretization, our code allows polynomial interpolations of degree 1, 2 or 3 (Worth et al.2012) in addition with various quadrature rules (Felippa 2004). Moreover, both linear and quadratic tetrahedra are enabled, leading to either iso-, super- or subparametric elements. The obtained linear system involving six right-hand side members, we solve it using a direct solver. Among several codes, PARDISO (Schenk & Gärtner 2006) is chosen here. Finally, the low-pass filter $$\mathcal {F}^{\varepsilon _0}_3$$ is applied in the space domain to obtain $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$. To perform the 3-D convolution, we use the quadrature associated with the mesh employed in finite element analysis. The periodicity imposed at the boundary ∂Ω of the elastic domain Ω involved in eq. (33) means that the medium is supposed to repeat itself periodically in the 3-D. When dealing with geological media, this condition is obviously not fulfilled. It is therefore replaced by either a homogeneous Neumann condition or a homogeneous Dirichlet condition. We choose the second option in our implementation so we impose χs = 0 on ∂Ω. The effect in the volume Ω of such an artificial condition decays exponentially (Dumontet 1990), so our numerical solution is meaningless in a thin layer from the border of the domain. Such meaningless values of χs do not matter, because the filtering process cannot be performed near ∂Ω anyway. Some elastic material to be convolved with the wavelet is actually missing there, so we are not able to compute the effective properties $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$ using (32) and (30). The layer ΩO in which the filter cannot be applied is called the outer domain. Its thickness is equal to half of the wavelet support. Its complement ΩI = Ω − ΩO is called the inner domain (or the interior of the domain). The solutions proposed in the present paper make sense in ΩI only. Further developments, such as those initiated by Capdeville & Marigo (2008, 2013), would be necessary to get meaningful effective properties in ΩO. When handling large models, the memory requirements for achieving the computation of the effective properties can be very large. This is mainly because solving large linear systems, even symmetric, is memory-demanding. This is also because high-order tensors are involved in the homogenization process. For these two reasons, a distributed-memory computation is necessary. A classic way of implementing it is presented in Fig. 4. It consists in Decomposing the domain Ω into n subdomains Ωk such that $$\Omega = \bigcup _{k=1}^n \Omega _k$$. Each subdomain being handled by a processor, n local stiffness matrices Kk are computed. Assembling the local stiffness matrices to build the global stiffness matrix K. Solving the linear system using a parallel solver to get the finite element solution χs of eq. (33). From χs, Gs is formed following (26). Using the domain decomposition again to filter Gs, Hs = C: Gs and ρ in each subdomain. From this filtering, the effective properties $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ are derived. In the outer domain ΩO (the dotted grey lines in Fig. 4), these properties cannot be computed because the filter cannot be applied. The same problem appears near the frontiers between the subdomains, so we have to either implement massive communications between the subdomains or enlarge each subdomain Ωk by a buffer layer $$\Omega _k^B$$ which is as thick as half of the filtering wavelet support. This second option is sketched in Fig. 4. It yields n overlapping parts$$P_k = \Omega _k \cup \Omega _k^B$$. Along with possible pieces of ΩO in Ωk, the buffer layer $$\Omega _k^B$$ acts as an outer part$$P_k^O = (\Omega ^O \cap \Omega _k) \cup \Omega _k^B$$. The effective properties can be calculated in the inner part$$P_k^I = P_k - P_k^O$$, so in the whole inner domain because $$\bigcup _{k=1}^n P_k^I = \Omega ^I$$. Communicating between processors to gather the results from the n subdomains. Figure 4. View largeDownload slide A classic distributed-memory workflow for solving the non-periodic homogenization problem. (a) The domain Ω is decomposed in n subdomains Ωk. As an example here n = 3. In each subdomain, a local stiffness matrix Kk is computed. (b) The local stiffness matrices are assembled to obtain the global stiffness matrix K. (c) The linear system is inverted using a parallel solver. From the gradient of the solution, Gs is formed. (d) Using the domain decomposition again, the effective properties $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ are computed in each subdomain by filtering Gs, Hs = C: Gs and ρ. In the outer domain (dotted grey), the filtering wavelet $$w^{\varepsilon _0}$$ cannot be applied so $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ cannot be computed. The same problem appears near the frontiers between the subdomains, so we enlarge these latter to build n overlapping parts Pk equipped with an outer part $$P_k^O$$ and an inner part $$P_k^I$$ in which the filter can be applied. (e) The result in the whole inner domain ΩI is obtained by merging the results from the n inner parts. Figure 4. View largeDownload slide A classic distributed-memory workflow for solving the non-periodic homogenization problem. (a) The domain Ω is decomposed in n subdomains Ωk. As an example here n = 3. In each subdomain, a local stiffness matrix Kk is computed. (b) The local stiffness matrices are assembled to obtain the global stiffness matrix K. (c) The linear system is inverted using a parallel solver. From the gradient of the solution, Gs is formed. (d) Using the domain decomposition again, the effective properties $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ are computed in each subdomain by filtering Gs, Hs = C: Gs and ρ. In the outer domain (dotted grey), the filtering wavelet $$w^{\varepsilon _0}$$ cannot be applied so $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ cannot be computed. The same problem appears near the frontiers between the subdomains, so we enlarge these latter to build n overlapping parts Pk equipped with an outer part $$P_k^O$$ and an inner part $$P_k^I$$ in which the filter can be applied. (e) The result in the whole inner domain ΩI is obtained by merging the results from the n inner parts. Such a workflow would certainly work, but it imposes the resolution of the whole linear system to fit the distributed-memory at once, meaning that a large stiffness matrix requires a large parallel computer. To overcome this limitation, we propose an alternative implementation which consists in working on the parts Pk all along the workflow (Fig. 5). Applying a homogeneous Dirichlet condition χs = 0 at the boundary of the parts, each of them becomes like a new domain which is totally independent from the other parts. Thanks to the outer parts $$P_k^O$$, we get the effective properties in all the inner parts $$P_k^I$$, so in the whole inner domain $$\Omega ^I = \bigcup _{k=1}^n P_k^I$$. This alternative implementation is embarrassingly parallel because the parts are independent from each other all along the workflow, which means that they can be treated sequentially. This allows for the homogenization of large models on small computers, possibly laptops, provided the available RAM covers the memory requirement of every single part. Such a distributed-memory implementation is also proposed by Capdeville et al. (2015). Figure 5. View largeDownload slide An alternative way of solving the non-periodic homogenization in parallel. (a) The domain Ω is decomposed in n parts Pk. As an example here n = 3. The parts overlap to equip them all with an outer part $$P_k^O$$ and an inner part $$P_k^I$$. Moreover, a homogeneous Dirichlet condition is applied at the boundary of the parts to make each of them like a new domain in which a linear system with a stiffness matrix Kk can be defined. (b) Solving the linear system in each part leads to n tensors Gsk. (c) Filtering Gsk, Hsk = C: Gsk and ρ, the effective properties are obtained in each inner part $$P_k^I$$. (d) The results from the n inner parts are gathered to obtain the effective properties in the whole inner domain ΩI. Figure 5. View largeDownload slide An alternative way of solving the non-periodic homogenization in parallel. (a) The domain Ω is decomposed in n parts Pk. As an example here n = 3. The parts overlap to equip them all with an outer part $$P_k^O$$ and an inner part $$P_k^I$$. Moreover, a homogeneous Dirichlet condition is applied at the boundary of the parts to make each of them like a new domain in which a linear system with a stiffness matrix Kk can be defined. (b) Solving the linear system in each part leads to n tensors Gsk. (c) Filtering Gsk, Hsk = C: Gsk and ρ, the effective properties are obtained in each inner part $$P_k^I$$. (d) The results from the n inner parts are gathered to obtain the effective properties in the whole inner domain ΩI. To speed up our code, multithreaded computations are performed whenever possible. Moreover, efficient algorithms based on k–d trees and stack data structures (e.g. Cormen et al.2009) are used to search for elements or points across the finite element mesh. Last but not least, the partitioning of the mesh is performed on ΩI instead of Ω to obtain well-balanced parts Pk (Fig. 6). The performance and capabilities of the code are illustrated through examples in the next part. Figure 6. View largeDownload slide Strategy for building well-balanced parts. (a) We partition ΩI instead of the whole domain Ω, then forming n subdomains $$\Omega ^I_k$$. As an example here n = 3. (b) Each subdomain is enlarged by a buffer layer $$\Omega ^B_k$$. At the border of the domain, the buffer is the outer domain ΩO (dotted grey). We end up with n parts $$P_k\!=\!\Omega ^I_k\cup \Omega ^B_k$$ which all have the same size. Moreover, $$P_k^I = \Omega _k^I$$ and $$P_k^O = \Omega ^B_k$$. Figure 6. View largeDownload slide Strategy for building well-balanced parts. (a) We partition ΩI instead of the whole domain Ω, then forming n subdomains $$\Omega ^I_k$$. As an example here n = 3. (b) Each subdomain is enlarged by a buffer layer $$\Omega ^B_k$$. At the border of the domain, the buffer is the outer domain ΩO (dotted grey). We end up with n parts $$P_k\!=\!\Omega ^I_k\cup \Omega ^B_k$$ which all have the same size. Moreover, $$P_k^I = \Omega _k^I$$ and $$P_k^O = \Omega ^B_k$$. 4 VALIDATION TESTS We here handle three different models to test the accuracy of the homogenization method in 3-D and challenge our code. First, the case of fine layers is investigated. For such a medium, analytical expressions of the effective elastic parameters exist so we can just compare the result of the homogenization with the result of these expressions to validate our code. The second model we study is made of small cubes. Because no reference solution for the effective properties is available in this case, we base our validation on the comparison of waveforms computed in the initial medium on the one hand, and in the homogenized medium on the other hand. Finally, a 3-D geological model is handled to emphasize the efficiency of the homogenization in a realistic case. 4.1 Homogenization of a finely layered medium The first model we consider for testing our 3-D homogenization code is a medium made of 60 isotropic layers randomly taken between 800 and 1200 m thick. Within each layer, the density is randomly chosen in the 2000–4000 kg m−3 range. Because the medium is isotropic, two parameters (e.g. the Lamé coefficients, or the S- and P-wave velocities) are sufficient to define the elastic tensor in each layer. We then randomly choose the S-wave velocity between 3000 and 5000 m s−1, and the P-wave velocity between 5000 and 8000 m s−1 (Fig. 7), with the constraint of having the Poisson’s ratio in the 0.1–0.45 range to make it geologically realistic. Figure 7. View largeDownload slide Comparison of effective properties computed with the 3-D non-periodic homogenization (dashed red) and the Backus theory (black) in the case of an original layered medium (blue). The two equivalent media both are anisotropic ($$\phi = \frac{V_{PV}^2}{V_{PH}^{2}}$$ and $$\xi = \frac{V_{SH}^{2}}{V_{SV}^{2}}$$ pointing out the rate of anisotropy) and agree very well with each other. Figure 7. View largeDownload slide Comparison of effective properties computed with the 3-D non-periodic homogenization (dashed red) and the Backus theory (black) in the case of an original layered medium (blue). The two equivalent media both are anisotropic ($$\phi = \frac{V_{PV}^2}{V_{PH}^{2}}$$ and $$\xi = \frac{V_{SH}^{2}}{V_{SV}^{2}}$$ pointing out the rate of anisotropy) and agree very well with each other. The homogenization of the model is performed using λ0 = 1600 m. In this case, the thickness of the outer domain (i.e. half of the filtering wavelet support) is 6400 m (Appendix C), so we take the extent of the layers equal to 15 km in order to get an inner domain in which a solution can be computed. Following a rule of thumb for the spatial discretization (Appendix E), each layer is meshed by a single layer of tetrahedral elements equipped with interpolation functions of degree 3. As expected from Backus (1962), the resulting effective medium is transversely isotropic. In Fig. 7, we represent it in terms of vertically- and horizontally polarized wave velocity:   \begin{eqnarray*} V_{PV} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{zzzz}}{\rho ^{\varepsilon _0\star }}}, \nonumber \\ V_{PH} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{xxxx}}{\rho ^{\varepsilon _0\star }}} = \sqrt{\frac{C^{\varepsilon _0\star }_{yyyy}}{\rho ^{\varepsilon _0\star }}}, \nonumber \\ V_{SV} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{xzxz}}{\rho ^{\varepsilon _0\star }}} = \sqrt{\frac{C^{\varepsilon _0\star }_{yzyz}}{\rho ^{\varepsilon _0\star }}}, \nonumber \\ V_{SH} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{xyxy}}{\rho ^{\varepsilon _0\star }}}, \nonumber \end{eqnarray*} z being the vertical, layered direction, and (x, y) defining the horizontal plane. Moreover, we emphasize the amount of anisotropy by plotting $$\phi = \frac{V_{PV}^2}{V_{PH}^2}$$ and $$\xi = \frac{V_{SH}^2}{V_{SV}^2}$$. The effective properties computed with the non-periodic homogenization in the layered case are expected to be the same as those proposed by Backus (1962). A comparison between the two solutions (Fig. 7) shows a very good match, meaning that our 3-D implementation of the non-periodic homogenization works properly in the present case. As mentioned in Part 3, no solution can be computed at the border of the domain because the convolution involved in the filtering operation is not possible there. Nonetheless, extending the domain with a relevant material to make the convolution possible and to get satisfying effective properties at the border would be possible. Simply extending the boundary values of density and elastic coefficients would lead to the zeroth-order solution (Capdeville & Marigo 2007). To reach the first order, a continuous periodic extension (in which the boundary actually acts as a mirror) would have to be used (Leptev 2005; Capdeville & Marigo 2007). We do not implement such extensions in the present work. 4.2 Homogenization of random cubes To push our validation further, we challenge our homogenization code to a highly heterogeneous medium made of small elastic cubes with random isotropic properties. Each cube is 1 km3 large. 100 cubes are considered in each direction, which gives rise to a large cubic volume made of 1 000 000 random cubes. As shown in Fig. 8(a), this cubic volume is embedded in a 13 km thick homogeneous medium. Such a thickness corresponds to the support of the filtering wavelet which will be used in the homogenization process. In each small cube, the properties are randomly picked between 2000 and 4000 kg m−3 for the density, 2500 and 5000 m s−1 for the S-wave velocity, and 4000 and 8000 m s−1 for the P-wave velocity. As in the previous example, the Poisson’s ratio is constrained in the 0.1–0.45 range. Figure 8. View largeDownload slide (a) Cut in the random cubes. The black lines emphasize the edges of the whole domain. (b) Cut in the homogenized medium. (c) Snapshot of the L2 norm of the wavefield u generated through the random cubes by a force along the z-axis at point xS marked by the blue star. (d) Snapshot of the L2 norm of the wavefield u0 generated through the homogenized medium by the exact same force. Figure 8. View largeDownload slide (a) Cut in the random cubes. The black lines emphasize the edges of the whole domain. (b) Cut in the homogenized medium. (c) Snapshot of the L2 norm of the wavefield u generated through the random cubes by a force along the z-axis at point xS marked by the blue star. (d) Snapshot of the L2 norm of the wavefield u0 generated through the homogenized medium by the exact same force. The homogenization of the random cubes is performed using λ0 = 1600 m (i.e. λm = 8000 m and ε0 = 0.2). The spatial discretization is ensured by dividing every single cube into six tetrahedra equipped with degree 3 interpolants. Using similar elements in the homogeneous region, we end up with 12 002 256 tetrahedra and 160 747 899 degrees of freedom (three components at 53 582 633 interpolation points). To achieve the computation, the domain is split into 100 overlapping subdomains. Calculating the effective properties in a single subdomain then requires about 116 GB and 4 hr on an Intel Xeon X5680 processor (6 cores, 3.33 GHz, 12 MB Cache, 6.4 GT s−1). The obtained homogenized model is shown in Fig. 8(b). To assess the quality of the effective medium calculated with the non-periodic homogenization technique, we perform a seismic wave simulation in it and we compare the obtained seismograms with traces computed in the initial model (i.e. the random cubes). In other words, we solve the initial problem   $$\left\lbrace \begin{array}{l}\rho \ddot{\boldsymbol {u}} - \boldsymbol {\nabla }\cdot \boldsymbol {\sigma } = \boldsymbol {f} \\ \boldsymbol {\sigma } = {\bf C}: \left[ \frac{1}{2} (\boldsymbol {\nabla }\boldsymbol {u} + {^t}\boldsymbol {\nabla }\boldsymbol {u}) \right] \end{array} \right.$$ (34)and the homogenized problem   $$\left\lbrace \begin{array}{l}\rho ^{\varepsilon _0\star } \ddot{\boldsymbol {u}}_0 - \boldsymbol {\nabla }\cdot \boldsymbol {\sigma }_0 = \boldsymbol {f} \\ \boldsymbol {\sigma }_0 = {\bf C}^{\varepsilon _0\star } : \left[ \frac{1}{2} (\boldsymbol {\nabla }\boldsymbol {u}_0 + {^t}\boldsymbol {\nabla }\boldsymbol {u}_0) \right] \end{array} \right.$$ (35)to estimate the quality of $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$ through the comparison of u0 to u. In eq. (34) and (35), $$\,\ddot{ }\,$$ represents the second time-derivative and f is the external force. This latter is chosen as a simple Ricker function R(t) along the z-axis at a given point xS in the homogeneous region: f(x, t) = R(t) δ(x − xS) ez. The dominant frequencies of R(t) is chosen to be equal to 0.2 Hz. The two simulations are performed using a spectral element method with PML-type absorbing boundaries (Cupillard et al.2012). Snapshots of the two obtained wavefields are shown in Figs 8(c) and (d). They look very similar, suggesting that our homogenized model is an accurate equivalent medium for the seismic wave propagation. We carefully compare u0 to u by looking at seismograms calculated at two particular receivers denoted by A and B in Figs 8(c) and (d). Receiver A is on a P-wave nodal plane, 92 km away from the source. The ballistic S-wave, which contains most of the seismic energy, appears on the z-component (Fig. 9). We see that this wave is very well-retrieved by the homogenized model, the difference between u0 and u (i.e. the residual) reaching 8   per  cent at most. Apart from the S-wave, a scattered wavefield is observed on the three components. Our homogenized model also reconstructs this wavefield very well. To emphasize the relevance of the homogenized solution, Fig. 9 also shows seismograms computed in a medium obtained by just filtering the initial density and elastic tensor. The displacement and the stress propagating in this medium are denoted by $$\boldsymbol {u}_\mathcal {F}$$ and $$\boldsymbol {\sigma }_\mathcal {F}$$, respectively. By definition, these two fields verify   $$\left\lbrace \begin{array}{l}\mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace \ddot{\boldsymbol {u}}_\mathcal {F} - \boldsymbol {\nabla }\cdot \boldsymbol {\sigma }_\mathcal {F} = \boldsymbol {f} \\ \boldsymbol {\sigma }_\mathcal {F} = \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf C}\rbrace : \left[ \frac{1}{2} (\boldsymbol {\nabla }\boldsymbol {u}_\mathcal {F} + {^t}\boldsymbol {\nabla }\boldsymbol {u}_\mathcal {F}) \right]. \end{array} \right.$$ (36)The waveforms of $$\boldsymbol {u}_\mathcal {F}$$ at receiver A do not fit the wavefield u in the initial model. Neither the phase nor the amplitude is properly reconstructed by the filtered medium, meaning that this latter does not hold the correct effective properties for the seismic wave propagation. Similar features are observed at receiver B (Fig. 10): $$\boldsymbol {u}_\mathcal {F}$$ is far from u whereas the homogenized solution accurately recovers the whole seismograms, including a ballistic P-wave which appears at this receiver location. Figure 9. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver A (Figs 8c and d). The maximum and dominant frequencies of the Ricker function R(t) emitted at the source are equal to 0.5 and 0.2 Hz, respectively. Figure 9. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver A (Figs 8c and d). The maximum and dominant frequencies of the Ricker function R(t) emitted at the source are equal to 0.5 and 0.2 Hz, respectively. Figure 10. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver B (Figs 8c and d). Contrary to receiver A, B records a ballistic P-wave. Figure 10. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver B (Figs 8c and d). Contrary to receiver A, B records a ballistic P-wave. The overall difference between u0 and u can be evaluated quantitatively using the error   $$E_0 = \frac{1}{50} \sum _{r=1}^{50} \sqrt{ \frac{\int _0^T (\boldsymbol {u}_0-\boldsymbol {u})^2(\boldsymbol {x}_r) \, \text{d}t}{\int _0^T \boldsymbol {u}^2(\boldsymbol {x}_r) \, \text{d}t} } \, ,$$ (37)where xr is the randomly chosen location of receiver r. In eq. (37), T is the duration of the simulated propagation. It is equal to 55 s here. For ε0 = 0.2 (i.e. λ0 = 1600 m), we obtain E0 ≃ 0.006. Computing the same kind of error for the naive solution $$\boldsymbol {u}_\mathcal {F}$$, we get $$E_\mathcal {F} \simeq 0.055 \simeq 9E_0$$. This result again emphasizes the much higher accuracy of the homogenized solution u0. We also compute E0 and $$E_\mathcal {F}$$ for ε0 = 0.4 (i.e. λ0 = 3200 m), ε0 = 0.8 (i.e. λ0 = 6400 m) and ε0 = 1.6 (i.e. λ0 = 12 800 m). This allows us to study how fast u0 and $$\boldsymbol {u}_\mathcal {F}$$ converge to u. From the homogenization theory, we expect E0 = O(ε0). Fig. 11 shows that we actually get $$E_0 \simeq O(\varepsilon _0^{3/2})$$. This unexpectedly fast convergence of u0 toward u could be due to the weak contribution of the higher-order terms in the present case study, as suggested by Capdeville et al. (2010b) in 2-D. Plotting also $$E_\mathcal {F}$$ in Fig. 11, it appears that the convergence of $$\boldsymbol {u}_\mathcal {F}$$ is way poorer. Figure 11. View largeDownload slide E0 (red) and $$E_\mathcal {F}$$ (green) as a function of ε0. It is clear from this plot that u0 converges much faster than $$\boldsymbol {u}_\mathcal {F}$$ toward the target wavefield u. Figure 11. View largeDownload slide E0 (red) and $$E_\mathcal {F}$$ (green) as a function of ε0. It is clear from this plot that u0 converges much faster than $$\boldsymbol {u}_\mathcal {F}$$ toward the target wavefield u. The spectral element simulation of the wavefield u in the initial model requires 126 hexahedra in each direction, that is, 2 000 376 elements in total. This is because a spectral element mesh has to honour all the physical discontinuities of the model in study to make the computation accurate, so each 1 km × 1 km × 1 km random cube has to be captured by a hexahedron. As a consequence, the obtained mesh highly oversamples the wavefield: because λm = 8 km, elements as large as 8 km × 8 km × 8 km with a polynomial order larger than 4 would be sufficient if the medium was smooth enough (i.e. if it only contained scales larger than λm). Because of the discontinuities, we here get a 512 times denser mesh, which yields a 4096 times higher numerical cost (a factor 512 in space times a factor 8 in time because of eq. 1). Computing a 55 s long simulation of u then takes about 6 hr on ten Intel Xeon X5680 processors. In homogenized media, such as those computed here using ε0 = 0.2, 0.4, 0.8 and 1.6, coarser spectral element meshes can be used, then decreasing the computation cost. When using ε0 = 0.4 for instance, we get an effective medium of the random cubes that only contains scales larger than λ0 = 3.2 km, so 3.2 km3 large hexahedra then suit an accurate spectral element simulation of u0, which is 105 times less numerically demanding than the simulation of u. Nevertheless, all the wave simulations presented in this paper are performed within the same fine hexahedral mesh to avoid possible numerical biases and thus make all our comparisons totally fair. 4.3 Homogenization of a realistic geological model In this section, we use a subsurface model of the Ribaute area in France (Caumon et al.2009) as an illustration of the 3-D non-periodic homogenization technique applied to a realistic geological medium (Fig. 12a). Simulating waves in the model as it is (i.e. composed of multiple faulted and folded horizons) is extremely challenging because of the fine grid required to accommodate thin layers and complex geometries formed by the discontinuities. Such a fine grid indeed makes the computation cost of any wave simulation enormous. Moreover, in the context of the spectral element method, generating a hexahedral mesh which honours all the horizons and faults, if possible, demands tremendous efforts. For these reasons, we do not compute reference waveforms here; we just perform a spectral element simulation in the homogenized model to show how convenient working with effective properties is, and we compare the results with a simulation in a model obtained by brutally filtering the initial density and elastic tensor to exemplify that the choice of the upscaling technique matters. Figure 12. View largeDownload slide (a) Structural model of a highly faulted and folded region near Ribaute, Southern France. (b) Adaptive tetrahedral mesh of the model. This mesh is used to perform the homogenization. The background colour corresponds to the (randomly chosen) isotropic S-wave velocity within each layer of the model. (c) One of the S-wave velocities resulting from the homogenization. Another S-wave velocity is shown in Fig. 13. The difference between these two velocities is plotted in (d), which emphasizes the seismic anisotropy produced by the structure. Figure 12. View largeDownload slide (a) Structural model of a highly faulted and folded region near Ribaute, Southern France. (b) Adaptive tetrahedral mesh of the model. This mesh is used to perform the homogenization. The background colour corresponds to the (randomly chosen) isotropic S-wave velocity within each layer of the model. (c) One of the S-wave velocities resulting from the homogenization. Another S-wave velocity is shown in Fig. 13. The difference between these two velocities is plotted in (d), which emphasizes the seismic anisotropy produced by the structure. The initial model is made of five homogeneous isotropic layers (Fig. 12b). To compute its effective properties using our finite-element code, it is meshed with VorteXLib from RINGMesh (Pellerin et al.2017). We choose λ0 = 30 m and polynomial interpolants of degree 2, so the optimal volume for the tetrahedra is about 200 m3 (Appendix E). While smaller elements are needed around the discontinuities for capturing their complex geometry, larger elements are allowed where the discontinuities have no influence, that is, at a distance larger than the size of the filtering wavelet support. We therefore take advantage of the refinement-derefinement technique proposed by Botella (2016) to generate an adaptive mesh which minimizes the finite-element computation cost. We end up with a 9 077 300 element mesh and 35 703 579 degrees of freedom. The homogenized model is obtained in 9 min on 120 Intel Xeon E5-2683 v4 (16 cores, 2.10 GHz, 40 MB Cache, 9.6 GT s−1). The result is shown in Fig. 12(c). We have plotted the $$\sqrt{\frac{C^{\varepsilon _0\star }_{xyxy}}{\rho ^{\varepsilon _0\star }}}$$ component there to make the figure comparable to Fig. 12(b) where the isotropic S-wave velocity of the initial medium is represented. As expected, the homogeneous areas (i.e. the interior of the layers) are not changed by the homogenization process, whereas all the discontinuities (i.e. the faults and the interfaces between the layers) are smoothed. To illustrate the S-anisotropy produced by these discontinuities, $$\frac{ \sqrt{C^{\varepsilon _0\star }_{xyxy}} - \sqrt{C^{\varepsilon _0\star }_{xzxz}} }{ \sqrt{\rho ^{\varepsilon _0\star }} }$$ is plotted in Fig. 12(d). We observe that the anisotropy can reach 20   per  cent where high velocity contrasts (i.e. at the bottom and the top of the fastest layer) occur. Because the homogenized model is smooth, seismic waves within it can be simulated using a coarse mesh. As an example, a coarse hexahedral mesh supporting a spectral element simulation (Cupillard et al.2012) is shown in Fig. 13. The elements there are 60 m3 large and hold degree-8 polynomials to capture all the variations of the model. We see that the direct seismic wave front is highly deformed by the effective structure and that strong reflected and diffracted waves are generated, even though the medium contains no discontinuities. Fig. 13 also shows waveforms obtained at two stations. At station A, near the source, several S-waves reflected on discontinuities are observed. Zooming in two of them emphasizes important discrepancies between our homogenized solution and the waveform computed in a medium obtained by just filtering the density and the elastic tensor of the initial model. These discrepancies can be explained by the lack of anisotropy in this last medium. By construction, it is fully isotropic, so it cannot hold the important structure-induced anisotropy observed in our homogenized model (e.g. Fig. 12d). As a consequence, the waveforms computed in the two models at station A show significant differences in phase and amplitude. Similar observations are also made far from the source, at station B. In the light of the results presented in Section 4.2 and in previous 2-D studies (Capdeville et al.2010b; Guillot et al.2010), we can reasonably think that the waveforms computed in the homogenized model are the most accurate, but we cannot assess it because we do not provide any reference solution here. The main point of this section is to show an application of our homogenization code to a 3-D realistic geological model and to put the importance of the structure-induced anisotropy in evidence. Figure 13. View largeDownload slide Left: three snapshots of the L2 norm of the wavefield u0 generated through the homogenized Ribaute model by a force along the x-axis acting at the point marked by the yellow star. The wave simulation is performed using a spectral element method on a regular hexahedral mesh. Top right: x-component of the displacement u0 recorded near the source, at station A. The displacement $$\boldsymbol {u}_\mathcal {F}$$ computed in a medium obtained by just filtering the density and the elastic tensor of the initial model is also plotted. As expected, the direct S-wave (arriving around t = 0.05 s) is identical in the two simulations. On the contrary, the reflected phases show significant differences. In particular, zooming in the waves reflected on the two major velocity contrasts (i.e. the bottom and the top of the fastest layer) shows that $$\boldsymbol {u}_\mathcal {F}$$ arrives earlier than u0 with a different amplitude. Bottom right: x- and z-components of both u0 and $$\boldsymbol {u}_\mathcal {F}$$ recorded far from the source, at station B. Because the direct S-wave has travelled through heterogeneities, it is now different in the two simulations. Obviously, discrepancies in the later, multiply reflected and scattered energy, are also observed. Figure 13. View largeDownload slide Left: three snapshots of the L2 norm of the wavefield u0 generated through the homogenized Ribaute model by a force along the x-axis acting at the point marked by the yellow star. The wave simulation is performed using a spectral element method on a regular hexahedral mesh. Top right: x-component of the displacement u0 recorded near the source, at station A. The displacement $$\boldsymbol {u}_\mathcal {F}$$ computed in a medium obtained by just filtering the density and the elastic tensor of the initial model is also plotted. As expected, the direct S-wave (arriving around t = 0.05 s) is identical in the two simulations. On the contrary, the reflected phases show significant differences. In particular, zooming in the waves reflected on the two major velocity contrasts (i.e. the bottom and the top of the fastest layer) shows that $$\boldsymbol {u}_\mathcal {F}$$ arrives earlier than u0 with a different amplitude. Bottom right: x- and z-components of both u0 and $$\boldsymbol {u}_\mathcal {F}$$ recorded far from the source, at station B. Because the direct S-wave has travelled through heterogeneities, it is now different in the two simulations. Obviously, discrepancies in the later, multiply reflected and scattered energy, are also observed. 5 DISCUSSION AND CONCLUSIONS Dealing with small-scale heterogeneities in seismic wave simulation is a difficult task because it usually involves enormous computation costs. To handle them, Graphics Processing Unit (GPU) and/or High Performance Computing (HPC) implementations of well-established numerical techniques such as the Spectral Element method (SEM), the Discontinuous Galerkin method (DGM) and the Finite Difference method (FDM), have been proposed (Komatitsch et al.2010; Peter et al.2011; Weiss & Shragge 2013; Gokhberg & Fichtner 2016; Remacle et al.2016; Rietmann et al.2017). Furthermore, the DGM allows local time-stepping and p-adaptivity (e.g. Dumbser et al.2007; Etienne et al.2010; Minisini et al.2013; Diaz & Grote 2015) which mitigate stability constraint (1) and therefore reduce the overall computation cost. In the context of the SEM and the FDM, Pelties et al. (2011) proposed empirical laws on the mesh size to relax the need of honouring geological discontinuities. Complementary to all these numerical advances, the use of effective properties for the seismic wave propagation can drastically reduce the computation cost related to small-scale features while preserving a good accuracy. In the recent years, the non-periodic homogenization method emerged as a general technique to compute such effective properties. We here applied it in 3-D for the first time. In part 2, we recalled the theory of the homogenization, skipping some technical details to focus on the main ideas and concepts. Then we described an efficient, embarrassingly parallel implementation of the method. In part 4, we challenged this implementation on various media, showing the high accuracy of the non-periodic homogenization and the ability of our code to handle large and complex 3-D models, with no restriction on the size and shape of the heterogeneities. The code is available at upon request. Homogenizing complex geological media to ease the numerical simulation of full seismic wavefields has straight-forward applications in seismic risk assessment (e.g. Chaljub et al.2010), survey design (e.g. Wei et al.2012), structural model validation (e.g. Irakarama et al.2017), seismic source characterization (e.g. Silwal & Tape 2016) and seismic tomography for keeping subwavelength geological details in the inversion (Fichtner et al.2013b; Capdeville & Cance 2015). Because it tells what the waves ‘see’, the non-periodic homogenization also opens important perspectives in the interpretation of full waveform inversion results. In particular, the structure-induced (or, equivalently, ‘extrinsic’ or ‘apparent’ or ‘geometric’) anisotropy (Fichtner et al.2013a; Wang et al.2013) and the structure-induced attenuation (Ricard et al.2014) could be downscaled, as recently initiated by Bodin et al. (2015), leading to probability distributions of small-scale features and to better estimations of the intrinsic anisotropy and attenuation at various scales, from the Earth’s core up to the subsurface. In addition, Capdeville et al. (2013) and Afanasiev et al. (2016) recently showed that the homogenization can help in regularizing full waveform inversion problems. The non-periodic homogenization relies on a two-scale asymptotic expansion of the displacement and the stress involved in the elastic wave equation. In this paper, we have investigated the zeroth-order term and the associated upscaled properties. As shown by Capdeville et al. (2010a,b) and Guillot et al. (2010), adding the first-order term allows retrieving site effects (i.e. small-amplitude high-frequency non-propagating signals) at the receivers. These authors also show that a correction can be applied to the external force term to take into account the effects of the local small-scale structure at the source. A recent application of this correction can be found in Burgos et al. (2016). We have claimed that our method does not require any constraint on the size and shape of the heterogeneities to be smoothed. In other words, the homogenization is able to upscale any media. The only limitation that has been noted so far occurs when trying to model subwavelength focusing in Helmholtz resonators. In this extreme case, the wavefield no longer holds a minimum wavelength, so the non-periodic homogenization fails (Zhao et al.2016). Moreover, we have mentioned that our homogenization theory is not able to compute long-wave equivalent properties of heterogeneities laying near Neumann or Dirichlet surfaces yet. A proper solution only exists for smoothing a free-surface topography on top of a homogeneous material (Capdeville & Marigo 2013) or fine horizontal layers below a flat free-surface (Capdeville & Marigo 2008). Relevant extensions of the model beyond its boundaries can also leads to accurate effective properties (Leptev 2005; Capdeville & Marigo 2007). Pushing these approaches would surely help the non-periodic homogenization theory in handling heterogeneities near various boundary shapes and conditions. Besides these physical and theoretical limitations, slight numerical weaknesses can be pointed out in our code. Even if it employs powerful algorithms and an efficient parallel scheme, it could be optimized using Basic Linear Algebra Subprograms and the Fortran column-major order more extensively, particularly when computing the stiffness matrix and when filtering the strain and the stress associated with the starting corrector. The performance of the code could also benefit from the nested homogenization technique proposed by Capdeville et al. (2015) and from an iterative solver which would be more efficient than a direct solver when dealing with very large stiffness matrices. Most of all, developing an adaptive homogenization, taking into account the fact that the minimum wavelength usually varies within the medium, would be a major advance. Acknowledgements We deeply thank Q. Liu and C. Boehm for their sharp reviews and constructive comments that significantly improve the initial manuscript. We also thank A. Fichtner, L. Guillot, P. Cance, Y. Ricard, F. Chambat, J.-P. Montagner, J. Virieux, A. Kutsenko, Y. Masson, G. Caumon and the European COST action TIDES (ES1401) for emulative remarks and discussions. We are also grateful to A. Botella for advising some powerful search algorithms and for providing 3-D adaptive meshes which significantly reduced the computation cost of the homogenization. Many thanks to A. Borghi as well, who completely reset the Cauchy cluster at GeoRessources. Part of the calculations has also been led on the S-CAPAD cluster at IPGP, on the Centre de Calcul Intensif des Pays de la Loire (CCIPL) computers, and on the EXPLOR centre hosted by the University of Lorraine. Financial supports from the RING-Gocad consortium and the ANR projects ANR-10-BLAN-613 and ANR-16-CE31-0022-01 also enabled this work. REFERENCES Abdulle A., Grote M.J., 2011. Finite element heterogeneous multiscale method for the wave equation, Multiscale Model. Simul. , 9( 2), 766– 792. https://doi.org/10.1137/100800488 Google Scholar CrossRef Search ADS   Afanasiev M., Boehm C., May D., Fichtner A., 2016. Using effective medium theory to better constrain full waveform inversion, First Break , 34( 7), 93– 94. Allaire G., 1992. Homogenization and two-scale convergence, SIAM J. Math. Anal. , 23, 1482– 1518. https://doi.org/10.1137/0523084 Google Scholar CrossRef Search ADS   Allaire G., Capdeboscq Y., 2000. Homogenization of a spectral problem in neutronic multigroup diffusion, Comput. Methods Appl. Mech. Eng. , 187( 1–2), 91– 117. https://doi.org/10.1016/S0045-7825(99)00112-7 Google Scholar CrossRef Search ADS   Allaire G., Habibi Z., 2013. Homogenization of a conductive, convective, and radiative heat transfer problem in a heterogeneous domain, SIAM J. Math. Anal. , 45( 3), 1136– 1178. https://doi.org/10.1137/110849821 Google Scholar CrossRef Search ADS   Bacigalupo A., Gambarotta L., 2014. Second-gradient homogenized model for wave propagation in heterogeneous periodic media, Int. J. Solids Struct. , 51( 5), 1052– 1065. https://doi.org/10.1016/j.ijsolstr.2013.12.001 Google Scholar CrossRef Search ADS   Backus G., 1962. Long-wave elastic anisotropy produced by horizontal layering, J. geophys. Res. , 67( 11), 4427– 4440. https://doi.org/10.1029/JZ067i011p04427 Google Scholar CrossRef Search ADS   Bensoussan A., Lions J.-L., Papanicolaou G., 1978. Asymptotic Analysis for Periodic Structures , North-Holland. Bodin T., Capdeville Y., Romanowicz B., Montagner J.-P., 2015. Interpreting radial anisotropy in global and regional tomographic models, in The Earth’s Heterogeneous Mantle , pp. 105– 144, eds Khan A., Deschamps F., Springer. Google Scholar CrossRef Search ADS   Botella A., 2016. Génération de maillages non structurés volumiques de modèles géologiques pour la simulation de phénomènes physiques, PhD thesis , Université de Lorraine. Boutin C., Auriault J.L., 1993. Rayleigh scattering in elastic composite materials, Int. J. Eng. Sci. , 31( 12), 1669– 1689. https://doi.org/10.1016/0020-7225(93)90082-6 Google Scholar CrossRef Search ADS   Brossier R., Operto S., Virieux J., 2015. Velocity model building from seismic reflection data by full waveform inversion, Geophys. Prospect. , 63( 2), 354– 367. https://doi.org/10.1111/1365-2478.12190 Google Scholar CrossRef Search ADS   Burgos G., Capdeville Y., Guillot L., 2016. Homogenized moment tensor and the effect of near-field heterogeneities on nonisotropic radiation in nuclear explosion, J. geophys. Res. , 121( 6), 4366– 4389. https://doi.org/10.1002/2015JB012744 Google Scholar CrossRef Search ADS   Cance P., Capdeville Y., 2015. Validity of the acoustic approximation for elastic waves in heterogeneous media, Geophysics , 80( 4), T161– T173. https://doi.org/10.1190/geo2014-0397.1 Google Scholar CrossRef Search ADS   Capdeville Y., Cance P., 2015. Residual homogenization for elastic wave propagation in complex media, Geophys. J. Int. , 200, 984– 997. https://doi.org/10.1093/gji/ggu452 Google Scholar CrossRef Search ADS   Capdeville Y., Marigo J., 2007. Second order homogenization of the elastic wave equation for non-periodic layered media, Geophys. J. Int. , 170, 823– 838. https://doi.org/10.1111/j.1365-246X.2007.03462.x Google Scholar CrossRef Search ADS   Capdeville Y., Marigo J., 2008. Shallow layer correction for spectral element like methods, Geophys. J. Int. , 172, 1135– 1150. https://doi.org/10.1111/j.1365-246X.2007.03703.x Google Scholar CrossRef Search ADS   Capdeville Y., Marigo J.-J., 2013. A non-periodic two scale asymptotic method to take account of rough topographies for 2-D elastic wave propagation, Geophys. J. Int. , 192( 1), 163– 189. https://doi.org/10.1093/gji/ggs001 Google Scholar CrossRef Search ADS   Capdeville Y., Guillot L., Marigo J., 2010a. 1-D non-periodic homogenization for the seismic wave equation, Geophys. J. Int. , 181, 907– 910. Capdeville Y., Guillot L., Marigo J., 2010b. 2-D non-periodic homogenization to upscale elastic media for P-SV waves, Geophys. J. Int. , 182, 903– 922. https://doi.org/10.1111/j.1365-246X.2010.04636.x Google Scholar CrossRef Search ADS   Capdeville Y., Stutzmann E., Wang N., Montagner J.-P., 2013. Residual homogenization for seismic forward and inverse problems in layered media, Geophys. J. Int. , 194( 1), 470– 487. https://doi.org/10.1093/gji/ggt102 Google Scholar CrossRef Search ADS   Capdeville Y., Zhao M., Cupillard P., 2015. Fast Fourier homogenization for elastic wave propagation in complex media, Wave Motion , 54, 170– 186. https://doi.org/10.1016/j.wavemoti.2014.12.006 Google Scholar CrossRef Search ADS   Carcione J., Picotti S., Santos J., 2012. Numerical experiments of fracture-induced velocity and attenuation anisotropy, Geophys. J. Int. , 191( 3), 1179– 1191. Casarotti E., Stupazzini M., Lee S., Komatitsch D., Piersanti A., Tromp J., 2008. CUBIT and seismic wave propagation based upon the spectral-element method: an advanced unstructured mesher for complex 3D geological media, in Proceedings of the 16th International Meshing Roundtable , pp. 579– 597, eds Brewer M.L., Marcum D., Springer, Berlin. Google Scholar CrossRef Search ADS   Caumon G., Collon P., Le Carlier de Veslud C., Sausse J., Viseur S., 2009. Surface-based 3D modeling of geological structures, Math. Geosci. , 41( 8), 927– 945. https://doi.org/10.1007/s11004-009-9244-2 Google Scholar CrossRef Search ADS   Chaljub E., Moczo P., Tsuno S., Bard P.-Y., Kristek J., Käser M., Stupazzini M., Kristekova M., 2010. Quantitative comparison of four numerical predictions of 3D ground motion in the Grenoble valley, France, Bull. seism. Soc. Am. , 100( 4), 1427– 1455. https://doi.org/10.1785/0120090052 Google Scholar CrossRef Search ADS   Cormen T.H., Leiserson C.E., Rivest R.L., Stein C., 2009. Introduction to Algorithms , 3rd edn, The MIT Press. Courant R., Friedrichs K., Lewy H., 1928. Uber die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann. , 100, 32– 74. https://doi.org/10.1007/BF01448839 Google Scholar CrossRef Search ADS   Cupillard P., Delavaud E., Burgos G., Festa G., Vilotte J.-P., Capdeville Y., Montagner J.-P., 2012. RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale, Geophys. J. Int. , 188, 1203– 1220. https://doi.org/10.1111/j.1365-246X.2011.05311.x Google Scholar CrossRef Search ADS   Dal Maso G., 1993. An Introduction to Γ-convergence , Birkhäuser. Google Scholar CrossRef Search ADS   Diaz J., Grote M.J., 2015. Multi-level explicit local time-stepping methods for second-order wave equations, Comput. Methods Appl. Mech. Eng. , 291, 240– 265. https://doi.org/10.1016/j.cma.2015.03.027 Google Scholar CrossRef Search ADS   Dumbser M., Käser M., Toro E.F., 2007. An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes—V. Local time stepping and p-adaptivity, Geophys. J. Int. , 171( 2), 695– 717. https://doi.org/10.1111/j.1365-246X.2007.03427.x Google Scholar CrossRef Search ADS   Dumontet H., 1990. Homogénéisation et effets de bords dans les matériaux composites, PhD thesis , Univertité Paris 6. Engquist B., Holst H., Runborg O., 2009. Multi-scale methods for wave propagation in heterogeneous media , preprint arXiv:0911.2638. Etienne V., Chaljub E., Virieux J., Glinsky N., 2010. An hp-adaptive discontinuous Galerkin Þnite-element method for 3-D elastic wave modelling, Geophys. J. Int. , 183, 941– 962. https://doi.org/10.1111/j.1365-246X.2010.04764.x Google Scholar CrossRef Search ADS   Felippa C., 2004. A compendium of FEM integration formulas for symbolic work, Eng. Comput. , 21( 8), 867– 890. https://doi.org/10.1108/02644400410554362 Google Scholar CrossRef Search ADS   Fichtner A., Igel H., 2008. Efficient numerical surface wave propagation through the optimization of discrete crustal models - a technique based on non-linear dispersion curve matching (DCM), Geophys. J. Int. , 173( 2), 519– 533. https://doi.org/10.1111/j.1365-246X.2008.03746.x Google Scholar CrossRef Search ADS   Fichtner A., Kennett B.L.N., Igel H., Bunge H.-P., 2008. Theoretical background for continental and global scale full-waveform inversion in the time-frequency domain, Geophys. J. Int. , 175( 8), 665– 685. https://doi.org/10.1111/j.1365-246X.2008.03923.x Google Scholar CrossRef Search ADS   Fichtner A., Kennett B.L., Trampert J., 2013a. Separating intrinsic and apparent anisotropy, Phys. Earth planet. Inter. , 219, 11– 20. https://doi.org/10.1016/j.pepi.2013.03.006 Google Scholar CrossRef Search ADS   Fichtner A., Trampert J., Cupillard P., Saygin E., Taymaz T., Capdeville Y., Villasenor A., 2013b. Multi-scale full waveform inversion, Geophys. J. Int. , 194, 534– 556. https://doi.org/10.1093/gji/ggt118 Google Scholar CrossRef Search ADS   Fish J., Chen W., 2001. Higher-order homogenization of initial/boundary-value problem, J. Eng. Mech. , 127( 12), 1223– 1230. https://doi.org/10.1061/(ASCE)0733-9399(2001)127:12(1223) Google Scholar CrossRef Search ADS   Fish J., Chen W., 2004. Space-time multiscale model for wave propagation in heterogeneous media, Comput. Methods Appl. Mech. Eng. , 193, 4837– 4856. https://doi.org/10.1016/j.cma.2004.05.006 Google Scholar CrossRef Search ADS   Gao K., Chung E., Gibson R., Fu S., Efendiev Y., 2015. A numerical homogenization method for heterogeneous, anisotropic elastic media based on multiscale theory, Geophysics , 80( 4), D385– D401. https://doi.org/10.1190/geo2014-0363.1 Google Scholar CrossRef Search ADS   Gokhberg A., Fichtner A., 2016. Full-waveform inversion on heterogeneous HPC systems, Comput. Geosci. , 89, 260– 268. https://doi.org/10.1016/j.cageo.2015.12.013 Google Scholar CrossRef Search ADS   Grechka V., 2003. Effective media: A forward modeling view, Geophysics , 68( 6), 2055– 2062. https://doi.org/10.1190/1.1635059 Google Scholar CrossRef Search ADS   Grechka V., 2007. Multiple cracks in vti rocks: Effective properties and fracture characterization, Geophysics , 72( 5), D81– D91. https://doi.org/10.1190/1.2751500 Google Scholar CrossRef Search ADS   Grechka V., Kachanov M., 2006. Effective elasticity of fractured rocks: A snapshot of the work in progress, Geophysics , 71( 6), W45– W58. https://doi.org/10.1190/1.2360212 Google Scholar CrossRef Search ADS   Guillot L., Capdeville Y., Marigo J., 2010. 2-D non-periodic homogenization of the elastic wave equation: SH case, Geophys. J. Int. , 182, 1438– 1454. https://doi.org/10.1111/j.1365-246X.2010.04688.x Google Scholar CrossRef Search ADS   Hashin Z., Shtrikman S., 1963. A variational approach to the elastic behavior of multiphase materials, J. Mech. Phys. Solids , 11, 127– 140. https://doi.org/10.1016/0022-5096(63)90060-7 Google Scholar CrossRef Search ADS   Hill R., 1965. A self–consistent mechanics of composit materials, J. Mech. Phys. Solids , 13, 213– 222. https://doi.org/10.1016/0022-5096(65)90010-4 Google Scholar CrossRef Search ADS   Hornung U., 1997. Homogenization and Porous Media , Vol. 6 of Interdisciplinary Applied Mathematics Series, Springer Verlag. Google Scholar CrossRef Search ADS   Hughes T.J., 2012. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis , Courier Corporation. Irakarama M., Cupillard P., Caumon G., Sava P., 2017. Appraising structural interpretations using seismic data misfit functionals, in 79th EAGE Conference and Exhibition , EAGE. Jordan T.H., 2015. An effective medium theory for three-dimensional elastic heterogeneities, Geophys. J. Int. , 203( 2), 1343– 1354. https://doi.org/10.1093/gji/ggv355 Google Scholar CrossRef Search ADS   Komatitsch D., Erlebacher G., Göddeke D., Michéa D., 2010. High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster, J. Comput. Phys. , 229( 20), 7692– 7714. https://doi.org/10.1016/j.jcp.2010.06.024 Google Scholar CrossRef Search ADS   Kutsenko A., Shuvalov A., Norris A., 2013. On the quasistatic effective elastic moduli for elastic waves in three-dimensional phononic crystals, J. Mech. Phys. Solids , 61( 11), 2240– 2259. https://doi.org/10.1016/j.jmps.2013.06.003 Google Scholar CrossRef Search ADS   Lekić V., Panning M., Romanowicz B., 2010. A simple method for improving crustal corrections in waveform tomography, Geophys. J. Int. , 182, 265– 278. Leptev V., 2005. Two-scale extensions for non-periodic coefficients, arXiv:math.AP/0512123. Lin C., Saleh R., Milkereit B., Liu Q., 2017. Effective media for transversely isotropic models based on homogenization theory: With applications to borehole sonic logs, Pure appl. Geophys. , 174( 7), 2631– 2647. https://doi.org/10.1007/s00024-017-1565-3 Google Scholar CrossRef Search ADS   Mauge C., Kachanov M., 1994. Effective elastic properties of an anisotropic material with arbitrarily oriented interacting cracks, J. Mech. Phys. Solids , 42, 561– 584. https://doi.org/10.1016/0022-5096(94)90052-3 Google Scholar CrossRef Search ADS   Métivier L., Bretaudeau F., Brossier R., Operto S., Virieux J., 2013. Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures, Geophys. Prospect. , 62( 6), 1353– 1375. Google Scholar CrossRef Search ADS   Minisini S., Zhebel E., Kononov A., Mulder W.A., 2013. Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media, Geophysics , 78( 3), T67– T77. https://doi.org/10.1190/geo2012-0252.1 Google Scholar CrossRef Search ADS   Moulinec H., Suquet P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Eng. , 157( 1), 69– 94. https://doi.org/10.1016/S0045-7825(97)00218-1 Google Scholar CrossRef Search ADS   Nguetseng G., 1989. A general convergence result for a functional related to the theory of homogenization, SIAM J. Math. Anal. , 20( 3), 608– 623. https://doi.org/10.1137/0520043 Google Scholar CrossRef Search ADS   Papanicolaou G.C., Varadhan S. R.S., 1981. Boundary value problems with rapidly oscillating random coefficients, in Random Fields, Vol. I, II (Esztergom, 1979) , Vol. 27 of Colloq. Math. Soc. János Bolyai, pp. 835– 873, North-Holland. Parnell W., Abrahams I., 2008. Homogenization for wave propagation in periodic fibre-reinforced media with complex microstructure. I - Theory, J. Mech. Phys. Solids , 56( 7), 2521– 2540. https://doi.org/10.1016/j.jmps.2008.02.003 Google Scholar CrossRef Search ADS   Pellerin J., Botella A., Bonneau F., Mazuyer A., Chauvin B., Levy B., Caumon G., 2017. RINGMesh: a programming library for developing mesh-based geomodeling applications, Comput. Geosci. , 104, 93– 100. https://doi.org/10.1016/j.cageo.2017.03.005 Google Scholar CrossRef Search ADS   Pelties C., Käser M., Hermann V., Castro C.E., 2011. Regular versus irregular meshing for complicated models and their effect on synthetic seismograms, Geophys. J. Int. , 183, 1031– 1051. https://doi.org/10.1111/j.1365-246X.2010.04777.x Google Scholar CrossRef Search ADS   Peter D. et al.  , 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, Geophys. J. Int. , 186, 721– 739. https://doi.org/10.1111/j.1365-246X.2011.05044.x Google Scholar CrossRef Search ADS   Plessix R.-E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophys. J. Int. , 167( 2), 495– 503. https://doi.org/10.1111/j.1365-246X.2006.02978.x Google Scholar CrossRef Search ADS   Pratt R., Shin C., Hicks G., 1998. Gauss–Newton and full Newton methods in frequency domain seismic waveform inversion, Geophys. J. Int. , 133, 341– 362. https://doi.org/10.1046/j.1365-246X.1998.00498.x Google Scholar CrossRef Search ADS   Remacle J.-F., Gandham R., Warburton T., 2016. GPU accelerated spectral finite elements on all-hex meshes, J. Comput. Phys. , 324, 246– 257. https://doi.org/10.1016/j.jcp.2016.08.005 Google Scholar CrossRef Search ADS   Reuss A., 1929. Berechnung der Fliessgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle, J. Appl. Math. Mech. , 9( 1), 49– 58. Ricard Y., Durand S., Montagner J.-P., Chambat F., 2014. Is there seismic attenuation in the mantle?, Earth planet. Sci. Lett. , 388, 257– 264. https://doi.org/10.1016/j.epsl.2013.12.008 Google Scholar CrossRef Search ADS   Rietmann M., Grote M., Peter D., Schenk O., 2017. Newmark local time stepping on high-performance computing architectures, J. Comput. Phys. , 334, 308– 326. https://doi.org/10.1016/j.jcp.2016.11.012 Google Scholar CrossRef Search ADS   Sanchez-Palencia E., 1980. Non-homogenous Media and Vibration Theory, Vol. 127 of Lecture notes in physics , Springer-Verlag. Santugini-Repiquet K., 2007. Homogenization of the demagnetization field operator in periodically perforated domains, J. Math. Anal. Appl. , 334( 1), 502– 516. https://doi.org/10.1016/j.jmaa.2007.01.001 Google Scholar CrossRef Search ADS   Sayers C., 1998. Long-wave seismic anisotropy of heterogeneous reservoirs, Geophys. J. Int. , 132, 667– 673. https://doi.org/10.1046/j.1365-246X.1998.00456.x Google Scholar CrossRef Search ADS   Sayers C., Kachanov M., 1991. A simple technique for finding effective elastic constants of cracked solids for arbitrary crack orientation statistics, Int. J. Solids Struct. , 27( 6), 671– 680. https://doi.org/10.1016/0020-7683(91)90027-D Google Scholar CrossRef Search ADS   Sayers C., Kachanov M., 1995. Microcrack-induced elastic wave anisotropy of brittle rocks, J. geophys. Res. , 100( B3), 4149– 4156. https://doi.org/10.1029/94JB03134 Google Scholar CrossRef Search ADS   Schenk O., Gärtner K., 2006. On fast factorization pivoting methods for symmetric indefinite systems, Electron. Trans. Numer. Anal. , 23, 158– 179. Schoenberg M., Helbig K., 1997. Orthorhombic media: modeling elastic wave behavior in a vertically fractured earth, Geophysics , 62( 6), 1954– 1974. https://doi.org/10.1190/1.1444297 Google Scholar CrossRef Search ADS   Schoenberg M., Muir F., 1989. A calculus for finely layered anisotropic media, Geophysics , 54, 581– 589. https://doi.org/10.1190/1.1442685 Google Scholar CrossRef Search ADS   Schoenberg M., Sayers C.M., 1995. Seismic anisotropy of fractured rock, Geophysics , 60( 1), 204– 211. https://doi.org/10.1190/1.1443748 Google Scholar CrossRef Search ADS   Silwal V., Tape C., 2016. Seismic moment tensors and estimated uncertainties in southern Alaska, J. geophys. Res. , 121( 4), 2772– 2797. https://doi.org/10.1002/2015JB012588 Google Scholar CrossRef Search ADS   Tromp J., Tape C., Liu Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, Geophys. J. Int. , 160, 195– 216. https://doi.org/10.1111/j.1365-246X.2004.02453.x Google Scholar CrossRef Search ADS   Virieux J., Calandra H., Plessix R.-E., 2011. A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging, Geophys. Prospect. , 59( 5), 794– 813. https://doi.org/10.1111/j.1365-2478.2011.00967.x Google Scholar CrossRef Search ADS   Voigt W., 1889. Uber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper, Ann. Phys. , 274( 12), 573– 587. https://doi.org/10.1002/andp.18892741206 Google Scholar CrossRef Search ADS   Wang N., Montagner J.-P., Fichtner A., Capdeville Y., 2013. Intrinsic versus extrinsic seismic anisotropy: The radial anisotropy in reference earth models, Geophys. Res. Lett. , 40( 16), 4284– 4288. https://doi.org/10.1002/grl.50873 Google Scholar CrossRef Search ADS   Warner M., Guasch L., 2016. Adaptive waveform inversion: Theory, Geophysics , 81( 6), R429– R445. https://doi.org/10.1190/geo2015-0387.1 Google Scholar CrossRef Search ADS   Wei W., Fu L., Blacquière G., 2012. Fast multifrequency focal beam analysis for 3D seismic acquisition geometry, Geophysics , 77( 2), P11– P21. https://doi.org/10.1190/geo2010-0327.1 Google Scholar CrossRef Search ADS   Weiss R.M., Shragge J., 2013. Solving 3D anisotropic elastic wave equations on parallel GPU devices, Geophysics , 78( 2), F7– F15. https://doi.org/10.1190/geo2012-0063.1 Google Scholar CrossRef Search ADS   Worth D., Greenough C., Chin S., 2012. Pragmatic software engineering for computational science, in Handbook of Research on Computational Science and Engineering , Vol. 1, pp. 119– 149, eds Leng J., Sharrock W., Information Science Reference. Zhao M., Capdeville Y., Zhang H., 2016. Direct numerical modeling of time-reversal acoustic subwavelength focusing, Wave Motion , 67, 102– 115. https://doi.org/10.1016/j.wavemoti.2016.07.010 Google Scholar CrossRef Search ADS   Zijl W., Hendriks M.A., ’t Hart C.M.P., 2002. Numerical homogenization of the rigidity tensor in Hooke’s law using the node-based finite element method, Math. Geol. , 34( 3), 291– 322. https://doi.org/10.1023/A:1014894923280 Google Scholar CrossRef Search ADS   APPENDIX A: THE LOW-PASS FILTER IN THE 1-D CASE To separate the macroscopic and the microscopic scales, the non-periodic homogenization method requires the use of a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$. Applied to any function $$f \! : \! \mathbb {R}\! \rightarrow \! \mathbb {R}$$, this filter removes all the scales smaller than λ0:   $$\mathcal {F}^{\varepsilon _0}\lbrace f\rbrace (x) = \int _\mathbb {R} f(x^{\prime }) w^{\varepsilon _0}(x-x^{\prime }) \, \text{d}x^{\prime }$$ (A1)where $$w^{\varepsilon _0}$$ is the wavelet whose spectrum is   $$\hat{w}^{\varepsilon _0}(k) = \left\lbrace \begin{array}{ll}1 \, {for }k \leqslant \frac{2\pi }{\lambda _0} \\ 0 \, {for} {k} > \frac{2 {\pi}}{\lambda _0}. \end{array} \right.$$ (A2)In this last expression, $$\, \hat{}\,$$ represents the Fourier transform. Note that the filter and the associated wavelet are indexed by ε0 because λ0 depends on the choice of ε0 (λ0 = ε0λm). If defined by spectrum (A2), the wavelet $$w^{\varepsilon _0}$$ has an infinite support, which is unmanageable in practice. For that reason, we introduce a cosine-taper to soften the cut-off:   $$\hat{w}^{\varepsilon _0}(k) = \left\lbrace \begin{array}{lll}1 \, {for }k \leqslant k_0 \\ \frac{1}{2} \left[ 1 + \cos \left( \pi \, \frac{k\,-\,k_0}{\Delta k} \right) \right] \, {for }k \in \left] k_0; \, k_0\!-\!\Delta k \right] \\ 0 \, {for }k > k_0 \end{array} \right.$$ (A3)where $$k_0=\frac{2\pi }{\lambda _0}$$. With such a definition, $$w^{\varepsilon _0}$$ can be truncated with a negligible error to get a compact support (Fig. A1). Figure A1. View largeDownload slide A 1-D spectrum $$\hat{w}^{\varepsilon _0}$$ (left) and its corresponding wavelet $$w^{\varepsilon _0}$$ (right). In this example, the length of the cosine-taper Δk is equal to 0.5 k0. Figure A1. View largeDownload slide A 1-D spectrum $$\hat{w}^{\varepsilon _0}$$ (left) and its corresponding wavelet $$w^{\varepsilon _0}$$ (right). In this example, the length of the cosine-taper Δk is equal to 0.5 k0. Figure A2. View largeDownload slide 3-D spectrum $$\hat{w}_3^{\varepsilon _0}$$ for Δk = 0.5 k0 (left) and its corresponding 3-D wavelet $$w_3^{\varepsilon _0}$$ (right). Because both $$\hat{w}_3^{\varepsilon _0}$$ and $$w_3^{\varepsilon _0}$$ are spherical, we represent them as a function of ‖k‖ and ‖x‖, respectively. Comparing the obtained plots with Fig. A1, we note that the 3-D wavelet is different from just having the wavelet of the 1-D case in all the directions. This is simply because the (inverse) spherical Fourier transform is not the (inverse) 1-D Fourier transform. Figure A2. View largeDownload slide 3-D spectrum $$\hat{w}_3^{\varepsilon _0}$$ for Δk = 0.5 k0 (left) and its corresponding 3-D wavelet $$w_3^{\varepsilon _0}$$ (right). Because both $$\hat{w}_3^{\varepsilon _0}$$ and $$w_3^{\varepsilon _0}$$ are spherical, we represent them as a function of ‖k‖ and ‖x‖, respectively. Comparing the obtained plots with Fig. A1, we note that the 3-D wavelet is different from just having the wavelet of the 1-D case in all the directions. This is simply because the (inverse) spherical Fourier transform is not the (inverse) 1-D Fourier transform. APPENDIX B: $${{\rho} ^{\varepsilon _0}}$$ MUST BELONG TO Ψ Inserting eqs (2)–(4) and properties $$E^{\varepsilon _0}(x,y)$$ and $$\rho ^{\varepsilon _0}(x,y)$$ instead of E(x) and ρ(x) into the initial elastodynamic problem yields a cascade of equations. For i = 0, it turns that   $$\rho ^{\varepsilon _0}(x,y)\ddot{u}_0(x) - \nabla _x\sigma _0(x) - \nabla _y\sigma _1(x,y) = f(x),$$ (B1)where f is the external force and $$\ddot{u}_0$$ is the second time-derivative of the zeroth-order displacement. Because σ1 has to lie in Ψ, $$\rho ^{\varepsilon _0}$$ must belong to Ψ as well. APPENDIX C: THE LOW-PASS FILTER IN THE 3-D CASE $$\mathcal {F}^{\varepsilon _0}_3$$ is the extension of the 1-D low-pass filter $$\mathcal {F}^{\varepsilon _0}$$ (Appendix A) to 3-D. Applied to any function $$f \! : \! \mathbb {R}^3 \! \rightarrow \! \mathbb {R}$$, it removes all the scales smaller than λ0:   $$\mathcal {F}^{\varepsilon _0}_3\lbrace f\rbrace (\boldsymbol {x}) = \int _{\mathbb {R}^3} f(\boldsymbol {x}^{\prime }) w_3^{\varepsilon _0}(\boldsymbol {x}-\boldsymbol {x}^{\prime }) \, \text{d}\boldsymbol {x}^{\prime },$$ (C1)where $$w_3^{\varepsilon _0}$$ is the wavelet whose 3-D spectrum is   \begin{eqnarray} \hat{w}_3^{\varepsilon _0}(\boldsymbol {k}) = \left\lbrace \begin{array}{lll} 1 \, {for }\Vert \boldsymbol {k}\Vert \leqslant k_0 \\ \frac{1}{2} \left[ 1 + \cos \left( \pi \, \frac{\Vert \boldsymbol {k}\Vert \,-\,k_0}{\Delta k} \right) \right] \, {for }\Vert \boldsymbol {k}\Vert \in \left] k_0; \, k_0\!-\!\Delta k \right] \\ 0 \, {for }\Vert \boldsymbol {k}\Vert > k_0. \end{array} \right. \end{eqnarray} (C2) k0 still is the cut-off wave-number $$\frac{2\pi }{\lambda _0}$$. For all the applications within the present paper (Section 4), the length of the cosine-taper Δk has been chosen to be equal to 0.5 k0. When using this value, a truncation of the wavelet at ‖x‖ = 4 λ0 is satisfactory (Fig. A2). APPENDIX D: $${\bf {H}}^{\varepsilon _0}$$ AND $${\chi}^{\varepsilon _0}$$ ARE PERIODIC IN y AND BELONG TO Ψ3 As a solution of a periodic boundary value problem, χs is periodic in y, so ∇yχs is periodic in y. Gs therefore is periodic in y, as well as Hs, consequently. Finally, $${\bf G}^{\varepsilon _0}$$ and $${\bf H}^{\varepsilon _0}$$ are both periodic in y. By construction (eqs 28 and 29), they also belong to Ψ3. Let us denote the cross product by ∧. Because ∇y∧Gs = 0 and $$\boldsymbol {\nabla }_{\!y} \wedge \mathcal {F}^{\varepsilon _0}_3\lbrace \boldsymbol {f}\rbrace = \mathcal {F}^{\varepsilon _0}_3\lbrace \boldsymbol {\nabla }_{\!y} \wedge \boldsymbol {f}\rbrace$$ for any tensorial function f, we have $$\boldsymbol {\nabla }_{\!y} \wedge {\bf G}^{\varepsilon _0} = \boldsymbol {0}$$, so $${\bf G}^{\varepsilon _0}$$ can be written as a constant plus a gradient of a displacement $$\boldsymbol {\chi }^{\varepsilon _0}$$. As a solution of eq. (23), $$\boldsymbol {\nabla }_{\!y} \boldsymbol {\chi }^{\varepsilon _0}$$ has no rotational component, so we end up with   $${\bf G}^{\varepsilon _0} = {\bf I}+\frac{1}{2} (\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}+{^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}).$$ (D1)By construction (eq. 28), $$\left\langle {{\bf G}^{\varepsilon _0}}\right\rangle _3 = {\bf I}$$ (that would not be the case if $${\bf G}^{\varepsilon _0}$$ was built as $$E^{\varepsilon _0}$$ and $$\rho ^{\varepsilon _0}$$ are—eqs 17 and 18–), so   $$\left\langle {\, \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} \right)\,}\right\rangle _3 = 0.$$ (D2)For any tensorial function f ∈ Ψ3 periodic in y,   \begin{eqnarray} \left\langle {\boldsymbol {f}}\right\rangle _3 &=& 0 \ \mbox{and} \ \boldsymbol {\nabla }_{\!y} \boldsymbol {g}\nonumber\\ &=& \boldsymbol {f}\ \Longrightarrow \ \boldsymbol {g} \ \mbox{is periodic in} \ \boldsymbol {y}\ \mbox{and belongs to} \ \Psi _3 , \end{eqnarray} (D3)so $$\boldsymbol {\chi }^{\varepsilon _0}$$ is periodic in y and belongs to Ψ3. This demonstration can also be found in the seminal paper by Capdeville et al. (2010b), section 3.5. APPENDIX E: HEURISTIC FOR SAMPLING A MODEL TO BE HOMOGENIZED When honouring discontinuities, finite element meshes enable a correct account for the geometry of the medium in which the phenomenon to be simulated occurs. Nevertheless, the element size and the interpolation degree which guarantee the accurate capture of the phenomenon cannot be assessed rigorously in most of finite element applications. To choose these two parameters, rules of thumb are often used. In the case of the non-periodic homogenization, we know that our outputs $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$ only contain scales larger than λ0, so we sample them using a space-step $$dx = \frac{\lambda _0}{4}$$. This choice corresponds to twice the Nyquist rate. Because $${\bf C}^{\varepsilon _0\star }$$ is obtained from the space-derivative of the solution χs of the finite element analysis, we required an additional degree for this solution, so we want χs to be sampled by a space-step $$\delta x = \frac{\lambda _0}{5}$$. Such a sampling is ensured by imposing   $$\frac{a}{N} \leqslant \frac{\lambda _0}{5},$$ (E1)where a is the length of the edges of the mesh and N is the interpolation degree. eq. (E1) is our heuristic for sampling a model to be homogenized by our finite element code. In 3-D, assuming the tetrahedral elements to be regular, we can rewrite this heuristic in terms of the volume of the elements v given a degree N:   $$v \leqslant \frac{\sqrt{2}}{12} \left( \frac{N\lambda _0}{5} \right)^3.$$ (E2) © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Non-periodic homogenization of 3-D elastic media for the seismic wave equation

, Volume 213 (2) – May 1, 2018
19 pages

/lp/ou_press/non-periodic-homogenization-of-3-d-elastic-media-for-the-seismic-wave-xgupNZ0vmG
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy032
Publisher site
See Article on Publisher Site

### Abstract

Summary Because seismic waves have a limited frequency spectrum, the velocity structure of the Earth that can be extracted from seismic records has a limited resolution. As a consequence, one obtains smooth images from waveform inversion, although the Earth holds discontinuities and small scales of various natures. Within the last decade, the non-periodic homogenization method shed light on how seismic waves interact with small geological heterogeneities and ‘see’ upscaled properties. This theory enables us to compute long-wave equivalent density and elastic coefficients of any media, with no constraint on the size, the shape and the contrast of the heterogeneities. In particular, the homogenization leads to the apparent, structure-induced anisotropy. In this paper, we implement this method in 3-D and show 3-D tests for the very first time. The non-periodic homogenization relies on an asymptotic expansion of the displacement and the stress involved in the elastic wave equation. Limiting ourselves to the order 0, we show that the practical computation of an upscaled elastic tensor basically requires (i) to solve an elastostatic problem and (ii) to low-pass filter the strain and the stress associated with the obtained solution. The elastostatic problem consists in finding the displacements due to local unit strains acting in all directions within the medium to upscale. This is solved using a parallel, highly optimized finite-element code. As for the filtering, we rely on the finite-element quadrature to perform the convolution in the space domain. We end up with an efficient numerical tool that we apply on various 3-D models to test the accuracy and the benefit of the homogenization. In the case of a finely layered model, our method agrees with results derived from Backus. In a more challenging model composed by a million of small cubes, waveforms computed in the homogenized medium fit reference waveforms very well. Both direct phases and complex diffracted waves are accurately retrieved in the upscaled model, although it is smooth. Finally, our upscaling method is applied to a realistic geological model. The obtained homogenized medium holds structure-induced anisotropy. Moreover, full seismic wavefields in this medium can be simulated with a coarse mesh (no matter what the numerical solver is), which significantly reduces computation costs usually associated with discontinuities and small heterogeneities. These three tests show that the non-periodic homogenization is both accurate and tractable in large 3-D cases, which opens the path to the correct account of the effect of small-scale features on seismic wave propagation for various applications and to a deeper understanding of the apparent anisotropy. Numerical solutions, Computational seismology, Seismic anisotropy, Theoretical seismology, Wave propagation, Wave scattering and diffraction 1 INTRODUCTION Seismic waves are one of the most powerful tools to image the Earth’s interior. Giving access to the geometry of geological structures and to the distribution of mechanical properties within our planet, they lead to a better understanding of geodynamic processes and resource potentials. In the last decades, the seismic tomography and imaging community took advantage of the increasing computational power and the development of efficient numerical methods to improve its techniques and results. Accurate solutions to the forward modelling (see Virieux et al.2011, for a review of the various numerical methods available to model the seismic wave propagation) and the inverse problem (e.g. Pratt et al.1998; Tromp et al.2005; Plessix 2006; Fichtner et al.2008; Métivier et al.2013; Brossier et al.2015; Warner & Guasch 2016) now allow accounting for full seismic waveforms and getting robust and well-resolved models of the Earth at different scales (from the subsurface in exploration geophysics to the entire globe in seismology). One of the important remaining challenges in seismic wave simulation and inversion is the understanding of and the correct account for the effect of small heterogeneities on wave propagation. By small heterogeneities we here mean structures which are smaller than the minimum seismic wavelength propagating in the medium. Seismic waves indeed have a finite frequency band which implies the existence of a minimum wavelength, whereas heterogeneities within the Earth can occur at all scales. When propagating through small heterogeneities, seismic waves naturally average the elastic properties of the medium. A deeper understanding of this averaging process could significantly improve the interpretation of seismic inversion results, pointing out what small-scale features can be ‘hidden’ behind the smooth images one usually obtains. A proper knowledge of the averaging is of major interest for the forward modelling as well, because replacing a given discontinuous and highly heterogeneous model by an effective (or, equivalently, ‘upscaled’ or ‘long-wave equivalent’) medium greatly eases the numerical simulation of wave propagation. When present, small-scale features indeed control the spatial sampling of the model and, consequently, the time sampling too. This is because of a stability condition that all the wave simulators have to ensure (Courant et al.1928):   $$\Delta t \le C \left[ \frac{\Delta x}{V_P} \right]_{\text{min}},$$ (1)Δt denoting the time-step of the simulation, C a constant and $$\left[ \frac{\Delta x}{V_P} \right]_{\text{min}}$$ the minimum ratio of grid spacing Δx and P-wave speed VP. Small heterogeneities therefore induce massive, possibly prohibitive, computation costs. In addition, an accurate simulation usually requires all the physical discontinuities to be honoured by the mesh, which can lead to enormous meshing efforts, especially when dealing with hexahedra (e.g. Casarotti et al.2008; Peter et al.2011). For these reasons, working with effective media is much preferable. Theoretical studies for understanding the effective properties of heterogeneous elastic media was initiated in the sixties by Hashin & Shtrikman (1963) and Hill (1965). These studies both use the pioneering ideas of Voigt (1889) and Reuss (1929) to derive upper and lower-bound effective constants of periodic materials. After these works, many analytical derivations followed (e.g. Kutsenko et al.2013, and the references within), handling more and more complex periodic units or improving accuracy and efficiency of the effective properties calculation. In the case of waves propagating in finely layered media, Backus (1962) derived formula for long-wave equivalent elastic coefficients which are still widely used in seismic exploration. Fruitfully, these formulae quantify the seismic anisotropy produced by fine layering, showing that a stack of isotropic layers can explain the anisotropy observed on wave measurements. Schoenberg & Muir (1989) extended the Backus theory by including any kind of anisotropy within each layer, which enables the account for several sets of fractures in the layers. Further studies on the upscaling of fractured media then came out, such as Sayers & Kachanov (1991); Mauge & Kachanov (1994); Sayers & Kachanov (1995); Schoenberg & Sayers (1995); Schoenberg & Helbig (1997); Sayers (1998); Grechka & Kachanov (2006); Grechka (2007); Carcione et al. (2012). At larger scale, particular efforts have been focused on smoothing the Earth’s crust to ease the simulation of long-period surface waves (Capdeville & Marigo 2008; Fichtner & Igel 2008; Lekić et al.2010). More recently, the understanding of the seismic anisotropy in the upper-mantle drew the attention (Fichtner et al.2013a; Wang et al.2013; Bodin et al.2015). In this context, Jordan (2015) developed an elegant effective medium theory for random media in which statistics of the local anisotropy (produced by lattice-preferred orientation of mineral grains for instance) are separated from those of ellipsoidal geometric heterogeneities (i.e. shape-preferred orientation of simple geological structures). In this paper, the upscaling method in consideration is the two-scale homogenization technique. This method emerged in the seventies from research in micromechanics for predicting the macroscopic response of composite and random materials to either static or dynamic excitations (Bensoussan et al.1978; Sanchez-Palencia 1980; Papanicolaou & Varadhan 1981). Since then, the technique has been successfully applied to many physical processes, such as heat transfer (e.g. Allaire & Habibi 2013), Stokes flow (e.g. Hornung 1997), neutronic diffusion (e.g. Allaire & Capdeboscq 2000), magnetization (e.g. Santugini-Repiquet 2007) and elastic wave propagation (e.g. Boutin & Auriault 1993; Fish & Chen 2001, 2004; Parnell & Abrahams 2008; Bacigalupo & Gambarotta 2014). In this last field, the two-scale homogenization was adapted to non-periodic media within the last decade (Capdeville & Marigo 2007; Capdeville et al.2010a,b; Guillot et al.2010; Cance & Capdeville 2015; Capdeville et al.2015), which opened the path to the upscaling of general elastic media, with no constraint on the shape and size of the heterogeneities. The method has been tested in 1-D (Capdeville & Marigo 2007; Capdeville et al.2010a) and 2-D (Capdeville et al.2010b; Guillot et al.2010; Capdeville et al.2015); the goal of the present paper is to show its efficiency (in terms of accuracy and computation speed) on 3-D cases. To solve the homogenization problem, we will rely on a finite-element method, but it is worth noting here that the algorithm proposed by Capdeville et al. (2015) could be used as well. This latter is up and running on continuous or pixel-based 3-D media. The non-periodic two-scale homogenization actually is different to what some authors called numerical homogenization (e.g. Zijl et al.2002; Grechka 2003; Gao et al.2015) or heterogeneous multi-scale methods (e.g. Engquist et al.2009; Abdulle & Grote 2011). In these classes of techniques, the effective elastic tensor is computed by solving the wave equation in the static regime (i.e. at zero frequency so that a simple elastostatic problem arises) with a set of unit stresses applied on a representative volume. Such a volume implicitly delineates the small and the large scales. It fully depends on the mesh and the numerical technique that will then be used for simulating seismic waves. Even though our non-periodic homogenization method also includes the numerical resolution of a static problem, this latter does not involve a representative volume. Our static problem is actually defined on the whole medium at once with volume forces equal to the divergence of the elastic tensor (eq. 33). Low-pass filtering the strain and the stress which result from these forces then leads to the effective properties (eq. 30). Our filter only depends on the minimum wavelength to be propagated in the medium and a precision factor, so that the scales are separated by the waves rather than imposed by a numerical method. This makes non-periodic homogenization results physically meaningful and independent from wave equation solvers. Homogenized properties actually reveal what the waves ‘see’. In particular, Capdeville et al. (2013) showed that models obtained by full-waveform inversion are homogenized models. Applying the non-periodic homogenization method in 3-D does not require further theoretical derivations than those developed in 2-D by Capdeville et al. (2010b) and Guillot et al. (2010). In a first part, the present paper recalls the ideas and concepts of the method. From this theoretical part emerge the practical issues which have to be tackled to get homogenized properties, namely (i) the resolution of an elastostatic problem and (ii) a filtering process. The second part of the article describes the implementation of a parallel finite-element method for solving the 3-D elastostatic problem and provides some details on the 3-D filtering. In a third part, the accuracy and performance of the resulting homogenization code are challenged on various elastic models: (i) a finely layered medium for which analytical expressions of the upscaled properties are available, (ii) a large and highly heterogeneous medium in which reference seismograms can be computed and (iii) a realistic geological model made of multiple folded and faulted horizons. Finally, we discuss the possible improvements of our code and the numerous perspectives opened by the 3-D homogenization method in terms of solving forward and inverse problems. 2 THE NON-PERIODIC HOMOGENIZATION THEORY IN A NUTSHELL We here attempt to give a synoptic overview of the non-periodic homogenization theory, starting from basics introduced in micromechanics for 1-D periodic materials (e.g. Bensoussan et al.1978), then extending these basics to non-periodic media (Capdeville et al.2010a), and finally moving to the general 2-D/3-D non-periodic case (Capdeville et al.2010b; Guillot et al.2010). Our goal is to provide a digest of the theory to make the non-periodic homogenization understandable by a large number of geophysicists, so we intentionally skip some technical details of the whole Capdeville et al. (2010a,b) derivation to focus on the main ideas and concepts of the method. 2.1 1-D periodic homogenization The homogenization theory relies on an ansatz for the solution of the physical problem in consideration. In elastodynamics, the displacement u(x; t) and the stress σ(x; t) involved in the 1-D case (x being the space variable and t being the time) are postulated to be two-scale asymptotic expansions:   $$u(x;t) = \sum _{i=0}^{+\infty } \varepsilon ^i u_i\left(x,\frac{x}{\varepsilon };t\right),$$ (2)  $$\sigma (x;t) = \sum _{i=0}^{+\infty } \varepsilon ^i \sigma _i\left(x,\frac{x}{\varepsilon };t\right),$$ (3)where $$\varepsilon =\frac{l}{\lambda _m}$$ is the ratio of the size of the periodic cell which constitutes the 1-D medium to the minimum wavelength propagating in this medium (Fig. 1). By definition, l is microscopic and λm is macroscopic. ε, which is called scaling parameter, is therefore much smaller than 1. It enables to explicitly separate the scales within coefficients ui and σi of series (2) and (3), x capturing the large-scale variations and $$y=\frac{x}{\varepsilon }$$ handling the small-scale variations. Because any change in y induces a very slight change in x, the two variables can be treated independently and the spatial derivative operator can be written   $$\nabla = \varepsilon ^{-1}\nabla _y + \nabla _x.$$ (4) Figure 1. View largeDownload slide 1-D periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in an infinite 1-D medium made of periodic cells whose size l is much smaller than λm. Figure 1. View largeDownload slide 1-D periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in an infinite 1-D medium made of periodic cells whose size l is much smaller than λm. The introduction of the small-scale variable $$y=\frac{x}{\varepsilon }$$ also allows us to rewrite the physical parameters E(x) (Young modulus) and ρ(x) (density) of the medium as λm-periodic quantities E(y) and ρ(y). In other words, the l-periodic bar in the initial one-variable problem is now seen as a medium containing only small scales which are repeated in space with a λm-periodicity. Furthermore, coefficients ui and σi are assumed to be λm-periodic in y as are E and ρ. This assumption imposes that small-scale variations of the displacement and stress fields are due to local small-scale structures of the medium. In seismology, such a phenomenon is commonly called site effect. Plugging eqs (2)–(4) into the elastodynamic problem (i.e. the wave equation and Hooke’s law) yields a cascade of coupled equations which can be solved for each i using the average over the periodic cell   $$\langle f \rangle (x) = \frac{1}{\lambda _m}\int _{-\frac{\lambda _m}{2}}^{\frac{\lambda _m}{2}} f(x,y) \, \text{d}y, \, \forall f \! : \! \mathbb {R}^2 \! \rightarrow \! \mathbb {R},$$ (5)and the periodicity in y of the problem. Doing so, it turns out that the zeroth-order terms u0 and σ0 do not depend on the microscopic variable y. This result formalizes the poor sensitivity of the wavefield to small heterogeneities. Going further, one shows that u0 and σ0 are the solution of the so-called homogenized problem, which is a classical elastodynamic problem with homogeneous effective properties E⋆ and ρ⋆ such that   $$\rho ^{\star } = \langle{\rho (y)}\rangle$$ (6)and   $$E^{\star } = \left\langle {E(y)[1+\nabla _y\chi (y)]}\right\rangle .$$ (7)In this last equation, χ is the so-called first-order corrector. It is the solution of the so-called cell problem, which is   $$\nabla _y\left[E(1+\nabla _y\chi )\right] = 0$$ (8)defined on the cell with periodic boundary conditions. Because an analytical solution exists for this problem, eq. (7) finally reduces to   $$E^{\star } = \left\langle {\,\frac{1}{E(y)}\,}\right\rangle ^{-1}.$$ (9)In conclusion, the quantity to average for obtaining the zeroth-order long-wave equivalent Young modulus E⋆ is the inverse of the initial Young modulus E (eq. 9), whereas the zeroth-order long-wave equivalent density ρ⋆ simply is the average of the initial density ρ (eq. 6). Contrary to u0 and σ0, the higher-order terms of series (2) and (3) contain small-scale variations. For instance, the first-order displacement term is y-dependent through the first-order corrector:   $$u_1(x,y) = \chi (y) \nabla _xu_0(x) + \langle {u_1}\rangle (x).$$ (10)Examples of the contribution of this last term to the whole wavefield can be found in Capdeville et al. (2010a). In the present work, we focus on the order 0 so we do not provide any further details on the higher-order terms. A crucial point in the homogenization theory is that the asymptotic convergence can be proved mathematically. One actually does not directly show the convergence of series (2) and (3) towards the exact wavefield; one uses the so-called Γ-convergence instead (Dal Maso 1993). Rather than studying a single problem for the physically relevant value of ε, one considers a sequence of problems indexed by ε which is now regarded as a small parameter going to zero. In other words, one builds a fictitious sequence of problems in which the periodic cell becomes smaller and smaller. In this context, one demonstrates that the exact solution converges to the solution of the homogenized problem u0 and σ0 (Nguetseng 1989; Allaire 1992). Such a demonstration provides rigorous mathematical foundations for the homogenization theory. 2.2 1-D non-periodic homogenization Let us now consider any given distribution of the properties E and ρ within a finite 1-D medium. The length of the medium is denoted by L. Again, our goal here is to find long-wave equivalent properties. We denote by λm the minimum wavelength of the wavefield propagating in the bar (Fig. 2). Inspired by the periodic homogenization, we will represent the displacement u(x; t) and the stress σ(x; t) as power series in ε (to be redefined) with coefficients ui and σi as functions of x and $$y=\frac{x}{\varepsilon }$$. Plugging these series into the wave equation and Hooke’s law will yield a cascade of equations. To get a solution for these equations, let us bring a periodicity into the problem by imposing periodic boundary conditions at the border of the domain. We end up with an infinite medium made of a periodic cell, so we can certainly benefit from the mathematical results of the previous section (eqs 6 and 9). Nevertheless, two big issues remain: Because the boundary conditions are forced to be periodic, the effective properties computed near the border of the domain are not meaningful for any kind of physically interesting conditions such as Dirichlet or Neumann conditions. In the present paper, we do not investigate this issue. Disregarding boundary effects, we focus on the computation of accurate effective properties in the interior of the domain. Precisions on what we exactly mean by the interior of the domain are given in Part 3. Contrary to the periodic case, the size of the periodic cell is not microscopic. The cell can contain various sizes of heterogeneity, including macroscopic scales (i.e. sizes equal to or larger than λm), so ε defined as the ratio of the size L of the periodic cell to the minimum wavelength λm no longer is a relevant scaling parameter to separate the scales through variables x and $$y=\frac{x}{\varepsilon }$$. The goal of the present section is (i) to redefine ε and (ii) to separate the scales within the mechanical properties of the bar in a way that allows for a solution of our homogenization problem. Figure 2. View largeDownload slide 1-D non-periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in a finite 1-D medium of size L. Contrary to the periodic case, no assumption is made on the distribution of the mechanical properties E(x) and ρ(x) within the medium. Figure 2. View largeDownload slide 1-D non-periodic homogenization framework: a wavefield having a minimum wavelength λm propagates in a finite 1-D medium of size L. Contrary to the periodic case, no assumption is made on the distribution of the mechanical properties E(x) and ρ(x) within the medium. 2.2.1 Redefining ε Because there is no quantity for defining the small scales yet, we introduce λ0 < λm: all the heterogeneities whose size is smaller than λ0 are considered as small. λ0 therefore is analogous to l in the periodic case. Using it, we can bring in a new scaling parameter $$\varepsilon _0=\frac{\lambda _0}{\lambda _m} < 1$$ and we can choose ε ≤ ε0. Similarly to the periodic case, ε is a formal quantity which can go to zero to prove the Γ-convergence of the asymptotic solution. It can be seen as the ratio of a small length λ ≤ λ0 to the minimum wavelength of the wavefield to be propagated λm, λ becoming smaller and smaller when studying the convergence. In practice, the only physically relevant value of ε is ε0. With this definition of ε, variables x and $$y=\frac{x}{\varepsilon }$$ can be treated independently. We denote by Ψ the set of functions $$f(x,y) \! : \mathbb {R}^2 \! \rightarrow \! \mathbb {R}$$ such that the x part of f carries the large-scale variations while the y part of f handles the small-scale variations. As in the periodic case, we impose coefficients ui and σi of the asymptotic expansions (2) and (3) to belong to Ψ. Moreover, we can write the spatial derivative operator ∇ as in eq. (4). 2.2.2 Separating the scales within the 1-D medium In the non-periodic case, the main issue we face is that the medium contains both microscopic and macroscopic scales, so we cannot write E and ρ as a function of y only. Using a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$ (Appendix A) to separate large-scale and small-scale variations, we will have to find a proper way to build $$E^{\varepsilon _0}(x,y)$$ and $$\rho ^{\varepsilon _0}(x,y)$$ from E(x) and ρ(x). By proper we here mean which allows for a solution to the cascade of equations that arise when plugging eqs (2)–(4) into the elastodynamic problem. Let us introduce $$\eta =\frac{L}{\lambda }$$ (which implies that $$\varepsilon =\frac{L}{\eta \lambda _m}$$) and assume that we properly built $$E^{\varepsilon _0}(x,y)$$ and $$\rho ^{\varepsilon _0}(x,y)$$ from E(x) and ρ(x). Because E and ρ are L-periodic in x, $$E^{\varepsilon _0}$$ and $$\rho ^{\varepsilon _0}$$ are ηλm-periodic in y. Assuming that coefficients ui(x, y) and σi(x, y) are also ηλm-periodic in y and using the average over the periodic cell   $$\langle{f}\rangle (x) = \frac{1}{\eta \lambda _m}\int _{-\frac{\eta \lambda _m}{2}}^{\frac{\eta \lambda _m}{2}} \! f(x,y) \, \text{d}y , \forall f \! : \! \mathbb {R}^2 \! \rightarrow \! \mathbb {R}\, ,$$ (11)we can easily derive the cascade of equations. Similarly to the periodic case, it turns out that the zeroth-order displacement u0 and stress σ0 do not depend on y. Again, these fields are the solution of the so-called homogenized problem, which is a classical elastodynamic problem involving effective properties $$\rho ^{\varepsilon _0\star }$$ and $$E^{\varepsilon _0\star }$$ such that   $$\rho ^{\varepsilon _0\star }(x) = \langle {\rho ^{\varepsilon _0}(x,y)}\rangle$$ (12)and   $$E^{\varepsilon _0\star }(x) = \left\langle {E^{\varepsilon _0}(x,y)[1+\nabla _y\chi ^{\varepsilon _0}(x,y)]}\right\rangle .$$ (13)Note that the effective properties no longer are homogeneous in the non-periodic case. Moreover, they depend on the choice of ε0. The first-order corrector $$\chi ^{\varepsilon _0}(x,y)$$ also depends on ε0. It necessarily belongs to Ψ and it is ηλm-periodic in y. Similarly to the periodic case, it is the solution of the cell problem, that is,   $$\nabla _y\left[E^{\varepsilon _0}(1+\nabla _y\chi ^{\varepsilon _0})\right] = 0$$ (14)with periodic boundary conditions. As shown by Capdeville et al. (2010a), function $$\nabla _y\chi ^{\varepsilon _0}(x,y)$$ therefore satisfies   $$\nabla _y\chi ^{\varepsilon _0} = -1 + \left\langle {\,\frac{1}{E^{\varepsilon _0}}\,}\right\rangle ^{-1} \frac{1}{E^{\varepsilon _0}}.$$ (15)This last equation implies that Eq. (13) reduces to   $$E^{\varepsilon _0\star }(x) = \left\langle {\,\frac{1}{E^{\varepsilon _0}(x,y)}\,}\right\rangle ^{-1},$$ (16)which along with eq. (12) tells the quantity to average to get the zeroth-order effective properties. $$\frac{1}{E^{\varepsilon _0}(x,y)}$$ belongs to Ψ, which is the condition to meet for building $$E^{\varepsilon _0}(x,y)$$ properly. Another equation from our cascade tells that $$\rho ^{\varepsilon _0}(x,y)$$ has to lie in Ψ (Appendix B), so we also have a condition for building $$\rho ^{\varepsilon _0}(x,y)$$ properly. Following these conditions, we easily form   $$E^{\varepsilon _0}(x,y) \, = \left[ \mathcal {F}^{\varepsilon _0}\left\lbrace \frac{1}{E}\right\rbrace \!(x) \, + \left[\frac{1}{E}-\mathcal {F}^{\varepsilon _0}\left\lbrace \frac{1}{E}\right\rbrace \right]\!(y) \right]^{-1}$$ (17)and   $$\rho ^{\varepsilon _0}(x,y) \, = \, \mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x) \, + \, [\rho -\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace ](y).$$ (18)Fig. 3 illustrates the mathematical construction. In this figure, symbol g represents either ρ or $$\frac{1}{E}$$, and FT stands for Fourier transform. From any distribution of g over a length L extended to $$\mathbb {R}$$ by periodicity (Fig. 3a), large scales (Fig. 3b) and small scales (Fig. 3c) are extracted using $$\mathcal {F}^{\varepsilon _0}$$. Within the small scales, space variable x is replaced by $$y=\frac{x}{\varepsilon }$$ (Fig. 3d), thus changing the L-periodicity in x to a ηλm-periodicity in y. Then, the two-variable quantity $$g^{\varepsilon _0}(x,y)$$ is formed by simply adding the large scales (which are a function of x) to the small scales (which are a function of y). Fig. 3(e) represents $$g^{\varepsilon _0}$$ for a particular value of $$x= \bar{x}$$; it shows the small scales oscillating around $$\mathcal {F}^{\varepsilon _0}\lbrace g\rbrace ( \bar{x})$$. Figure 3. View largeDownload slide Separation of the scales within a 1-D medium of length L. g here represents either ρ or $$\frac{1}{E}$$. From any distribution of g extended to $$\mathbb {R}$$ by periodicity (top row), we build $$g^{\varepsilon _0}(x,y)$$ (bottom row) using a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$. In this figure, $$g^{\varepsilon _0}(x,y)$$ is represented for a given $$x= \bar{x}$$. See the text for more details. Figure 3. View largeDownload slide Separation of the scales within a 1-D medium of length L. g here represents either ρ or $$\frac{1}{E}$$. From any distribution of g extended to $$\mathbb {R}$$ by periodicity (top row), we build $$g^{\varepsilon _0}(x,y)$$ (bottom row) using a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$. In this figure, $$g^{\varepsilon _0}(x,y)$$ is represented for a given $$x= \bar{x}$$. See the text for more details. 2.2.3 Final results Including eq. (18) into (12) yields   \begin{eqnarray} \rho ^{\varepsilon _0\star }(x) & =& \left\langle {\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x) + [\rho -\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace ](y)}\right\rangle \nonumber \\ &=& \left\langle {\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x)}\right\rangle + \langle {[\rho -\mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace ](y)}\rangle \nonumber \\ & =& \mathcal {F}^{\varepsilon _0}\lbrace \rho \rbrace (x). \end{eqnarray} (19) Similarly, including eq. (17) into (16) leads to   $$E^{\varepsilon _0\star }(x) = \left[ \mathcal {F}^{\varepsilon _0}\left\lbrace \frac{1}{E}\right\rbrace (x) \right]^{-1}.$$ (20)Eqs (19) and (20) tell that the zeroth-order effective properties of any 1-D medium can be computed by filtering the density and the inverse of the Young modulus. This result can be easily intuited from the solution inferred in the periodic case (eqs 6 and 9). We have here given a rigorous demonstration for it. To derive a zeroth-order solution for our non-periodic homogenization problem, we introduced a new scaling parameter ε0. Because the obtained effective medium depends on this parameter (eqs 19 and 20), u0 and σ0 also depend on it. For sake of simplicity, we did not index these two fields by ε0 as Capdeville et al. (2010a) did. Nevertheless, we will study the ε0-convergence of these fields in a specific case (Section 4.2). 2.3 3-D non-periodic homogenization In the 3-D case, Hooke’s law involves a fourth-order tensor  C to relate the stress and the strain, so the linear elastic behaviour of a given medium no longer is fully described by E only. Nevertheless, the homogenization theory developed for 1-D non-periodic media is still valid up to the zeroth-order effective properties (12) and (13) which are now written   $$\rho ^{\varepsilon _0\star }(\boldsymbol {x}) = \left\langle {\rho ^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})}\right\rangle _3$$ (21)and   \begin{eqnarray} {\bf C}^{\varepsilon _0\star }(\boldsymbol {x}) = \left\langle {\, {\bf C}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) : \left[ {\bf I}+ \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) \right) \right] \,}\right\rangle _3\!,\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (22)where I is the fourth-order identity tensor, t is the transpose operator, : is the tensor contraction [A:  B]ijkl = AijmnBmnkl and 〈〉3 is the average on y of any function $$f(\boldsymbol {x},\boldsymbol {y}) \! : \! \mathbb {R}^{3\times 2} \! \rightarrow \! \mathbb {R}$$ over the periodic cell. The first-order corrector $$\boldsymbol {\chi }^{\varepsilon _0}$$ now is a third-order tensor. It is periodic in y and it necessarily belongs to Ψ3 (the extension of Ψ to 3-D). It is the solution of the following cell problem:   $$\boldsymbol {\nabla }_y \cdot \left\lbrace {\bf C}^{\varepsilon _0} : \left[ {\bf I}+ \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} \right) \right] \right\rbrace = \boldsymbol {0}$$ (23)with periodic boundary conditions. Contrary to the 1-D case, there is no analytical solution for this problem here (unless the medium is layered transverse isotropic, as demonstrated by Guillot et al. (2010) and Lin et al. (2017)), so we are not able (i) to write $${\bf C}^{\varepsilon _0\star }$$ as a function of $${\bf C}^{\varepsilon _0}$$ only and (ii) to get a mathematical condition on $${\bf C}^{\varepsilon _0}$$ for properly separating the scales within the elastic properties of the medium. To overcome this issue, a procedure inspired by Papanicolaou & Varadhan (1981) is proposed. Introducing $${\bf G}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})$$ and $${\bf H}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})$$ such that   \begin{eqnarray} {\bf H}^{\varepsilon _0} = {\bf C}^{\varepsilon _0}:{\bf G}^{\varepsilon _0} \end{eqnarray} (24)  \begin{eqnarray} \phantom{{\bf H}^{\varepsilon _0}}= {\bf C}^{\varepsilon _0} : \left[{\bf I}+\frac{1}{2} (\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}+{^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0})\right], \end{eqnarray} (25)this procedure enables an implicit construction of $${\bf C}^{\varepsilon _0\star }$$ such that both $${\bf H}^{\varepsilon _0}$$ and $$\boldsymbol {\chi }^{\varepsilon _0}$$ are periodic in y and belong to Ψ3, which are the conditions to meet to get a proper solution to our homogenization problem. Here are the steps of the procedure: Eq. (23) with periodic boundary conditions is solved using C(y) instead of $${\bf C}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})$$. C(y) is just the original elastic tensor in which the space variable x has been changed by y = ε0x, meaning that all the scales within the medium are considered as small in this first step. This yields a cell problem that we can solve numerically (Part 3). The solution of such a cell problem is called the starting corrector χs(y). It can be seen as the static response of the medium to local unit strains expressed by the identity tensor I. From χs, two tensors Gs and Hs are built:   $${\bf G}_s(\boldsymbol {y}) = {\bf I}+ \frac{1}{2} \left[ \boldsymbol {\nabla }_y\boldsymbol {\chi }_s(\boldsymbol {y}) + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }_s(\boldsymbol {y}) \right],$$ (26)  $${\bf H}_s(\boldsymbol {y}) = {\bf C}(\boldsymbol {y}) : {\bf G}_s(\boldsymbol {y}).$$ (27)Gs can be seen as the unit strains plus the strains associated with the starting corrector. It is called strain concentrator. In the same way, Hs is called stress concentrator. Using $$\mathcal {F}^{\varepsilon _0}_3$$ (the extension of $$\mathcal {F}^{\varepsilon _0}$$ to 3-D, see Appendix C), the large-scale and small-scale variations are separated in Gs and Hs to form $${\bf G}^{\varepsilon _0}$$ and $${\bf H}^{\varepsilon _0}$$:   $${\bf G}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) = {\bf I}+ \left[{\bf G}_s - \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace \right ](\boldsymbol {y}) : \left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace (\boldsymbol {x}) \right]^{-1} ,$$ (28)  \begin{eqnarray} {{\bf H}^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y})}\nonumber\\ {\quad = \left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf H}_s\rbrace (\boldsymbol {x}) + [{\bf H}_s - \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf H}_s\rbrace ](\boldsymbol {y}) \right] :\left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace (\boldsymbol {x}) \right]^{-1}.} \end{eqnarray} (29)Through definitions (24) and (25), $${\bf G}^{\varepsilon _0}$$ and $${\bf H}^{\varepsilon _0}$$ are the strain and the stress concentrators associated with the first-order corrector $$\boldsymbol {\chi }^{\varepsilon _0}$$. They are here constructed by separating the scales in Gs and Hs (eqs 28 and 29). From such a construction, one can demonstrate that $${\bf H}^{\varepsilon _0}$$ and $$\boldsymbol {\chi }^{\varepsilon _0}$$ are indeed periodic in y and belong to Ψ3 (Appendix D). Using (25) in (22), we note that $${\bf C}^{\varepsilon _0\star } = \left\langle {{\bf H}^{\varepsilon _0}}\right\rangle_3$$. Introducing (29) into this latter equation, we end up with   $${\bf C}^{\varepsilon _0\star }(\boldsymbol {x}) = \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf H}_s\rbrace (\boldsymbol {x}) : \left[ \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf G}_s\rbrace (\boldsymbol {x}) \right]^{-1}.$$ (30) The zeroth-order effective elastic tensor given by the procedure therefore is obtained by filtering the stress and the strain associated with the starting corrector. This means that one just needs to go through steps (i), (ii) and (iv) to get this tensor in practice. Step (iii) would be necessary to obtain $$\boldsymbol {\chi }^{\varepsilon _0}$$ (through the construction of $${\bf C}^{\varepsilon _0} = {\bf H}^{\varepsilon _0} : {\bf G}^{\varepsilon _0^{\ -1}}$$ and the resolution of eq. 23 for instance) in order to access the first-order displacement. Because we focus on the zeroth-order solution in the present paper, we will not deal with this latter aspect. Furthermore, there is no demonstration for the minor and major symmetries of $${\bf C}^{\varepsilon _0\star }$$ yet. We do not address this point here, but we have checked the symmetries (up to a certain precision) in all the applications that we present in Part 4. As regards the density, an equation similar to (B1) also emerges in 3-D, implying that $$\rho ^{\varepsilon _0}$$ must lie in Ψ3, so the scales can be separated as they were in 1-D:   $$\rho ^{\varepsilon _0}(\boldsymbol {x},\boldsymbol {y}) \, = \, \mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace (\boldsymbol {x}) \, + \, \left[\rho -\mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace \right ](\boldsymbol {y}).$$ (31)Using (31) in (21), the effective density comes out:   $$\rho ^{\varepsilon _0\star } = \mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace.$$ (32)As in the 1-D case, $$\rho ^{\varepsilon _0\star }$$ is obtained by simply filtering the initial density ρ. 3 IMPLEMENTATION OF THE 3-D NON-PERIODIC HOMOGENIZATION While the non-periodic homogenization theory involves many specific concepts and quantities, its zeroth-order result can be stated quite shortly. From a given 3-D medium described by its density ρ and its elastic tensor C, a long-wave equivalent density $$\rho ^{\varepsilon _0\star }$$ can be computed by just filtering ρ (eq. 32), and a long-wave elastic tensor $${\bf C}^{\varepsilon _0\star }$$ can be obtained by Calculate the starting corrector χs from the cell problem (23) using C with periodic boundary conditions instead of $${\bf C}^{\varepsilon _0}$$, that is,   $$\boldsymbol {\nabla }_y \cdot \left\lbrace {\bf C}: \left[ {\bf I}+ \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }_s + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }_s \right) \right] \right\rbrace = \boldsymbol {0}.$$ (33) Build Gs and Hs from χs using eqs (26) and (27), and then filter these two quantities to get $${\bf C}^{\varepsilon _0\star }$$ (eq. 30). Differential eq. (33) can be seen as a classic elastostatic problem with a specific load consisting in the divergence of the elastic tensor C. Such a divergence yields a third-order tensor ∂iCijkl. Thanks to the symmetry of C (Cijkl = Cijlk), this tensor reduces to six force vectors. To determine the six corresponding displacements, eq. (33) can be solved in two different ways: using a strong-form iterative scheme based on Fast Fourier Transform (Moulinec & Suquet 1998; Capdeville et al.2015) or using a more classic weak-form finite-element approach (e.g. Hughes 2012). The latter is adapted to strongly discontinuous media. It is the method we implement here, relying on tetrahedral meshes. Because we want to investigate the behaviour of the solution with respect to the discretization, our code allows polynomial interpolations of degree 1, 2 or 3 (Worth et al.2012) in addition with various quadrature rules (Felippa 2004). Moreover, both linear and quadratic tetrahedra are enabled, leading to either iso-, super- or subparametric elements. The obtained linear system involving six right-hand side members, we solve it using a direct solver. Among several codes, PARDISO (Schenk & Gärtner 2006) is chosen here. Finally, the low-pass filter $$\mathcal {F}^{\varepsilon _0}_3$$ is applied in the space domain to obtain $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$. To perform the 3-D convolution, we use the quadrature associated with the mesh employed in finite element analysis. The periodicity imposed at the boundary ∂Ω of the elastic domain Ω involved in eq. (33) means that the medium is supposed to repeat itself periodically in the 3-D. When dealing with geological media, this condition is obviously not fulfilled. It is therefore replaced by either a homogeneous Neumann condition or a homogeneous Dirichlet condition. We choose the second option in our implementation so we impose χs = 0 on ∂Ω. The effect in the volume Ω of such an artificial condition decays exponentially (Dumontet 1990), so our numerical solution is meaningless in a thin layer from the border of the domain. Such meaningless values of χs do not matter, because the filtering process cannot be performed near ∂Ω anyway. Some elastic material to be convolved with the wavelet is actually missing there, so we are not able to compute the effective properties $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$ using (32) and (30). The layer ΩO in which the filter cannot be applied is called the outer domain. Its thickness is equal to half of the wavelet support. Its complement ΩI = Ω − ΩO is called the inner domain (or the interior of the domain). The solutions proposed in the present paper make sense in ΩI only. Further developments, such as those initiated by Capdeville & Marigo (2008, 2013), would be necessary to get meaningful effective properties in ΩO. When handling large models, the memory requirements for achieving the computation of the effective properties can be very large. This is mainly because solving large linear systems, even symmetric, is memory-demanding. This is also because high-order tensors are involved in the homogenization process. For these two reasons, a distributed-memory computation is necessary. A classic way of implementing it is presented in Fig. 4. It consists in Decomposing the domain Ω into n subdomains Ωk such that $$\Omega = \bigcup _{k=1}^n \Omega _k$$. Each subdomain being handled by a processor, n local stiffness matrices Kk are computed. Assembling the local stiffness matrices to build the global stiffness matrix K. Solving the linear system using a parallel solver to get the finite element solution χs of eq. (33). From χs, Gs is formed following (26). Using the domain decomposition again to filter Gs, Hs = C: Gs and ρ in each subdomain. From this filtering, the effective properties $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ are derived. In the outer domain ΩO (the dotted grey lines in Fig. 4), these properties cannot be computed because the filter cannot be applied. The same problem appears near the frontiers between the subdomains, so we have to either implement massive communications between the subdomains or enlarge each subdomain Ωk by a buffer layer $$\Omega _k^B$$ which is as thick as half of the filtering wavelet support. This second option is sketched in Fig. 4. It yields n overlapping parts$$P_k = \Omega _k \cup \Omega _k^B$$. Along with possible pieces of ΩO in Ωk, the buffer layer $$\Omega _k^B$$ acts as an outer part$$P_k^O = (\Omega ^O \cap \Omega _k) \cup \Omega _k^B$$. The effective properties can be calculated in the inner part$$P_k^I = P_k - P_k^O$$, so in the whole inner domain because $$\bigcup _{k=1}^n P_k^I = \Omega ^I$$. Communicating between processors to gather the results from the n subdomains. Figure 4. View largeDownload slide A classic distributed-memory workflow for solving the non-periodic homogenization problem. (a) The domain Ω is decomposed in n subdomains Ωk. As an example here n = 3. In each subdomain, a local stiffness matrix Kk is computed. (b) The local stiffness matrices are assembled to obtain the global stiffness matrix K. (c) The linear system is inverted using a parallel solver. From the gradient of the solution, Gs is formed. (d) Using the domain decomposition again, the effective properties $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ are computed in each subdomain by filtering Gs, Hs = C: Gs and ρ. In the outer domain (dotted grey), the filtering wavelet $$w^{\varepsilon _0}$$ cannot be applied so $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ cannot be computed. The same problem appears near the frontiers between the subdomains, so we enlarge these latter to build n overlapping parts Pk equipped with an outer part $$P_k^O$$ and an inner part $$P_k^I$$ in which the filter can be applied. (e) The result in the whole inner domain ΩI is obtained by merging the results from the n inner parts. Figure 4. View largeDownload slide A classic distributed-memory workflow for solving the non-periodic homogenization problem. (a) The domain Ω is decomposed in n subdomains Ωk. As an example here n = 3. In each subdomain, a local stiffness matrix Kk is computed. (b) The local stiffness matrices are assembled to obtain the global stiffness matrix K. (c) The linear system is inverted using a parallel solver. From the gradient of the solution, Gs is formed. (d) Using the domain decomposition again, the effective properties $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ are computed in each subdomain by filtering Gs, Hs = C: Gs and ρ. In the outer domain (dotted grey), the filtering wavelet $$w^{\varepsilon _0}$$ cannot be applied so $${\bf C}^{\varepsilon _0\star }$$ and $$\rho ^{\varepsilon _0\star }$$ cannot be computed. The same problem appears near the frontiers between the subdomains, so we enlarge these latter to build n overlapping parts Pk equipped with an outer part $$P_k^O$$ and an inner part $$P_k^I$$ in which the filter can be applied. (e) The result in the whole inner domain ΩI is obtained by merging the results from the n inner parts. Such a workflow would certainly work, but it imposes the resolution of the whole linear system to fit the distributed-memory at once, meaning that a large stiffness matrix requires a large parallel computer. To overcome this limitation, we propose an alternative implementation which consists in working on the parts Pk all along the workflow (Fig. 5). Applying a homogeneous Dirichlet condition χs = 0 at the boundary of the parts, each of them becomes like a new domain which is totally independent from the other parts. Thanks to the outer parts $$P_k^O$$, we get the effective properties in all the inner parts $$P_k^I$$, so in the whole inner domain $$\Omega ^I = \bigcup _{k=1}^n P_k^I$$. This alternative implementation is embarrassingly parallel because the parts are independent from each other all along the workflow, which means that they can be treated sequentially. This allows for the homogenization of large models on small computers, possibly laptops, provided the available RAM covers the memory requirement of every single part. Such a distributed-memory implementation is also proposed by Capdeville et al. (2015). Figure 5. View largeDownload slide An alternative way of solving the non-periodic homogenization in parallel. (a) The domain Ω is decomposed in n parts Pk. As an example here n = 3. The parts overlap to equip them all with an outer part $$P_k^O$$ and an inner part $$P_k^I$$. Moreover, a homogeneous Dirichlet condition is applied at the boundary of the parts to make each of them like a new domain in which a linear system with a stiffness matrix Kk can be defined. (b) Solving the linear system in each part leads to n tensors Gsk. (c) Filtering Gsk, Hsk = C: Gsk and ρ, the effective properties are obtained in each inner part $$P_k^I$$. (d) The results from the n inner parts are gathered to obtain the effective properties in the whole inner domain ΩI. Figure 5. View largeDownload slide An alternative way of solving the non-periodic homogenization in parallel. (a) The domain Ω is decomposed in n parts Pk. As an example here n = 3. The parts overlap to equip them all with an outer part $$P_k^O$$ and an inner part $$P_k^I$$. Moreover, a homogeneous Dirichlet condition is applied at the boundary of the parts to make each of them like a new domain in which a linear system with a stiffness matrix Kk can be defined. (b) Solving the linear system in each part leads to n tensors Gsk. (c) Filtering Gsk, Hsk = C: Gsk and ρ, the effective properties are obtained in each inner part $$P_k^I$$. (d) The results from the n inner parts are gathered to obtain the effective properties in the whole inner domain ΩI. To speed up our code, multithreaded computations are performed whenever possible. Moreover, efficient algorithms based on k–d trees and stack data structures (e.g. Cormen et al.2009) are used to search for elements or points across the finite element mesh. Last but not least, the partitioning of the mesh is performed on ΩI instead of Ω to obtain well-balanced parts Pk (Fig. 6). The performance and capabilities of the code are illustrated through examples in the next part. Figure 6. View largeDownload slide Strategy for building well-balanced parts. (a) We partition ΩI instead of the whole domain Ω, then forming n subdomains $$\Omega ^I_k$$. As an example here n = 3. (b) Each subdomain is enlarged by a buffer layer $$\Omega ^B_k$$. At the border of the domain, the buffer is the outer domain ΩO (dotted grey). We end up with n parts $$P_k\!=\!\Omega ^I_k\cup \Omega ^B_k$$ which all have the same size. Moreover, $$P_k^I = \Omega _k^I$$ and $$P_k^O = \Omega ^B_k$$. Figure 6. View largeDownload slide Strategy for building well-balanced parts. (a) We partition ΩI instead of the whole domain Ω, then forming n subdomains $$\Omega ^I_k$$. As an example here n = 3. (b) Each subdomain is enlarged by a buffer layer $$\Omega ^B_k$$. At the border of the domain, the buffer is the outer domain ΩO (dotted grey). We end up with n parts $$P_k\!=\!\Omega ^I_k\cup \Omega ^B_k$$ which all have the same size. Moreover, $$P_k^I = \Omega _k^I$$ and $$P_k^O = \Omega ^B_k$$. 4 VALIDATION TESTS We here handle three different models to test the accuracy of the homogenization method in 3-D and challenge our code. First, the case of fine layers is investigated. For such a medium, analytical expressions of the effective elastic parameters exist so we can just compare the result of the homogenization with the result of these expressions to validate our code. The second model we study is made of small cubes. Because no reference solution for the effective properties is available in this case, we base our validation on the comparison of waveforms computed in the initial medium on the one hand, and in the homogenized medium on the other hand. Finally, a 3-D geological model is handled to emphasize the efficiency of the homogenization in a realistic case. 4.1 Homogenization of a finely layered medium The first model we consider for testing our 3-D homogenization code is a medium made of 60 isotropic layers randomly taken between 800 and 1200 m thick. Within each layer, the density is randomly chosen in the 2000–4000 kg m−3 range. Because the medium is isotropic, two parameters (e.g. the Lamé coefficients, or the S- and P-wave velocities) are sufficient to define the elastic tensor in each layer. We then randomly choose the S-wave velocity between 3000 and 5000 m s−1, and the P-wave velocity between 5000 and 8000 m s−1 (Fig. 7), with the constraint of having the Poisson’s ratio in the 0.1–0.45 range to make it geologically realistic. Figure 7. View largeDownload slide Comparison of effective properties computed with the 3-D non-periodic homogenization (dashed red) and the Backus theory (black) in the case of an original layered medium (blue). The two equivalent media both are anisotropic ($$\phi = \frac{V_{PV}^2}{V_{PH}^{2}}$$ and $$\xi = \frac{V_{SH}^{2}}{V_{SV}^{2}}$$ pointing out the rate of anisotropy) and agree very well with each other. Figure 7. View largeDownload slide Comparison of effective properties computed with the 3-D non-periodic homogenization (dashed red) and the Backus theory (black) in the case of an original layered medium (blue). The two equivalent media both are anisotropic ($$\phi = \frac{V_{PV}^2}{V_{PH}^{2}}$$ and $$\xi = \frac{V_{SH}^{2}}{V_{SV}^{2}}$$ pointing out the rate of anisotropy) and agree very well with each other. The homogenization of the model is performed using λ0 = 1600 m. In this case, the thickness of the outer domain (i.e. half of the filtering wavelet support) is 6400 m (Appendix C), so we take the extent of the layers equal to 15 km in order to get an inner domain in which a solution can be computed. Following a rule of thumb for the spatial discretization (Appendix E), each layer is meshed by a single layer of tetrahedral elements equipped with interpolation functions of degree 3. As expected from Backus (1962), the resulting effective medium is transversely isotropic. In Fig. 7, we represent it in terms of vertically- and horizontally polarized wave velocity:   \begin{eqnarray*} V_{PV} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{zzzz}}{\rho ^{\varepsilon _0\star }}}, \nonumber \\ V_{PH} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{xxxx}}{\rho ^{\varepsilon _0\star }}} = \sqrt{\frac{C^{\varepsilon _0\star }_{yyyy}}{\rho ^{\varepsilon _0\star }}}, \nonumber \\ V_{SV} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{xzxz}}{\rho ^{\varepsilon _0\star }}} = \sqrt{\frac{C^{\varepsilon _0\star }_{yzyz}}{\rho ^{\varepsilon _0\star }}}, \nonumber \\ V_{SH} &=& \sqrt{\frac{C^{\varepsilon _0\star }_{xyxy}}{\rho ^{\varepsilon _0\star }}}, \nonumber \end{eqnarray*} z being the vertical, layered direction, and (x, y) defining the horizontal plane. Moreover, we emphasize the amount of anisotropy by plotting $$\phi = \frac{V_{PV}^2}{V_{PH}^2}$$ and $$\xi = \frac{V_{SH}^2}{V_{SV}^2}$$. The effective properties computed with the non-periodic homogenization in the layered case are expected to be the same as those proposed by Backus (1962). A comparison between the two solutions (Fig. 7) shows a very good match, meaning that our 3-D implementation of the non-periodic homogenization works properly in the present case. As mentioned in Part 3, no solution can be computed at the border of the domain because the convolution involved in the filtering operation is not possible there. Nonetheless, extending the domain with a relevant material to make the convolution possible and to get satisfying effective properties at the border would be possible. Simply extending the boundary values of density and elastic coefficients would lead to the zeroth-order solution (Capdeville & Marigo 2007). To reach the first order, a continuous periodic extension (in which the boundary actually acts as a mirror) would have to be used (Leptev 2005; Capdeville & Marigo 2007). We do not implement such extensions in the present work. 4.2 Homogenization of random cubes To push our validation further, we challenge our homogenization code to a highly heterogeneous medium made of small elastic cubes with random isotropic properties. Each cube is 1 km3 large. 100 cubes are considered in each direction, which gives rise to a large cubic volume made of 1 000 000 random cubes. As shown in Fig. 8(a), this cubic volume is embedded in a 13 km thick homogeneous medium. Such a thickness corresponds to the support of the filtering wavelet which will be used in the homogenization process. In each small cube, the properties are randomly picked between 2000 and 4000 kg m−3 for the density, 2500 and 5000 m s−1 for the S-wave velocity, and 4000 and 8000 m s−1 for the P-wave velocity. As in the previous example, the Poisson’s ratio is constrained in the 0.1–0.45 range. Figure 8. View largeDownload slide (a) Cut in the random cubes. The black lines emphasize the edges of the whole domain. (b) Cut in the homogenized medium. (c) Snapshot of the L2 norm of the wavefield u generated through the random cubes by a force along the z-axis at point xS marked by the blue star. (d) Snapshot of the L2 norm of the wavefield u0 generated through the homogenized medium by the exact same force. Figure 8. View largeDownload slide (a) Cut in the random cubes. The black lines emphasize the edges of the whole domain. (b) Cut in the homogenized medium. (c) Snapshot of the L2 norm of the wavefield u generated through the random cubes by a force along the z-axis at point xS marked by the blue star. (d) Snapshot of the L2 norm of the wavefield u0 generated through the homogenized medium by the exact same force. The homogenization of the random cubes is performed using λ0 = 1600 m (i.e. λm = 8000 m and ε0 = 0.2). The spatial discretization is ensured by dividing every single cube into six tetrahedra equipped with degree 3 interpolants. Using similar elements in the homogeneous region, we end up with 12 002 256 tetrahedra and 160 747 899 degrees of freedom (three components at 53 582 633 interpolation points). To achieve the computation, the domain is split into 100 overlapping subdomains. Calculating the effective properties in a single subdomain then requires about 116 GB and 4 hr on an Intel Xeon X5680 processor (6 cores, 3.33 GHz, 12 MB Cache, 6.4 GT s−1). The obtained homogenized model is shown in Fig. 8(b). To assess the quality of the effective medium calculated with the non-periodic homogenization technique, we perform a seismic wave simulation in it and we compare the obtained seismograms with traces computed in the initial model (i.e. the random cubes). In other words, we solve the initial problem   $$\left\lbrace \begin{array}{l}\rho \ddot{\boldsymbol {u}} - \boldsymbol {\nabla }\cdot \boldsymbol {\sigma } = \boldsymbol {f} \\ \boldsymbol {\sigma } = {\bf C}: \left[ \frac{1}{2} (\boldsymbol {\nabla }\boldsymbol {u} + {^t}\boldsymbol {\nabla }\boldsymbol {u}) \right] \end{array} \right.$$ (34)and the homogenized problem   $$\left\lbrace \begin{array}{l}\rho ^{\varepsilon _0\star } \ddot{\boldsymbol {u}}_0 - \boldsymbol {\nabla }\cdot \boldsymbol {\sigma }_0 = \boldsymbol {f} \\ \boldsymbol {\sigma }_0 = {\bf C}^{\varepsilon _0\star } : \left[ \frac{1}{2} (\boldsymbol {\nabla }\boldsymbol {u}_0 + {^t}\boldsymbol {\nabla }\boldsymbol {u}_0) \right] \end{array} \right.$$ (35)to estimate the quality of $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$ through the comparison of u0 to u. In eq. (34) and (35), $$\,\ddot{ }\,$$ represents the second time-derivative and f is the external force. This latter is chosen as a simple Ricker function R(t) along the z-axis at a given point xS in the homogeneous region: f(x, t) = R(t) δ(x − xS) ez. The dominant frequencies of R(t) is chosen to be equal to 0.2 Hz. The two simulations are performed using a spectral element method with PML-type absorbing boundaries (Cupillard et al.2012). Snapshots of the two obtained wavefields are shown in Figs 8(c) and (d). They look very similar, suggesting that our homogenized model is an accurate equivalent medium for the seismic wave propagation. We carefully compare u0 to u by looking at seismograms calculated at two particular receivers denoted by A and B in Figs 8(c) and (d). Receiver A is on a P-wave nodal plane, 92 km away from the source. The ballistic S-wave, which contains most of the seismic energy, appears on the z-component (Fig. 9). We see that this wave is very well-retrieved by the homogenized model, the difference between u0 and u (i.e. the residual) reaching 8   per  cent at most. Apart from the S-wave, a scattered wavefield is observed on the three components. Our homogenized model also reconstructs this wavefield very well. To emphasize the relevance of the homogenized solution, Fig. 9 also shows seismograms computed in a medium obtained by just filtering the initial density and elastic tensor. The displacement and the stress propagating in this medium are denoted by $$\boldsymbol {u}_\mathcal {F}$$ and $$\boldsymbol {\sigma }_\mathcal {F}$$, respectively. By definition, these two fields verify   $$\left\lbrace \begin{array}{l}\mathcal {F}^{\varepsilon _0}_3\lbrace \rho \rbrace \ddot{\boldsymbol {u}}_\mathcal {F} - \boldsymbol {\nabla }\cdot \boldsymbol {\sigma }_\mathcal {F} = \boldsymbol {f} \\ \boldsymbol {\sigma }_\mathcal {F} = \mathcal {F}^{\varepsilon _0}_3\lbrace {\bf C}\rbrace : \left[ \frac{1}{2} (\boldsymbol {\nabla }\boldsymbol {u}_\mathcal {F} + {^t}\boldsymbol {\nabla }\boldsymbol {u}_\mathcal {F}) \right]. \end{array} \right.$$ (36)The waveforms of $$\boldsymbol {u}_\mathcal {F}$$ at receiver A do not fit the wavefield u in the initial model. Neither the phase nor the amplitude is properly reconstructed by the filtered medium, meaning that this latter does not hold the correct effective properties for the seismic wave propagation. Similar features are observed at receiver B (Fig. 10): $$\boldsymbol {u}_\mathcal {F}$$ is far from u whereas the homogenized solution accurately recovers the whole seismograms, including a ballistic P-wave which appears at this receiver location. Figure 9. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver A (Figs 8c and d). The maximum and dominant frequencies of the Ricker function R(t) emitted at the source are equal to 0.5 and 0.2 Hz, respectively. Figure 9. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver A (Figs 8c and d). The maximum and dominant frequencies of the Ricker function R(t) emitted at the source are equal to 0.5 and 0.2 Hz, respectively. Figure 10. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver B (Figs 8c and d). Contrary to receiver A, B records a ballistic P-wave. Figure 10. View largeDownload slide Comparison of the wavefields u (black), u0 (dashed red) and $$\boldsymbol {u}_\mathcal {F}$$ (green) computed at receiver B (Figs 8c and d). Contrary to receiver A, B records a ballistic P-wave. The overall difference between u0 and u can be evaluated quantitatively using the error   $$E_0 = \frac{1}{50} \sum _{r=1}^{50} \sqrt{ \frac{\int _0^T (\boldsymbol {u}_0-\boldsymbol {u})^2(\boldsymbol {x}_r) \, \text{d}t}{\int _0^T \boldsymbol {u}^2(\boldsymbol {x}_r) \, \text{d}t} } \, ,$$ (37)where xr is the randomly chosen location of receiver r. In eq. (37), T is the duration of the simulated propagation. It is equal to 55 s here. For ε0 = 0.2 (i.e. λ0 = 1600 m), we obtain E0 ≃ 0.006. Computing the same kind of error for the naive solution $$\boldsymbol {u}_\mathcal {F}$$, we get $$E_\mathcal {F} \simeq 0.055 \simeq 9E_0$$. This result again emphasizes the much higher accuracy of the homogenized solution u0. We also compute E0 and $$E_\mathcal {F}$$ for ε0 = 0.4 (i.e. λ0 = 3200 m), ε0 = 0.8 (i.e. λ0 = 6400 m) and ε0 = 1.6 (i.e. λ0 = 12 800 m). This allows us to study how fast u0 and $$\boldsymbol {u}_\mathcal {F}$$ converge to u. From the homogenization theory, we expect E0 = O(ε0). Fig. 11 shows that we actually get $$E_0 \simeq O(\varepsilon _0^{3/2})$$. This unexpectedly fast convergence of u0 toward u could be due to the weak contribution of the higher-order terms in the present case study, as suggested by Capdeville et al. (2010b) in 2-D. Plotting also $$E_\mathcal {F}$$ in Fig. 11, it appears that the convergence of $$\boldsymbol {u}_\mathcal {F}$$ is way poorer. Figure 11. View largeDownload slide E0 (red) and $$E_\mathcal {F}$$ (green) as a function of ε0. It is clear from this plot that u0 converges much faster than $$\boldsymbol {u}_\mathcal {F}$$ toward the target wavefield u. Figure 11. View largeDownload slide E0 (red) and $$E_\mathcal {F}$$ (green) as a function of ε0. It is clear from this plot that u0 converges much faster than $$\boldsymbol {u}_\mathcal {F}$$ toward the target wavefield u. The spectral element simulation of the wavefield u in the initial model requires 126 hexahedra in each direction, that is, 2 000 376 elements in total. This is because a spectral element mesh has to honour all the physical discontinuities of the model in study to make the computation accurate, so each 1 km × 1 km × 1 km random cube has to be captured by a hexahedron. As a consequence, the obtained mesh highly oversamples the wavefield: because λm = 8 km, elements as large as 8 km × 8 km × 8 km with a polynomial order larger than 4 would be sufficient if the medium was smooth enough (i.e. if it only contained scales larger than λm). Because of the discontinuities, we here get a 512 times denser mesh, which yields a 4096 times higher numerical cost (a factor 512 in space times a factor 8 in time because of eq. 1). Computing a 55 s long simulation of u then takes about 6 hr on ten Intel Xeon X5680 processors. In homogenized media, such as those computed here using ε0 = 0.2, 0.4, 0.8 and 1.6, coarser spectral element meshes can be used, then decreasing the computation cost. When using ε0 = 0.4 for instance, we get an effective medium of the random cubes that only contains scales larger than λ0 = 3.2 km, so 3.2 km3 large hexahedra then suit an accurate spectral element simulation of u0, which is 105 times less numerically demanding than the simulation of u. Nevertheless, all the wave simulations presented in this paper are performed within the same fine hexahedral mesh to avoid possible numerical biases and thus make all our comparisons totally fair. 4.3 Homogenization of a realistic geological model In this section, we use a subsurface model of the Ribaute area in France (Caumon et al.2009) as an illustration of the 3-D non-periodic homogenization technique applied to a realistic geological medium (Fig. 12a). Simulating waves in the model as it is (i.e. composed of multiple faulted and folded horizons) is extremely challenging because of the fine grid required to accommodate thin layers and complex geometries formed by the discontinuities. Such a fine grid indeed makes the computation cost of any wave simulation enormous. Moreover, in the context of the spectral element method, generating a hexahedral mesh which honours all the horizons and faults, if possible, demands tremendous efforts. For these reasons, we do not compute reference waveforms here; we just perform a spectral element simulation in the homogenized model to show how convenient working with effective properties is, and we compare the results with a simulation in a model obtained by brutally filtering the initial density and elastic tensor to exemplify that the choice of the upscaling technique matters. Figure 12. View largeDownload slide (a) Structural model of a highly faulted and folded region near Ribaute, Southern France. (b) Adaptive tetrahedral mesh of the model. This mesh is used to perform the homogenization. The background colour corresponds to the (randomly chosen) isotropic S-wave velocity within each layer of the model. (c) One of the S-wave velocities resulting from the homogenization. Another S-wave velocity is shown in Fig. 13. The difference between these two velocities is plotted in (d), which emphasizes the seismic anisotropy produced by the structure. Figure 12. View largeDownload slide (a) Structural model of a highly faulted and folded region near Ribaute, Southern France. (b) Adaptive tetrahedral mesh of the model. This mesh is used to perform the homogenization. The background colour corresponds to the (randomly chosen) isotropic S-wave velocity within each layer of the model. (c) One of the S-wave velocities resulting from the homogenization. Another S-wave velocity is shown in Fig. 13. The difference between these two velocities is plotted in (d), which emphasizes the seismic anisotropy produced by the structure. The initial model is made of five homogeneous isotropic layers (Fig. 12b). To compute its effective properties using our finite-element code, it is meshed with VorteXLib from RINGMesh (Pellerin et al.2017). We choose λ0 = 30 m and polynomial interpolants of degree 2, so the optimal volume for the tetrahedra is about 200 m3 (Appendix E). While smaller elements are needed around the discontinuities for capturing their complex geometry, larger elements are allowed where the discontinuities have no influence, that is, at a distance larger than the size of the filtering wavelet support. We therefore take advantage of the refinement-derefinement technique proposed by Botella (2016) to generate an adaptive mesh which minimizes the finite-element computation cost. We end up with a 9 077 300 element mesh and 35 703 579 degrees of freedom. The homogenized model is obtained in 9 min on 120 Intel Xeon E5-2683 v4 (16 cores, 2.10 GHz, 40 MB Cache, 9.6 GT s−1). The result is shown in Fig. 12(c). We have plotted the $$\sqrt{\frac{C^{\varepsilon _0\star }_{xyxy}}{\rho ^{\varepsilon _0\star }}}$$ component there to make the figure comparable to Fig. 12(b) where the isotropic S-wave velocity of the initial medium is represented. As expected, the homogeneous areas (i.e. the interior of the layers) are not changed by the homogenization process, whereas all the discontinuities (i.e. the faults and the interfaces between the layers) are smoothed. To illustrate the S-anisotropy produced by these discontinuities, $$\frac{ \sqrt{C^{\varepsilon _0\star }_{xyxy}} - \sqrt{C^{\varepsilon _0\star }_{xzxz}} }{ \sqrt{\rho ^{\varepsilon _0\star }} }$$ is plotted in Fig. 12(d). We observe that the anisotropy can reach 20   per  cent where high velocity contrasts (i.e. at the bottom and the top of the fastest layer) occur. Because the homogenized model is smooth, seismic waves within it can be simulated using a coarse mesh. As an example, a coarse hexahedral mesh supporting a spectral element simulation (Cupillard et al.2012) is shown in Fig. 13. The elements there are 60 m3 large and hold degree-8 polynomials to capture all the variations of the model. We see that the direct seismic wave front is highly deformed by the effective structure and that strong reflected and diffracted waves are generated, even though the medium contains no discontinuities. Fig. 13 also shows waveforms obtained at two stations. At station A, near the source, several S-waves reflected on discontinuities are observed. Zooming in two of them emphasizes important discrepancies between our homogenized solution and the waveform computed in a medium obtained by just filtering the density and the elastic tensor of the initial model. These discrepancies can be explained by the lack of anisotropy in this last medium. By construction, it is fully isotropic, so it cannot hold the important structure-induced anisotropy observed in our homogenized model (e.g. Fig. 12d). As a consequence, the waveforms computed in the two models at station A show significant differences in phase and amplitude. Similar observations are also made far from the source, at station B. In the light of the results presented in Section 4.2 and in previous 2-D studies (Capdeville et al.2010b; Guillot et al.2010), we can reasonably think that the waveforms computed in the homogenized model are the most accurate, but we cannot assess it because we do not provide any reference solution here. The main point of this section is to show an application of our homogenization code to a 3-D realistic geological model and to put the importance of the structure-induced anisotropy in evidence. Figure 13. View largeDownload slide Left: three snapshots of the L2 norm of the wavefield u0 generated through the homogenized Ribaute model by a force along the x-axis acting at the point marked by the yellow star. The wave simulation is performed using a spectral element method on a regular hexahedral mesh. Top right: x-component of the displacement u0 recorded near the source, at station A. The displacement $$\boldsymbol {u}_\mathcal {F}$$ computed in a medium obtained by just filtering the density and the elastic tensor of the initial model is also plotted. As expected, the direct S-wave (arriving around t = 0.05 s) is identical in the two simulations. On the contrary, the reflected phases show significant differences. In particular, zooming in the waves reflected on the two major velocity contrasts (i.e. the bottom and the top of the fastest layer) shows that $$\boldsymbol {u}_\mathcal {F}$$ arrives earlier than u0 with a different amplitude. Bottom right: x- and z-components of both u0 and $$\boldsymbol {u}_\mathcal {F}$$ recorded far from the source, at station B. Because the direct S-wave has travelled through heterogeneities, it is now different in the two simulations. Obviously, discrepancies in the later, multiply reflected and scattered energy, are also observed. Figure 13. View largeDownload slide Left: three snapshots of the L2 norm of the wavefield u0 generated through the homogenized Ribaute model by a force along the x-axis acting at the point marked by the yellow star. The wave simulation is performed using a spectral element method on a regular hexahedral mesh. Top right: x-component of the displacement u0 recorded near the source, at station A. The displacement $$\boldsymbol {u}_\mathcal {F}$$ computed in a medium obtained by just filtering the density and the elastic tensor of the initial model is also plotted. As expected, the direct S-wave (arriving around t = 0.05 s) is identical in the two simulations. On the contrary, the reflected phases show significant differences. In particular, zooming in the waves reflected on the two major velocity contrasts (i.e. the bottom and the top of the fastest layer) shows that $$\boldsymbol {u}_\mathcal {F}$$ arrives earlier than u0 with a different amplitude. Bottom right: x- and z-components of both u0 and $$\boldsymbol {u}_\mathcal {F}$$ recorded far from the source, at station B. Because the direct S-wave has travelled through heterogeneities, it is now different in the two simulations. Obviously, discrepancies in the later, multiply reflected and scattered energy, are also observed. 5 DISCUSSION AND CONCLUSIONS Dealing with small-scale heterogeneities in seismic wave simulation is a difficult task because it usually involves enormous computation costs. To handle them, Graphics Processing Unit (GPU) and/or High Performance Computing (HPC) implementations of well-established numerical techniques such as the Spectral Element method (SEM), the Discontinuous Galerkin method (DGM) and the Finite Difference method (FDM), have been proposed (Komatitsch et al.2010; Peter et al.2011; Weiss & Shragge 2013; Gokhberg & Fichtner 2016; Remacle et al.2016; Rietmann et al.2017). Furthermore, the DGM allows local time-stepping and p-adaptivity (e.g. Dumbser et al.2007; Etienne et al.2010; Minisini et al.2013; Diaz & Grote 2015) which mitigate stability constraint (1) and therefore reduce the overall computation cost. In the context of the SEM and the FDM, Pelties et al. (2011) proposed empirical laws on the mesh size to relax the need of honouring geological discontinuities. Complementary to all these numerical advances, the use of effective properties for the seismic wave propagation can drastically reduce the computation cost related to small-scale features while preserving a good accuracy. In the recent years, the non-periodic homogenization method emerged as a general technique to compute such effective properties. We here applied it in 3-D for the first time. In part 2, we recalled the theory of the homogenization, skipping some technical details to focus on the main ideas and concepts. Then we described an efficient, embarrassingly parallel implementation of the method. In part 4, we challenged this implementation on various media, showing the high accuracy of the non-periodic homogenization and the ability of our code to handle large and complex 3-D models, with no restriction on the size and shape of the heterogeneities. The code is available at upon request. Homogenizing complex geological media to ease the numerical simulation of full seismic wavefields has straight-forward applications in seismic risk assessment (e.g. Chaljub et al.2010), survey design (e.g. Wei et al.2012), structural model validation (e.g. Irakarama et al.2017), seismic source characterization (e.g. Silwal & Tape 2016) and seismic tomography for keeping subwavelength geological details in the inversion (Fichtner et al.2013b; Capdeville & Cance 2015). Because it tells what the waves ‘see’, the non-periodic homogenization also opens important perspectives in the interpretation of full waveform inversion results. In particular, the structure-induced (or, equivalently, ‘extrinsic’ or ‘apparent’ or ‘geometric’) anisotropy (Fichtner et al.2013a; Wang et al.2013) and the structure-induced attenuation (Ricard et al.2014) could be downscaled, as recently initiated by Bodin et al. (2015), leading to probability distributions of small-scale features and to better estimations of the intrinsic anisotropy and attenuation at various scales, from the Earth’s core up to the subsurface. In addition, Capdeville et al. (2013) and Afanasiev et al. (2016) recently showed that the homogenization can help in regularizing full waveform inversion problems. The non-periodic homogenization relies on a two-scale asymptotic expansion of the displacement and the stress involved in the elastic wave equation. In this paper, we have investigated the zeroth-order term and the associated upscaled properties. As shown by Capdeville et al. (2010a,b) and Guillot et al. (2010), adding the first-order term allows retrieving site effects (i.e. small-amplitude high-frequency non-propagating signals) at the receivers. These authors also show that a correction can be applied to the external force term to take into account the effects of the local small-scale structure at the source. A recent application of this correction can be found in Burgos et al. (2016). We have claimed that our method does not require any constraint on the size and shape of the heterogeneities to be smoothed. In other words, the homogenization is able to upscale any media. The only limitation that has been noted so far occurs when trying to model subwavelength focusing in Helmholtz resonators. In this extreme case, the wavefield no longer holds a minimum wavelength, so the non-periodic homogenization fails (Zhao et al.2016). Moreover, we have mentioned that our homogenization theory is not able to compute long-wave equivalent properties of heterogeneities laying near Neumann or Dirichlet surfaces yet. A proper solution only exists for smoothing a free-surface topography on top of a homogeneous material (Capdeville & Marigo 2013) or fine horizontal layers below a flat free-surface (Capdeville & Marigo 2008). Relevant extensions of the model beyond its boundaries can also leads to accurate effective properties (Leptev 2005; Capdeville & Marigo 2007). Pushing these approaches would surely help the non-periodic homogenization theory in handling heterogeneities near various boundary shapes and conditions. Besides these physical and theoretical limitations, slight numerical weaknesses can be pointed out in our code. Even if it employs powerful algorithms and an efficient parallel scheme, it could be optimized using Basic Linear Algebra Subprograms and the Fortran column-major order more extensively, particularly when computing the stiffness matrix and when filtering the strain and the stress associated with the starting corrector. The performance of the code could also benefit from the nested homogenization technique proposed by Capdeville et al. (2015) and from an iterative solver which would be more efficient than a direct solver when dealing with very large stiffness matrices. Most of all, developing an adaptive homogenization, taking into account the fact that the minimum wavelength usually varies within the medium, would be a major advance. Acknowledgements We deeply thank Q. Liu and C. Boehm for their sharp reviews and constructive comments that significantly improve the initial manuscript. We also thank A. Fichtner, L. Guillot, P. Cance, Y. Ricard, F. Chambat, J.-P. Montagner, J. Virieux, A. Kutsenko, Y. Masson, G. Caumon and the European COST action TIDES (ES1401) for emulative remarks and discussions. We are also grateful to A. Botella for advising some powerful search algorithms and for providing 3-D adaptive meshes which significantly reduced the computation cost of the homogenization. Many thanks to A. Borghi as well, who completely reset the Cauchy cluster at GeoRessources. Part of the calculations has also been led on the S-CAPAD cluster at IPGP, on the Centre de Calcul Intensif des Pays de la Loire (CCIPL) computers, and on the EXPLOR centre hosted by the University of Lorraine. Financial supports from the RING-Gocad consortium and the ANR projects ANR-10-BLAN-613 and ANR-16-CE31-0022-01 also enabled this work. REFERENCES Abdulle A., Grote M.J., 2011. Finite element heterogeneous multiscale method for the wave equation, Multiscale Model. Simul. , 9( 2), 766– 792. https://doi.org/10.1137/100800488 Google Scholar CrossRef Search ADS   Afanasiev M., Boehm C., May D., Fichtner A., 2016. Using effective medium theory to better constrain full waveform inversion, First Break , 34( 7), 93– 94. Allaire G., 1992. Homogenization and two-scale convergence, SIAM J. Math. Anal. , 23, 1482– 1518. https://doi.org/10.1137/0523084 Google Scholar CrossRef Search ADS   Allaire G., Capdeboscq Y., 2000. Homogenization of a spectral problem in neutronic multigroup diffusion, Comput. Methods Appl. Mech. Eng. , 187( 1–2), 91– 117. https://doi.org/10.1016/S0045-7825(99)00112-7 Google Scholar CrossRef Search ADS   Allaire G., Habibi Z., 2013. Homogenization of a conductive, convective, and radiative heat transfer problem in a heterogeneous domain, SIAM J. Math. Anal. , 45( 3), 1136– 1178. https://doi.org/10.1137/110849821 Google Scholar CrossRef Search ADS   Bacigalupo A., Gambarotta L., 2014. Second-gradient homogenized model for wave propagation in heterogeneous periodic media, Int. J. Solids Struct. , 51( 5), 1052– 1065. https://doi.org/10.1016/j.ijsolstr.2013.12.001 Google Scholar CrossRef Search ADS   Backus G., 1962. Long-wave elastic anisotropy produced by horizontal layering, J. geophys. Res. , 67( 11), 4427– 4440. https://doi.org/10.1029/JZ067i011p04427 Google Scholar CrossRef Search ADS   Bensoussan A., Lions J.-L., Papanicolaou G., 1978. Asymptotic Analysis for Periodic Structures , North-Holland. Bodin T., Capdeville Y., Romanowicz B., Montagner J.-P., 2015. Interpreting radial anisotropy in global and regional tomographic models, in The Earth’s Heterogeneous Mantle , pp. 105– 144, eds Khan A., Deschamps F., Springer. Google Scholar CrossRef Search ADS   Botella A., 2016. Génération de maillages non structurés volumiques de modèles géologiques pour la simulation de phénomènes physiques, PhD thesis , Université de Lorraine. Boutin C., Auriault J.L., 1993. Rayleigh scattering in elastic composite materials, Int. J. Eng. Sci. , 31( 12), 1669– 1689. https://doi.org/10.1016/0020-7225(93)90082-6 Google Scholar CrossRef Search ADS   Brossier R., Operto S., Virieux J., 2015. Velocity model building from seismic reflection data by full waveform inversion, Geophys. Prospect. , 63( 2), 354– 367. https://doi.org/10.1111/1365-2478.12190 Google Scholar CrossRef Search ADS   Burgos G., Capdeville Y., Guillot L., 2016. Homogenized moment tensor and the effect of near-field heterogeneities on nonisotropic radiation in nuclear explosion, J. geophys. Res. , 121( 6), 4366– 4389. https://doi.org/10.1002/2015JB012744 Google Scholar CrossRef Search ADS   Cance P., Capdeville Y., 2015. Validity of the acoustic approximation for elastic waves in heterogeneous media, Geophysics , 80( 4), T161– T173. https://doi.org/10.1190/geo2014-0397.1 Google Scholar CrossRef Search ADS   Capdeville Y., Cance P., 2015. Residual homogenization for elastic wave propagation in complex media, Geophys. J. Int. , 200, 984– 997. https://doi.org/10.1093/gji/ggu452 Google Scholar CrossRef Search ADS   Capdeville Y., Marigo J., 2007. Second order homogenization of the elastic wave equation for non-periodic layered media, Geophys. J. Int. , 170, 823– 838. https://doi.org/10.1111/j.1365-246X.2007.03462.x Google Scholar CrossRef Search ADS   Capdeville Y., Marigo J., 2008. Shallow layer correction for spectral element like methods, Geophys. J. Int. , 172, 1135– 1150. https://doi.org/10.1111/j.1365-246X.2007.03703.x Google Scholar CrossRef Search ADS   Capdeville Y., Marigo J.-J., 2013. A non-periodic two scale asymptotic method to take account of rough topographies for 2-D elastic wave propagation, Geophys. J. Int. , 192( 1), 163– 189. https://doi.org/10.1093/gji/ggs001 Google Scholar CrossRef Search ADS   Capdeville Y., Guillot L., Marigo J., 2010a. 1-D non-periodic homogenization for the seismic wave equation, Geophys. J. Int. , 181, 907– 910. Capdeville Y., Guillot L., Marigo J., 2010b. 2-D non-periodic homogenization to upscale elastic media for P-SV waves, Geophys. J. Int. , 182, 903– 922. https://doi.org/10.1111/j.1365-246X.2010.04636.x Google Scholar CrossRef Search ADS   Capdeville Y., Stutzmann E., Wang N., Montagner J.-P., 2013. Residual homogenization for seismic forward and inverse problems in layered media, Geophys. J. Int. , 194( 1), 470– 487. https://doi.org/10.1093/gji/ggt102 Google Scholar CrossRef Search ADS   Capdeville Y., Zhao M., Cupillard P., 2015. Fast Fourier homogenization for elastic wave propagation in complex media, Wave Motion , 54, 170– 186. https://doi.org/10.1016/j.wavemoti.2014.12.006 Google Scholar CrossRef Search ADS   Carcione J., Picotti S., Santos J., 2012. Numerical experiments of fracture-induced velocity and attenuation anisotropy, Geophys. J. Int. , 191( 3), 1179– 1191. Casarotti E., Stupazzini M., Lee S., Komatitsch D., Piersanti A., Tromp J., 2008. CUBIT and seismic wave propagation based upon the spectral-element method: an advanced unstructured mesher for complex 3D geological media, in Proceedings of the 16th International Meshing Roundtable , pp. 579– 597, eds Brewer M.L., Marcum D., Springer, Berlin. Google Scholar CrossRef Search ADS   Caumon G., Collon P., Le Carlier de Veslud C., Sausse J., Viseur S., 2009. Surface-based 3D modeling of geological structures, Math. Geosci. , 41( 8), 927– 945. https://doi.org/10.1007/s11004-009-9244-2 Google Scholar CrossRef Search ADS   Chaljub E., Moczo P., Tsuno S., Bard P.-Y., Kristek J., Käser M., Stupazzini M., Kristekova M., 2010. Quantitative comparison of four numerical predictions of 3D ground motion in the Grenoble valley, France, Bull. seism. Soc. Am. , 100( 4), 1427– 1455. https://doi.org/10.1785/0120090052 Google Scholar CrossRef Search ADS   Cormen T.H., Leiserson C.E., Rivest R.L., Stein C., 2009. Introduction to Algorithms , 3rd edn, The MIT Press. Courant R., Friedrichs K., Lewy H., 1928. Uber die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann. , 100, 32– 74. https://doi.org/10.1007/BF01448839 Google Scholar CrossRef Search ADS   Cupillard P., Delavaud E., Burgos G., Festa G., Vilotte J.-P., Capdeville Y., Montagner J.-P., 2012. RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale, Geophys. J. Int. , 188, 1203– 1220. https://doi.org/10.1111/j.1365-246X.2011.05311.x Google Scholar CrossRef Search ADS   Dal Maso G., 1993. An Introduction to Γ-convergence , Birkhäuser. Google Scholar CrossRef Search ADS   Diaz J., Grote M.J., 2015. Multi-level explicit local time-stepping methods for second-order wave equations, Comput. Methods Appl. Mech. Eng. , 291, 240– 265. https://doi.org/10.1016/j.cma.2015.03.027 Google Scholar CrossRef Search ADS   Dumbser M., Käser M., Toro E.F., 2007. An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes—V. Local time stepping and p-adaptivity, Geophys. J. Int. , 171( 2), 695– 717. https://doi.org/10.1111/j.1365-246X.2007.03427.x Google Scholar CrossRef Search ADS   Dumontet H., 1990. Homogénéisation et effets de bords dans les matériaux composites, PhD thesis , Univertité Paris 6. Engquist B., Holst H., Runborg O., 2009. Multi-scale methods for wave propagation in heterogeneous media , preprint arXiv:0911.2638. Etienne V., Chaljub E., Virieux J., Glinsky N., 2010. An hp-adaptive discontinuous Galerkin Þnite-element method for 3-D elastic wave modelling, Geophys. J. Int. , 183, 941– 962. https://doi.org/10.1111/j.1365-246X.2010.04764.x Google Scholar CrossRef Search ADS   Felippa C., 2004. A compendium of FEM integration formulas for symbolic work, Eng. Comput. , 21( 8), 867– 890. https://doi.org/10.1108/02644400410554362 Google Scholar CrossRef Search ADS   Fichtner A., Igel H., 2008. Efficient numerical surface wave propagation through the optimization of discrete crustal models - a technique based on non-linear dispersion curve matching (DCM), Geophys. J. Int. , 173( 2), 519– 533. https://doi.org/10.1111/j.1365-246X.2008.03746.x Google Scholar CrossRef Search ADS   Fichtner A., Kennett B.L.N., Igel H., Bunge H.-P., 2008. Theoretical background for continental and global scale full-waveform inversion in the time-frequency domain, Geophys. J. Int. , 175( 8), 665– 685. https://doi.org/10.1111/j.1365-246X.2008.03923.x Google Scholar CrossRef Search ADS   Fichtner A., Kennett B.L., Trampert J., 2013a. Separating intrinsic and apparent anisotropy, Phys. Earth planet. Inter. , 219, 11– 20. https://doi.org/10.1016/j.pepi.2013.03.006 Google Scholar CrossRef Search ADS   Fichtner A., Trampert J., Cupillard P., Saygin E., Taymaz T., Capdeville Y., Villasenor A., 2013b. Multi-scale full waveform inversion, Geophys. J. Int. , 194, 534– 556. https://doi.org/10.1093/gji/ggt118 Google Scholar CrossRef Search ADS   Fish J., Chen W., 2001. Higher-order homogenization of initial/boundary-value problem, J. Eng. Mech. , 127( 12), 1223– 1230. https://doi.org/10.1061/(ASCE)0733-9399(2001)127:12(1223) Google Scholar CrossRef Search ADS   Fish J., Chen W., 2004. Space-time multiscale model for wave propagation in heterogeneous media, Comput. Methods Appl. Mech. Eng. , 193, 4837– 4856. https://doi.org/10.1016/j.cma.2004.05.006 Google Scholar CrossRef Search ADS   Gao K., Chung E., Gibson R., Fu S., Efendiev Y., 2015. A numerical homogenization method for heterogeneous, anisotropic elastic media based on multiscale theory, Geophysics , 80( 4), D385– D401. https://doi.org/10.1190/geo2014-0363.1 Google Scholar CrossRef Search ADS   Gokhberg A., Fichtner A., 2016. Full-waveform inversion on heterogeneous HPC systems, Comput. Geosci. , 89, 260– 268. https://doi.org/10.1016/j.cageo.2015.12.013 Google Scholar CrossRef Search ADS   Grechka V., 2003. Effective media: A forward modeling view, Geophysics , 68( 6), 2055– 2062. https://doi.org/10.1190/1.1635059 Google Scholar CrossRef Search ADS   Grechka V., 2007. Multiple cracks in vti rocks: Effective properties and fracture characterization, Geophysics , 72( 5), D81– D91. https://doi.org/10.1190/1.2751500 Google Scholar CrossRef Search ADS   Grechka V., Kachanov M., 2006. Effective elasticity of fractured rocks: A snapshot of the work in progress, Geophysics , 71( 6), W45– W58. https://doi.org/10.1190/1.2360212 Google Scholar CrossRef Search ADS   Guillot L., Capdeville Y., Marigo J., 2010. 2-D non-periodic homogenization of the elastic wave equation: SH case, Geophys. J. Int. , 182, 1438– 1454. https://doi.org/10.1111/j.1365-246X.2010.04688.x Google Scholar CrossRef Search ADS   Hashin Z., Shtrikman S., 1963. A variational approach to the elastic behavior of multiphase materials, J. Mech. Phys. Solids , 11, 127– 140. https://doi.org/10.1016/0022-5096(63)90060-7 Google Scholar CrossRef Search ADS   Hill R., 1965. A self–consistent mechanics of composit materials, J. Mech. Phys. Solids , 13, 213– 222. https://doi.org/10.1016/0022-5096(65)90010-4 Google Scholar CrossRef Search ADS   Hornung U., 1997. Homogenization and Porous Media , Vol. 6 of Interdisciplinary Applied Mathematics Series, Springer Verlag. Google Scholar CrossRef Search ADS   Hughes T.J., 2012. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis , Courier Corporation. Irakarama M., Cupillard P., Caumon G., Sava P., 2017. Appraising structural interpretations using seismic data misfit functionals, in 79th EAGE Conference and Exhibition , EAGE. Jordan T.H., 2015. An effective medium theory for three-dimensional elastic heterogeneities, Geophys. J. Int. , 203( 2), 1343– 1354. https://doi.org/10.1093/gji/ggv355 Google Scholar CrossRef Search ADS   Komatitsch D., Erlebacher G., Göddeke D., Michéa D., 2010. High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster, J. Comput. Phys. , 229( 20), 7692– 7714. https://doi.org/10.1016/j.jcp.2010.06.024 Google Scholar CrossRef Search ADS   Kutsenko A., Shuvalov A., Norris A., 2013. On the quasistatic effective elastic moduli for elastic waves in three-dimensional phononic crystals, J. Mech. Phys. Solids , 61( 11), 2240– 2259. https://doi.org/10.1016/j.jmps.2013.06.003 Google Scholar CrossRef Search ADS   Lekić V., Panning M., Romanowicz B., 2010. A simple method for improving crustal corrections in waveform tomography, Geophys. J. Int. , 182, 265– 278. Leptev V., 2005. Two-scale extensions for non-periodic coefficients, arXiv:math.AP/0512123. Lin C., Saleh R., Milkereit B., Liu Q., 2017. Effective media for transversely isotropic models based on homogenization theory: With applications to borehole sonic logs, Pure appl. Geophys. , 174( 7), 2631– 2647. https://doi.org/10.1007/s00024-017-1565-3 Google Scholar CrossRef Search ADS   Mauge C., Kachanov M., 1994. Effective elastic properties of an anisotropic material with arbitrarily oriented interacting cracks, J. Mech. Phys. Solids , 42, 561– 584. https://doi.org/10.1016/0022-5096(94)90052-3 Google Scholar CrossRef Search ADS   Métivier L., Bretaudeau F., Brossier R., Operto S., Virieux J., 2013. Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures, Geophys. Prospect. , 62( 6), 1353– 1375. Google Scholar CrossRef Search ADS   Minisini S., Zhebel E., Kononov A., Mulder W.A., 2013. Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media, Geophysics , 78( 3), T67– T77. https://doi.org/10.1190/geo2012-0252.1 Google Scholar CrossRef Search ADS   Moulinec H., Suquet P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Eng. , 157( 1), 69– 94. https://doi.org/10.1016/S0045-7825(97)00218-1 Google Scholar CrossRef Search ADS   Nguetseng G., 1989. A general convergence result for a functional related to the theory of homogenization, SIAM J. Math. Anal. , 20( 3), 608– 623. https://doi.org/10.1137/0520043 Google Scholar CrossRef Search ADS   Papanicolaou G.C., Varadhan S. R.S., 1981. Boundary value problems with rapidly oscillating random coefficients, in Random Fields, Vol. I, II (Esztergom, 1979) , Vol. 27 of Colloq. Math. Soc. János Bolyai, pp. 835– 873, North-Holland. Parnell W., Abrahams I., 2008. Homogenization for wave propagation in periodic fibre-reinforced media with complex microstructure. I - Theory, J. Mech. Phys. Solids , 56( 7), 2521– 2540. https://doi.org/10.1016/j.jmps.2008.02.003 Google Scholar CrossRef Search ADS   Pellerin J., Botella A., Bonneau F., Mazuyer A., Chauvin B., Levy B., Caumon G., 2017. RINGMesh: a programming library for developing mesh-based geomodeling applications, Comput. Geosci. , 104, 93– 100. https://doi.org/10.1016/j.cageo.2017.03.005 Google Scholar CrossRef Search ADS   Pelties C., Käser M., Hermann V., Castro C.E., 2011. Regular versus irregular meshing for complicated models and their effect on synthetic seismograms, Geophys. J. Int. , 183, 1031– 1051. https://doi.org/10.1111/j.1365-246X.2010.04777.x Google Scholar CrossRef Search ADS   Peter D. et al.  , 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, Geophys. J. Int. , 186, 721– 739. https://doi.org/10.1111/j.1365-246X.2011.05044.x Google Scholar CrossRef Search ADS   Plessix R.-E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophys. J. Int. , 167( 2), 495– 503. https://doi.org/10.1111/j.1365-246X.2006.02978.x Google Scholar CrossRef Search ADS   Pratt R., Shin C., Hicks G., 1998. Gauss–Newton and full Newton methods in frequency domain seismic waveform inversion, Geophys. J. Int. , 133, 341– 362. https://doi.org/10.1046/j.1365-246X.1998.00498.x Google Scholar CrossRef Search ADS   Remacle J.-F., Gandham R., Warburton T., 2016. GPU accelerated spectral finite elements on all-hex meshes, J. Comput. Phys. , 324, 246– 257. https://doi.org/10.1016/j.jcp.2016.08.005 Google Scholar CrossRef Search ADS   Reuss A., 1929. Berechnung der Fliessgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle, J. Appl. Math. Mech. , 9( 1), 49– 58. Ricard Y., Durand S., Montagner J.-P., Chambat F., 2014. Is there seismic attenuation in the mantle?, Earth planet. Sci. Lett. , 388, 257– 264. https://doi.org/10.1016/j.epsl.2013.12.008 Google Scholar CrossRef Search ADS   Rietmann M., Grote M., Peter D., Schenk O., 2017. Newmark local time stepping on high-performance computing architectures, J. Comput. Phys. , 334, 308– 326. https://doi.org/10.1016/j.jcp.2016.11.012 Google Scholar CrossRef Search ADS   Sanchez-Palencia E., 1980. Non-homogenous Media and Vibration Theory, Vol. 127 of Lecture notes in physics , Springer-Verlag. Santugini-Repiquet K., 2007. Homogenization of the demagnetization field operator in periodically perforated domains, J. Math. Anal. Appl. , 334( 1), 502– 516. https://doi.org/10.1016/j.jmaa.2007.01.001 Google Scholar CrossRef Search ADS   Sayers C., 1998. Long-wave seismic anisotropy of heterogeneous reservoirs, Geophys. J. Int. , 132, 667– 673. https://doi.org/10.1046/j.1365-246X.1998.00456.x Google Scholar CrossRef Search ADS   Sayers C., Kachanov M., 1991. A simple technique for finding effective elastic constants of cracked solids for arbitrary crack orientation statistics, Int. J. Solids Struct. , 27( 6), 671– 680. https://doi.org/10.1016/0020-7683(91)90027-D Google Scholar CrossRef Search ADS   Sayers C., Kachanov M., 1995. Microcrack-induced elastic wave anisotropy of brittle rocks, J. geophys. Res. , 100( B3), 4149– 4156. https://doi.org/10.1029/94JB03134 Google Scholar CrossRef Search ADS   Schenk O., Gärtner K., 2006. On fast factorization pivoting methods for symmetric indefinite systems, Electron. Trans. Numer. Anal. , 23, 158– 179. Schoenberg M., Helbig K., 1997. Orthorhombic media: modeling elastic wave behavior in a vertically fractured earth, Geophysics , 62( 6), 1954– 1974. https://doi.org/10.1190/1.1444297 Google Scholar CrossRef Search ADS   Schoenberg M., Muir F., 1989. A calculus for finely layered anisotropic media, Geophysics , 54, 581– 589. https://doi.org/10.1190/1.1442685 Google Scholar CrossRef Search ADS   Schoenberg M., Sayers C.M., 1995. Seismic anisotropy of fractured rock, Geophysics , 60( 1), 204– 211. https://doi.org/10.1190/1.1443748 Google Scholar CrossRef Search ADS   Silwal V., Tape C., 2016. Seismic moment tensors and estimated uncertainties in southern Alaska, J. geophys. Res. , 121( 4), 2772– 2797. https://doi.org/10.1002/2015JB012588 Google Scholar CrossRef Search ADS   Tromp J., Tape C., Liu Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, Geophys. J. Int. , 160, 195– 216. https://doi.org/10.1111/j.1365-246X.2004.02453.x Google Scholar CrossRef Search ADS   Virieux J., Calandra H., Plessix R.-E., 2011. A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging, Geophys. Prospect. , 59( 5), 794– 813. https://doi.org/10.1111/j.1365-2478.2011.00967.x Google Scholar CrossRef Search ADS   Voigt W., 1889. Uber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper, Ann. Phys. , 274( 12), 573– 587. https://doi.org/10.1002/andp.18892741206 Google Scholar CrossRef Search ADS   Wang N., Montagner J.-P., Fichtner A., Capdeville Y., 2013. Intrinsic versus extrinsic seismic anisotropy: The radial anisotropy in reference earth models, Geophys. Res. Lett. , 40( 16), 4284– 4288. https://doi.org/10.1002/grl.50873 Google Scholar CrossRef Search ADS   Warner M., Guasch L., 2016. Adaptive waveform inversion: Theory, Geophysics , 81( 6), R429– R445. https://doi.org/10.1190/geo2015-0387.1 Google Scholar CrossRef Search ADS   Wei W., Fu L., Blacquière G., 2012. Fast multifrequency focal beam analysis for 3D seismic acquisition geometry, Geophysics , 77( 2), P11– P21. https://doi.org/10.1190/geo2010-0327.1 Google Scholar CrossRef Search ADS   Weiss R.M., Shragge J., 2013. Solving 3D anisotropic elastic wave equations on parallel GPU devices, Geophysics , 78( 2), F7– F15. https://doi.org/10.1190/geo2012-0063.1 Google Scholar CrossRef Search ADS   Worth D., Greenough C., Chin S., 2012. Pragmatic software engineering for computational science, in Handbook of Research on Computational Science and Engineering , Vol. 1, pp. 119– 149, eds Leng J., Sharrock W., Information Science Reference. Zhao M., Capdeville Y., Zhang H., 2016. Direct numerical modeling of time-reversal acoustic subwavelength focusing, Wave Motion , 67, 102– 115. https://doi.org/10.1016/j.wavemoti.2016.07.010 Google Scholar CrossRef Search ADS   Zijl W., Hendriks M.A., ’t Hart C.M.P., 2002. Numerical homogenization of the rigidity tensor in Hooke’s law using the node-based finite element method, Math. Geol. , 34( 3), 291– 322. https://doi.org/10.1023/A:1014894923280 Google Scholar CrossRef Search ADS   APPENDIX A: THE LOW-PASS FILTER IN THE 1-D CASE To separate the macroscopic and the microscopic scales, the non-periodic homogenization method requires the use of a low-pass filter $$\mathcal {F}^{\varepsilon _0}$$. Applied to any function $$f \! : \! \mathbb {R}\! \rightarrow \! \mathbb {R}$$, this filter removes all the scales smaller than λ0:   $$\mathcal {F}^{\varepsilon _0}\lbrace f\rbrace (x) = \int _\mathbb {R} f(x^{\prime }) w^{\varepsilon _0}(x-x^{\prime }) \, \text{d}x^{\prime }$$ (A1)where $$w^{\varepsilon _0}$$ is the wavelet whose spectrum is   $$\hat{w}^{\varepsilon _0}(k) = \left\lbrace \begin{array}{ll}1 \, {for }k \leqslant \frac{2\pi }{\lambda _0} \\ 0 \, {for} {k} > \frac{2 {\pi}}{\lambda _0}. \end{array} \right.$$ (A2)In this last expression, $$\, \hat{}\,$$ represents the Fourier transform. Note that the filter and the associated wavelet are indexed by ε0 because λ0 depends on the choice of ε0 (λ0 = ε0λm). If defined by spectrum (A2), the wavelet $$w^{\varepsilon _0}$$ has an infinite support, which is unmanageable in practice. For that reason, we introduce a cosine-taper to soften the cut-off:   $$\hat{w}^{\varepsilon _0}(k) = \left\lbrace \begin{array}{lll}1 \, {for }k \leqslant k_0 \\ \frac{1}{2} \left[ 1 + \cos \left( \pi \, \frac{k\,-\,k_0}{\Delta k} \right) \right] \, {for }k \in \left] k_0; \, k_0\!-\!\Delta k \right] \\ 0 \, {for }k > k_0 \end{array} \right.$$ (A3)where $$k_0=\frac{2\pi }{\lambda _0}$$. With such a definition, $$w^{\varepsilon _0}$$ can be truncated with a negligible error to get a compact support (Fig. A1). Figure A1. View largeDownload slide A 1-D spectrum $$\hat{w}^{\varepsilon _0}$$ (left) and its corresponding wavelet $$w^{\varepsilon _0}$$ (right). In this example, the length of the cosine-taper Δk is equal to 0.5 k0. Figure A1. View largeDownload slide A 1-D spectrum $$\hat{w}^{\varepsilon _0}$$ (left) and its corresponding wavelet $$w^{\varepsilon _0}$$ (right). In this example, the length of the cosine-taper Δk is equal to 0.5 k0. Figure A2. View largeDownload slide 3-D spectrum $$\hat{w}_3^{\varepsilon _0}$$ for Δk = 0.5 k0 (left) and its corresponding 3-D wavelet $$w_3^{\varepsilon _0}$$ (right). Because both $$\hat{w}_3^{\varepsilon _0}$$ and $$w_3^{\varepsilon _0}$$ are spherical, we represent them as a function of ‖k‖ and ‖x‖, respectively. Comparing the obtained plots with Fig. A1, we note that the 3-D wavelet is different from just having the wavelet of the 1-D case in all the directions. This is simply because the (inverse) spherical Fourier transform is not the (inverse) 1-D Fourier transform. Figure A2. View largeDownload slide 3-D spectrum $$\hat{w}_3^{\varepsilon _0}$$ for Δk = 0.5 k0 (left) and its corresponding 3-D wavelet $$w_3^{\varepsilon _0}$$ (right). Because both $$\hat{w}_3^{\varepsilon _0}$$ and $$w_3^{\varepsilon _0}$$ are spherical, we represent them as a function of ‖k‖ and ‖x‖, respectively. Comparing the obtained plots with Fig. A1, we note that the 3-D wavelet is different from just having the wavelet of the 1-D case in all the directions. This is simply because the (inverse) spherical Fourier transform is not the (inverse) 1-D Fourier transform. APPENDIX B: $${{\rho} ^{\varepsilon _0}}$$ MUST BELONG TO Ψ Inserting eqs (2)–(4) and properties $$E^{\varepsilon _0}(x,y)$$ and $$\rho ^{\varepsilon _0}(x,y)$$ instead of E(x) and ρ(x) into the initial elastodynamic problem yields a cascade of equations. For i = 0, it turns that   $$\rho ^{\varepsilon _0}(x,y)\ddot{u}_0(x) - \nabla _x\sigma _0(x) - \nabla _y\sigma _1(x,y) = f(x),$$ (B1)where f is the external force and $$\ddot{u}_0$$ is the second time-derivative of the zeroth-order displacement. Because σ1 has to lie in Ψ, $$\rho ^{\varepsilon _0}$$ must belong to Ψ as well. APPENDIX C: THE LOW-PASS FILTER IN THE 3-D CASE $$\mathcal {F}^{\varepsilon _0}_3$$ is the extension of the 1-D low-pass filter $$\mathcal {F}^{\varepsilon _0}$$ (Appendix A) to 3-D. Applied to any function $$f \! : \! \mathbb {R}^3 \! \rightarrow \! \mathbb {R}$$, it removes all the scales smaller than λ0:   $$\mathcal {F}^{\varepsilon _0}_3\lbrace f\rbrace (\boldsymbol {x}) = \int _{\mathbb {R}^3} f(\boldsymbol {x}^{\prime }) w_3^{\varepsilon _0}(\boldsymbol {x}-\boldsymbol {x}^{\prime }) \, \text{d}\boldsymbol {x}^{\prime },$$ (C1)where $$w_3^{\varepsilon _0}$$ is the wavelet whose 3-D spectrum is   \begin{eqnarray} \hat{w}_3^{\varepsilon _0}(\boldsymbol {k}) = \left\lbrace \begin{array}{lll} 1 \, {for }\Vert \boldsymbol {k}\Vert \leqslant k_0 \\ \frac{1}{2} \left[ 1 + \cos \left( \pi \, \frac{\Vert \boldsymbol {k}\Vert \,-\,k_0}{\Delta k} \right) \right] \, {for }\Vert \boldsymbol {k}\Vert \in \left] k_0; \, k_0\!-\!\Delta k \right] \\ 0 \, {for }\Vert \boldsymbol {k}\Vert > k_0. \end{array} \right. \end{eqnarray} (C2) k0 still is the cut-off wave-number $$\frac{2\pi }{\lambda _0}$$. For all the applications within the present paper (Section 4), the length of the cosine-taper Δk has been chosen to be equal to 0.5 k0. When using this value, a truncation of the wavelet at ‖x‖ = 4 λ0 is satisfactory (Fig. A2). APPENDIX D: $${\bf {H}}^{\varepsilon _0}$$ AND $${\chi}^{\varepsilon _0}$$ ARE PERIODIC IN y AND BELONG TO Ψ3 As a solution of a periodic boundary value problem, χs is periodic in y, so ∇yχs is periodic in y. Gs therefore is periodic in y, as well as Hs, consequently. Finally, $${\bf G}^{\varepsilon _0}$$ and $${\bf H}^{\varepsilon _0}$$ are both periodic in y. By construction (eqs 28 and 29), they also belong to Ψ3. Let us denote the cross product by ∧. Because ∇y∧Gs = 0 and $$\boldsymbol {\nabla }_{\!y} \wedge \mathcal {F}^{\varepsilon _0}_3\lbrace \boldsymbol {f}\rbrace = \mathcal {F}^{\varepsilon _0}_3\lbrace \boldsymbol {\nabla }_{\!y} \wedge \boldsymbol {f}\rbrace$$ for any tensorial function f, we have $$\boldsymbol {\nabla }_{\!y} \wedge {\bf G}^{\varepsilon _0} = \boldsymbol {0}$$, so $${\bf G}^{\varepsilon _0}$$ can be written as a constant plus a gradient of a displacement $$\boldsymbol {\chi }^{\varepsilon _0}$$. As a solution of eq. (23), $$\boldsymbol {\nabla }_{\!y} \boldsymbol {\chi }^{\varepsilon _0}$$ has no rotational component, so we end up with   $${\bf G}^{\varepsilon _0} = {\bf I}+\frac{1}{2} (\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}+{^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0}).$$ (D1)By construction (eq. 28), $$\left\langle {{\bf G}^{\varepsilon _0}}\right\rangle _3 = {\bf I}$$ (that would not be the case if $${\bf G}^{\varepsilon _0}$$ was built as $$E^{\varepsilon _0}$$ and $$\rho ^{\varepsilon _0}$$ are—eqs 17 and 18–), so   $$\left\langle {\, \frac{1}{2} \left( \boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} + {^t}\boldsymbol {\nabla }_y\boldsymbol {\chi }^{\varepsilon _0} \right)\,}\right\rangle _3 = 0.$$ (D2)For any tensorial function f ∈ Ψ3 periodic in y,   \begin{eqnarray} \left\langle {\boldsymbol {f}}\right\rangle _3 &=& 0 \ \mbox{and} \ \boldsymbol {\nabla }_{\!y} \boldsymbol {g}\nonumber\\ &=& \boldsymbol {f}\ \Longrightarrow \ \boldsymbol {g} \ \mbox{is periodic in} \ \boldsymbol {y}\ \mbox{and belongs to} \ \Psi _3 , \end{eqnarray} (D3)so $$\boldsymbol {\chi }^{\varepsilon _0}$$ is periodic in y and belongs to Ψ3. This demonstration can also be found in the seminal paper by Capdeville et al. (2010b), section 3.5. APPENDIX E: HEURISTIC FOR SAMPLING A MODEL TO BE HOMOGENIZED When honouring discontinuities, finite element meshes enable a correct account for the geometry of the medium in which the phenomenon to be simulated occurs. Nevertheless, the element size and the interpolation degree which guarantee the accurate capture of the phenomenon cannot be assessed rigorously in most of finite element applications. To choose these two parameters, rules of thumb are often used. In the case of the non-periodic homogenization, we know that our outputs $$\rho ^{\varepsilon _0\star }$$ and $${\bf C}^{\varepsilon _0\star }$$ only contain scales larger than λ0, so we sample them using a space-step $$dx = \frac{\lambda _0}{4}$$. This choice corresponds to twice the Nyquist rate. Because $${\bf C}^{\varepsilon _0\star }$$ is obtained from the space-derivative of the solution χs of the finite element analysis, we required an additional degree for this solution, so we want χs to be sampled by a space-step $$\delta x = \frac{\lambda _0}{5}$$. Such a sampling is ensured by imposing   $$\frac{a}{N} \leqslant \frac{\lambda _0}{5},$$ (E1)where a is the length of the edges of the mesh and N is the interpolation degree. eq. (E1) is our heuristic for sampling a model to be homogenized by our finite element code. In 3-D, assuming the tetrahedral elements to be regular, we can rewrite this heuristic in terms of the volume of the elements v given a degree N:   $$v \leqslant \frac{\sqrt{2}}{12} \left( \frac{N\lambda _0}{5} \right)^3.$$ (E2) © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations