Discrete wave equation upscaling

S U M M A R Y We present homogenization technique for the uniformly discretized wave equation, based on the derivation of an effective equation for the low-wavenumber component of the solution. The method produces a down-sampled, effective medium, thus making the solution of the effective equation less computationally expensive. Advantages of the method include its conceptual simplicity and ease of implementation, the applicability to any uniformly discretized wave equation in 1-D, 2-D or 3-D, and the absence of any constraints on the medium properties. We illustrate our method with a numerical example of wave propagation through a 1-D multiscale medium and demonstrate the accurate reproduction of the original wavefield for sufficiently low frequencies.


I N T RO D U C T I O N
Small-scale structure affects waves in a similar way as an effective structure with smooth variations.This fortunate fact allows us to understand and model their propagation without knowing or describing subwavelength details.The relation between the detailed medium and the coarse effective version is the subject of homogenization or upscaling theory.In the context of seismic wave propagation, it has been used to study apparent anisotropy induced by fine-scale heterogeneities (e.g.Backus 1962;Levshin & Ratnikova 1984;Fichtner et al. 2013), and to coarsen numerical meshes to reduce computational costs (e.g.Fichtner & Igel 2008;Capdeville et al. 2010b).Most upscaling methods employ a separation of fast and slow variables (Backus 1962), possibly in addition with a scale expansion (e.g.Capdeville et al. 2010aCapdeville et al. , 2015) ) or the Born approximation (Jordan 2015).They operate on the level of the partial differential equation (PDE), that is, on the differential form of the wave equation, given in the frequency domain and in 1-D by with position x, circular frequency ω, displacement u, external force f, density ρ and shear modulus μ.Our objective is to develop a discrete upscaling formulation that is conceptually simpler than PDE-based approaches, applicable to any type of wave equation, and free of restrictions on the medium properties.For this, we work directly on a suitably discretized version of (1), generically written in the form In eq. ( 2), u and f denote the discretized N-dimensional wavefield and force vectors, respectively.The N × N matrix D approximates a spatial derivative, and the sampled material properties ρ and μ are contained in the diagonal N × N matrices R and M, respectively.The discretizations of all quantities are assumed to be equidistant.In the presence of small-scale heterogeneities, the number of grid points, N, may need to be chosen much larger than required by the wavelength of interest, in order to capture all structural details.
To reduce the number of required grid points and the computational costs, we seek in the following an effective wave equation for the long-wavelength part of u that (i) has the same functional form as the discrete wave eq. ( 2), (ii) only depends on effective material parameters without small-scale detail and (iii) is sampled on K < N grid points.

D I S C R E T E U P S C A L I N G
We adopt the strategy of Hanasoge (2016) and start by transforming eq. ( 2) into the wavenumber domain using the discrete N-point Fourier transform

354
A. Fichtner and S.M. Hanasoge Multiplying eq. ( 2) by F N and inserting the identity matrix Defining the projections â = F N a and Â = F N AF −1 N for a generic vector a and a generic matrix A, and expanding eq. ( 3), we obtain (4) Eq. ( 4) is the wavenumber-domain version of the discrete wave eq. ( 2).The matrix D is nearly diagonal with components D jk ≈ iκ δ jk , where the continuous wavenumber κ is related to the grid spacing x via κ = 2π k/(N x).The solution of eq. ( 4) is the vector û containing all wavenumbers of the displacement field from 0 to N − 1.Our goal is to derive an equation for only the low wavenumbers of u, from 0 to some K < N − 1.For this, we decompose the N × N Fourier transform F N into a K × N matrix L that yields the low wavenumbers from 0 to K − 1, and an (N − K) × N matrix H that yields the high wavenumbers from K to N − 1.Using this definition, we deduce the following relation between L and H and their Hermitian conjugates L H and H H : We now repeat the steps that led from eq. ( 2) to (4), but using eq.( 5) instead of Defining the low-and high-wavenumber projections âL = La and âH = Ha for a generic vector a, and ÂLL = LAL H , ÂLH = LAH H , ÂHL = HAL H , and ÂHH = HAH H for a generic matrix A, simplifies eq. ( 6) to where we also used the approximate diagonality of D. Since off-diagonal elements of D are the result of boundary conditions and naturally imperfect finite-differences, this approximation corresponds to the assumption of a full space and perfect finite-differences.Eq. ( 7) represents two coupled vector-matrix equations for the low-and high-wavenumber components of the wavefield, ûL and ûH , respectively.Solving for ûL yields with the effective source Eqs ( 8) and ( 9) are separate equations for the low-wavenumber displacement field ûL .Their solution is equivalent to first solving eq. ( 4) for the complete û and then separating the low wavenumbers.While the solution of eq. ( 8) is the wavefield ûL , it is not a wave equation.Time derivatives represented by ω 2 are mixed with spatial derivatives, and products of density and the elastic parameter appear.However, eq. ( 8) approximates a wave equation when ω is sufficiently small.To see this, we estimate the L 1 -norm of DLL and DHH : where we assumed that the cut-off wavenumber K is significantly smaller than the maximum wavenumber N. Eq. ( 10) implies that the highwavenumber projection DHH is numerically much larger than its low-wavenumber counterpart DLL .It follows that small enough frequencies ω can be chosen such that terms involving ω 2 are insignificant in terms II, III and IV of eq. ( 8) that involve DHH .Even lower frequencies would be needed to eliminate ω 2 RLL from term I because DLL MLL DLL may not be dominant when the user-defined K is small.Following this argument, we omit ω-dependent expressions in terms II, III and IV, which condenses eq. ( 8) to The simplification from eq. ( 8) to ( 11) is ultimately based on a plausibility argument for low enough frequencies.Since the significance of the ω-dependent terms in (8) also depends on the properties of the medium that we wish to homogenize, the precise meaning of 'low enough' can only be assessed on a case-by-case basis-either by comparison to frequency-independent terms, or by numerically solving the wave equation for the original and the upscaled media.The latter will be done in the examples of Section 3. Rearranging terms in (11) leads to The quantities in eq. ( 14) are either K-dimensional vectors or K × K matrices.Thus, to return to the space domain, we insert K × K identity matrices I K = F K F −1 K expressed in terms of K-point Fourier transforms into eq.( 13): The application of F −1 K to ûe = ûL has the effect of a spatial resampling that transforms the low-wavenumber part of the original N-dimensional displacement into a smaller K-dimensional effective wavefield u e .Similarly, R e and M e are effective, that is, coarsened, material parameters.Omitting the leading Fourier transform of each term gives the effective version of the full space-domain wave eq. ( 2 (16) Eqs ( 13) and ( 15) reveal that upscaling density corresponds to the naive procedure of low-pass filtering (R → RLL = Re ) followed by downsampling ( Re → F −1 K Re F K = R e ).Upscaling the elastic modulus, in contrast, requires the additional correction term MLH M−1 HH MHL in eq. ( 13), called the corrector in classical homogenization theory (e.g.Capdeville et al. 2010aCapdeville et al. , 2015)).
Using the definitions of effective quantities from eq. ( 15), we can study upscaling for the simplest case of a homogeneous medium: Letting the elastic parameter M be the identity matrix Similarly, we find R e = I K for an originally homogeneous density R = I N .In accord with our intuitive expectation, upscaling a homogeneous medium with N grid points produces a homogeneous medium with K < N grid points.
More complex scenarios are most conveniently studied numerically, which we will do in the following section, where we will find that the matrices R e and M e representing the effective medium may become non-diagonal, meaning that the effective wave eq. ( 14) has a non-local rheology.Stress M e D e u e at a given grid point then depends on the strain De ûe at and around the grid point.

N U M E R I C A L E X A M P L E S A N D F U RT H E R S I M P L I F I C AT I O N S
We illustrate our developments with a model sampled at N = 3000 grid points with a spacing of x = 1 m.The model, shown in Fig. 1(a), features oscillations with 4 m period, a single 2 m wide spike and non-periodic variations with length scales ranging from 10 to 100 m.The upscaled model is sampled at K = 750 = 3000/4 grid points with spacing x e = 4 m.It results from extracting the diagonal elements of the non-diagonal matrices R e and M e , that is, by simplifying a non-local to a local rheology.The non-diagonal M e is visualized in Fig. 1(b).In the discrete effective medium, the periodic oscillations are replaced by constant properties, the spike is widened and reduced in amplitude, and sharp jumps appear smoothed.To assess upscaling quality, we compare wavefields computed by time-domain finite-differencing in the original 3000-point multiscale medium and its smooth 750-point version.We approximate the spatial derivatives in eq. ( 1) by the second-order stencil (e.g.Moczo et al. 2014) with the discrete stress σ i .First-order Neumann conditions are implemented at the boundaries using the stress imaging method (e.g.Levander 1988).Eq. ( 17) together with the boundary condition defines the discrete derivative operator D. The upscaled derivative operator D e is a densely populated matrix, meaning that the computation of a derivative at grid point i involves all K grid points, instead of only 2 as in eq. ( 17).The exact implementation of the upscaled eq. ( 16) would therefore lead to a very inefficient finite-difference method.Instead, we apply a further simplification and re-use the original finite-difference approximation (17), merely replacing x with x e .While this introduces another level of approximation, it leads to a numerical solution that is significantly less expensive; a factor 16 in our case because both the spatial and the temporal grid spacing are multiplied by 4.
In the interest of simplicity, we place the source at x = 400 m where the medium is locally homogeneous, thus avoiding the need to upscale the right-hand side.Original and upscaled wavefields agree well below ∼40 Hz, as exemplified by snapshots for various peak frequencies of the radiated Ricker wavelet (Figs 2a-e).This visual impression is quantified by the upscaling error in the form of the normalized L 2 misfit (Fig. 2f), that has various contributions: (i) the assumption of a diagonal D, (ii) the low-frequency approximation, (iii) the diagonalization of the effective modulus, and (iv) the increased grid spacing.Since the precise quantification of these contributions and their convergence properties requires lengthy and case-specific analytical efforts, we suggest to assess the quality of the upscaling numerically for concrete applications.In our example, Fig. 2(f) suggests that the increased grid spacing and the resulting insufficient number of grid points per wavelength ( 10) are dominant.The effect of omitting the correction term MLH M−1 HH MHL in eq. ( 13) in the upscaling of the elastic modulus-that is, replacing correct by naive upscaling in the form of low-pass filtering followed by downsampling-is illustrated in Fig. 3.The naively upscaled modulus is generally larger than the correctly upscaled modulus.While the corresponding velocity differences are only at the order of 1 per cent relative to the original medium, the resulting waveform differences are large, locally exceeding a quarter of a wavelength.

D I S C U S S I O N A N D C O N C L U S I O N S
We proposed a novel homogenization procedure for the wave equation, based on its equidistant discretization and the algebraic isolation of low-wavenumber components.Advantages of the method include its conceptual simplicity, the applicability to any type of wave equation and any type of medium.While we develop the concept in 1-D, it may be extended to higher dimensions using higher-dimensional Fourier transforms and modern parallel matrix solvers to invert the potentially large MHH .
Our method is similar to wavelet-based homogenization (e.g.Brewster & Beylkin 1995;Dorobantu & Engquist 1998;Engquist & Runborg 2002) that constructs a coarse-grid version of the operator L = −ω 2 R − DMD from eq. ( 2).In addition to using Fourier instead of wavelet transforms, our method is different mostly through the explicit incorporation of the finite-difference matrix D. This allows us to (i) construct an effective medium and an effective operator, (ii) use this medium with standard finite-difference stencils that operate on few grid points, and (iii) recover an effective wave equation that can be solved with the same code used for the original medium.While the use of the Fourier transform is largely responsible for our method's simplicity, it also imposes its greatest disadvantage, namely the limitation to equidistant discretizations.Work to eliminate this restriction is in progress.Engquist & Runborg (2002) noted similarities between wavelet-based homogenization and two-scale homogenization on the continuous PDE level (e.g.Capdeville et al. 2010aCapdeville et al. , 2015) ) that also translate to our method.The homogenized coefficients are generally a sum of two contributions: (i) a simple arithmetic average or downsampling, and (ii) a corrector for the interaction between high and low wavenumbers.
An interesting difference between two-scale homogenization and our method is the explicit appearance of a non-local rheology in the effective eq.( 16); well-known in studies of elastic composites and granular flow (e.g.Drugan & Willis 1995;Pouliquen & Forterre 2009).Obeying an uncertainty principle, the material properties de-localize and expand in the space domain as more and more components are omitted in the wavenumber domain (Hanasoge 2016).To obtain a numerically tractable scheme, the non-local rheology needs to be localized, which involves a diagonalization of the upscaled material properties R e and M e .We diagonalized by extracting the diagonal and omitting the off-diagonal terms.At least in the examples considered here, this simplistic approach produced significantly more accurate results than mass-lumping off-diagonal into diagonal elements.As in finite-element analysis, the effect of mass lumping is hard to predict and not generally positive.Finally, we note that our notion of non-locality should not be confused with the one used in multi-scale homogenization where it refers to the appearance of higher spatial derivatives (e.g.Fish et al. 2002).

A C K N O W L E D G E M E N T S
The authors would like to thank Michael Afanasiev, Yann Capdeville, Joerg Renner and an anonymous reviewer for discussions and suggestions on how to improve the manuscript.This study received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement number 714069).

Figure 2 .Figure 3 .
Figure 2. Wavefield snapshots and convergence.(a-e) Snapshots of the original wavefield (black) on 3000 grid points and the upscaled wavefield (red) on 750 grid points for various propagation times t and peak frequencies f peak .Dots mark every 100th grid point.(f) L 2 distance between original and upscaled wavefields for propagation time t = 1.0 s, normalized by the maximum amplitude.The peak frequency f peak and the minimum number of grid points per wavelength (g.p.p.w.) are on the horizontal axes.(Colour version available online.)