Yann Capdeville
email: yann.capdeville@univ-nantes.fr

Homogenization of seismic point and extended sources

Keywords: Theoretical seismology, Computational seismology, Earthquake ground motions, Wave propagation, Numerical modelling

Seismic sources are mostly modelled as point sources: moment tensors associated with the gradient of a Dirac distribution. Such sources contain an infinite range of scales and induce a discontinuity in the displacement wavefield. This makes the near-source wavefield expensive to model and the event location complex to invert, in particular for large events for which many point sources are required. In this work, we propose to apply the non-periodic two-scale homogenization method to the wave equation source term for both force and couple-sources. We show it is possible to replace the Dirac point source with a smooth source term, valid in a given seismic signal frequency band. The discontinuous wavefield near-source wavefield can be recovered using a corrector that needs to be added to the solution obtained solving the wave equation with the smooth source term. We show that, compared to classical applications of the two-scale homogenization method to heterogeneous media, the source term homogenization has some interesting particularities: for couple-sources, the leading term of the homogenization asymptotic expansion is dependent on the fine spatial scale, depending on the source type, only one or two first terms of the expansion are non-zero and there is no periodic case equivalent (the source term cannot be made spatially periodic). For heterogeneous media, two options are developed.

In the first one, only the source is homogenized while the medium itself remains the same, 2 Y. CAPDEVILLE including its discontinuities. In the second one, both the source and the medium are homogenized successively: first the medium and then the source. We present a set of tests in 1-D and 2-D, showing accurate results both in the far-source and near-source wavefields, before discussing the interest of this work in the forward and inverse problem contexts.

INTRODUCTION

Seismic waves are excited by atmospheric turbulences, landslides, meteorite impacts, natural and artificial explosions and earthquakes. Most of these sources include complex processes but they are often mathematically represented by simple point sources [START_REF] Aki | Quantitative Seismology: Theory and Methods[END_REF]. This approximation is widely used in seismology and is valid as long as the minimum wavelength is large compared to the spatial extent of the event. Two types of point sources are commonly used:

• single force:

f (x, t) = Fδ(x -x 0 )g(t)
(1)

• double-couple:

f (x, t) = -M • ∇δ(x -x 0 )g(t) (2)
where f is the force vector to be used in the wave equation, δ(xx 0 ) the is Dirac distribution centered on x 0 , F is a force vector, M the moment tensor density, ∇ the gradient operator and g(t)

the source time function. Classical earthquakes correspond to double-couple sources and more exotic sources, such as landslides, to single force sources. More complex sources can also be treated in that framework using a sum over many point sources, with the notion of moment tensor density [START_REF] Aki | Quantitative Seismology: Theory and Methods[END_REF][START_REF] Madariaga | 4.02 -seismic source theory[END_REF]:

f (x, t) = S m(x , t) • ∇δ(x -x ) dx , ( 3
)
where S is the fault surface and m is the moment tensor density on the fault. The moment tensor density can be related to the displacement jump through the fault [u] using the following relation

m ij (x, t) = c klij n k [u l ](x, t) , ( 4
)
where c is the elastic tensor and n is normal vector to the fault surface. Using a numerical quadrature, the integral in (3) can be expressed as a sum over many points, leading to a sum over distributed point sources.

Effective seismic sources 3

The spatial extent of these point sources is assumed to be infinitely small compared to the farfield minimum wavelength associated with the maximum frequency of the source wavelet g(t). This difference of spatial scales between the source size and minimum wavelength of the wavefield is often not considered as a serious difficulty for both forward and inverse problems in seismology. Source implementation is not difficult for most forward modelling solvers and source inversions are routinely done by most seismological observatories. However, as we are about to show, this requires a closer look.

In the following, we use the notion of near-source field and far-source field. The near-source field is here the displacement field, solution of the wave equation, close to the point source location within a distance smaller than a fraction of the minimum wavelength. It is different from the near-field as defined by [START_REF] Aki | Quantitative Seismology: Theory and Methods[END_REF] which is the static part of the soil response to the source. The near-field is not necessarily near the source. The two fields have nevertheless a connection as the near-source field is dominated by the near-field part of the wave equation solution.

In the forward modelling context, to discuss the implementation of point sources, we need to first distinguish strong-form from weak-form of the wave equation. For methods based the strongform of the wave equation, such as the finite difference (FD) method, point sources are problematic.

This difficulty comes from the fact that point sources are defined with a Dirac distribution. Here, the mathematical specificity of distributions is important: Dirac distributions are generalized functions that only make sense once integrated against any regular function (a test function). Nevertheless, for most strong-form methods, the point source is seen as a regular function. If the source is located on a grid point, this has little consequences and the source force is added to the stress only at this specific grid point [START_REF] Alterman | Source of finite extent, applied force and couple in an elastic half-space[END_REF][START_REF] Virieux | P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF][START_REF] Coutant | Numerical source implementation in a 2d finite difference scheme for wave propagation[END_REF]. This solution is not available when the source is not located on a grid point, or when the mechanical properties at the source location are heterogeneous. For the issue related to the source location, a classical solution is to use initial conditions computed from an analytical solution in a small area around the source [START_REF] Alterman | Propagation of elastic waves in layered media by finite difference methods[END_REF]. This method suffers from several issues: the source cannot be located too close to material discontinuities and an analytical solution must exist, limiting this solution to very simple media in the source area. Moreover, backscattered energy passing through the source region can also be a problem. Nevertheless, when this solution is possible, it has the advantage to be accurate also near the source and not only in the far-field. The other solution is to extend the spatial extent of the source to several grid points [START_REF] Mittet | General source and receiver positions in coarse-grid finitedifference schemes[END_REF][START_REF] Hicks | Arbitrary source and receiver positioning in finite-difference schemes using kaiser windowed sinc functions[END_REF], basically by spatially low-pass filtering the Dirac distribution. This solution also requires the medium mechanical properties to be smoothly varying near the source and the wavefield is only accurate in the far-field.

For methods based on the weak-form of the wave equation, such as, for example, the normal mode summations method (NMS) [START_REF] Gilbert | An application of normal mode theory to the retrieval of structural parameters and source mechanisms for seismic spectra[END_REF] or the spectral element method (SEM) [START_REF] Komatitsch | The spectral element method: an effective tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Chaljub | Spectral element analysis in seismology[END_REF], point sources are less of a problem. For such methods, the source term contributes through the integral of the product of test functions with the point source. In this case, the Dirac distribution is used as it should be, i.e. integrated against a smooth function. The source term can therefore easily and naturally implemented into any weak-form method.

Nevertheless, it has been observed that, while the solution is accurate in the far-field, it is not the case near the source. For the NMS, for example, the test functions are based on spherical harmonics and have a global support at the whole earth scale. NMS is truncated based on the maximum frequency of the source, making it possible to accurately compute the far-field for a low numerical cost. But in the near-source field, the solution is not accurate. Depending on how the truncation is performed, the source can spread all over the sphere and damage the solution for the source duration time because the test functions are global. Similarly, for the SEM, it is been observed that, in the source element, the discretized source term spreads all over the element, leading to an inaccurate solution in this element [START_REF] Faccioli | 2d and 3d elastic wave propagation by a pseudo-spectral domain decomposition method[END_REF][START_REF] Chaljub | Modèlisation numérique de la propagation d'ondes sismiques à l'échelle du globe[END_REF][START_REF] Nissen-Meyer | A two-dimensional spectral-element method for computing spherical-earth seismograms-i. moment-tensor source[END_REF][START_REF] Fichtner | Full seismic waveform modelling and inversion[END_REF][START_REF] Igel | Computational seismology: a practical introduction[END_REF]). This effect is often explained by the fact that the Dirac source cannot be accurately represented on the test function basis, and consequently leads to a source that is not a point source in practice. However, it is rather the wavefield that cannot be represented on the test function basis near the source. Indeed, the wavefield is singular at the source location and test functions do not allow such a discontinuity. Interestingly, the missing near-source field part does not carry energy and this problem does not impact the solution away from the source point. If one needs an accurate solution near the source, the solution is to use a denser mesh (for SEM) or a higher spherical harmonic angular order truncation (for NMS), keeping in mind the solution will still be inaccurate close to the source, but in a smaller region. The drawback of such a solution is a higher numerical cost. In the end, for both forms, a point source is a difficulty to which we look for a general solution.

The inverse problem for point sources is also a very important topic in seismology. The centroid moment tensor (CMT) is obtained with such an inverse problem [START_REF] Dziewonski | Determination of earthquake source parameters from waveform data for studies of global and regional seismicity[END_REF]). If the moment tensor inversion is not a problem, the inversion of the event location is: it is strongly non-linear and a good initial guess is necessary. For large sources, the point source approximation is challenged and the extension from a point source to something more general is not trivial. The existing solutions, such as higher order moment tensors [START_REF] Clévédé | New constraints on the rupture process of the 1999 august 17 izmit earthquake deduced from estimates of stress glut rate moments[END_REF][START_REF] Jordan | Representation of complex seismic sources by orthogonal moment-tensor fields[END_REF] or finite source inversions [START_REF] Mai | The earthquake-source inversion validation (siv) project[END_REF], are difficult to implement in tomographic full waveform inversion in a automated way for example.

These difficulties in forward and inverse problems related to point sources are similar to the difficulties encountered for small-scale present in the earth model: inhomogeneities of size much smaller

Effective seismic sources 5 than the minimum wavelength make the forward modelling very expensive and the inverse problem unstable when we attempt to retrieve them from seismic data. These can be addressed with the notion of effective media and homogenization for both the forward (Capdeville et al. 2010b) and the inverse problem [START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-d numerical illustrations[END_REF]. In the context of the homogenization technique, the seismic point source interaction with nearby small-scale heterogeneities has been dealt with in previous works (Capdeville et al. 2010b;[START_REF] Burgos | Homogenized moment tensor and the effect of near-field heterogeneities on nonisotropic radiation in nuclear explosion[END_REF]) but without homogenizing the source. Indeed, in those works, the effect of homogenizing the media on the point source is accounted for through an effective moment tensor but the source is still a point source. The idea of present work is to apply the nonperiodic homogenization method to point sources, allowing one to find an equivalent smooth effective source and a local corrector to model the sharp near-source field, both valid for a limited frequency

band. An effective source replacing a point source, or an ensemble of point sources of an extended fault system, would solve all the difficulties mentioned above, both in the forward modeling and the inverse problem contexts.

The paper is organized as follows: in the first part, the theory and examples are developed in 1-D.

In the second part, the higher dimension theory is given and some examples in 2-D are shown. We then discuss the results and conclude our observations.

1-D THEORETICAL DEVELOPMENT

In this part, we consider a wave in a 1-D elastic bar of length L with free boundaries at both ends. The displacement along the bar u(x, t) is driven by the equations

ρ∂ tt u -∂ x σ = f , (5)
σ = E∂ x u -τ , (6)
where x is the position in the bar, ρ(x) is the density along the bar, E(x) the elastic property, and σ(x, t) the stress. ρ and E are both positive quantities. The free stress boundary condition at both ends imply ∂ x u(0, t) = ∂ x u(L, t) = 0. The wave speed in the bar is α(x) = E(x)/ρ(x). The source terms f (x, t) and τ (x, t) are point sources:

• single force:

f (x, t) = F δ(x -x 0 )g(t) , (7)
• couple-force or stress-source:

τ (x, t) = M δ(x -x 0 )g(t) , (8)
where F and M are scalars related to the amplitude of the source. Note that the stress source term is often inputted in the dynamic equations (5) using a force term f τ (x, t) = -∂ x τ (x, t) such that

ρ∂ tt u -∂ x σ = f -∂ x τ (x, t) , (9)
σ = E∂ x u , (10)

Numerical preliminary observations

Before starting the core theoretical development of this paper, it is useful to visualize some examples of wavefields in the vicinity of a point source. To do so, we choose a very heterogeneous bar of length 40 m made up of 1000 material pieces of constant velocity and density materials. In each piece, the mechanical properties are generated randomly (but only once). A sample of the velocities in the bar is displayed in Fig. 1a. The source time function g(t) is a Ricker wavelet (second derivative of a

Gaussian) with a central frequency of 0.5 kHhz (maximum frequency about 1.5 kHz) leading to an averaged minimum wavelength (λ min ) of 3.3 m in the far-field. This wavelength is much larger than the individual bar pieces, which are of length 0.04 m. The source position is x 0 = 15 m. We solve the wave equation (5-6) numerically using SEM. To obtain an accurate solution, each velocity and density discontinuity is honored by an edge of an element mesh. Doing so, the mesh is very fine compared to a mesh that would otherwise be used in a homogeneous bar. With such a mesh, the solution can be considered as a reference solution everywhere except at the element containing the source. If the source is near an element edge, the solution can also be less accurate in the element right next to it. This is expected: knowing that the solution u has an order 0 or 1 discontinuity (depending on the type of source), it cannot be accurately represented on the SEM Lagrange polynomial basis. Interestingly, as mentioned in the introduction, the fact that the solution is not accurate in the element of the source does not harm the quality of the solution in the other elements. This is one of the interesting properties of SEM (and any weak-form numerical method): because it is based on an energy formulation, the wavefield is accurately computed almost-everywhere, even if it is not true in the element containing the source.

For the couple-source (8), Fig. 1b displays the snapshot for displacement response u(x, t) at a time t taken just after the source time function wavelet g maximum. It shows a sharp discontinuity at the source position and a smooth variation otherwise. As expected, the displacement is discontinuous at the source location. Of course, a real discontinuity cannot be accommodated by SEM (the displacement is explicitly continuous in the method formulation) and this creates an error in the element containing the source. This error is poorly visible at the scale of the plot, but a small glitch can be observed at the lower part of the displacement curve at x = 15 m. The corresponding strain (Fig. 1c)

shows a singularity at the source location. The imprint of the elastic structure is visible, as expected (the strain is discontinuous in discontinuous models [START_REF] Capdeville | Chapter six -an introduction to the two-scale homogenization method for seismology[END_REF]).

For the single force (7), Fig. 1d displays the displacement response just after the source wavelet maximum. It shows a kink (discontinuity in the derivative) at the source position and a smooth variation otherwise. The corresponding strain (Fig. 1e) shows a step function at the source location. The imprint of the elastic structure is once again clearly visible.

From the point of view of numerical modelling, this example shows two difficulties: the small size of heterogeneities and the discontinuity of the solution at the source. The size of the heterogeneities is a problem because it implies a fine mesh, here about 90 elements per λ min where, in a smooth medium, a single degree 5 element would have been enough. This leads to a high numerical cost.

Removing the small-scale using homogenization is a solution to this difficulty. The homogenization technique is presented in the next section. For weak-form methods, the point source discontinuity is not seen as difficulty (because it does not need any specific treatment), unless one is interested in the solution very close to the source. In this case, the only option so far is to have a fine mesh near the source, also leading to a high numerical cost. In the following, we will show that the source can also be homogenized and this can help remove this difficulty.

Two-scale homogenization of the mechanical properties with no scale separation

Two-scale homogenization is a good mathematical framework to deal with scale problems. It has been developed for periodic and stochastic problems (Sanchez-Palencia 1980) and has been extended

to multi-scale heterogeneous media for wave propagation (Capdeville et al. 2010b;[START_REF] Cupillard | Non-periodic homogenization of 3-d elastic media for the seismic wave equation[END_REF]. In this section, we summarize the results of Capdeville et al. (2010a) on the twoscale homogenization method applied to the elastic wave equation when the mechanical property heterogeneities do not present any natural scale separation. Many mathematical subtleties have been omitted here but can be found in Capdeville et al. (2010a) and [START_REF] Capdeville | Chapter six -an introduction to the two-scale homogenization method for seismology[END_REF].

Two-scale homogenization is a mathematical process that makes it possible to explicitly extract large-scale (effective) effects from small-scale ones. It is an asymptotic theory based on a small parameter

ε 0 = λ 0 λmin , (11)
where λ min is the minimum wavelength and λ 0 is a scale below which scales in the medium are considered as small. The 0 underscore of ε 0 is meant to differentiate it from the ε of classical periodic two-scale homogenization which is tied to the periodicity of the heterogeneities. In the general case, where periodicity is not assumed, λ 0 is a user-defined scale, usually set to be smaller than λ min . To separate the scales, a low-pass filter F ε 0 is introduced. In 1-D, when applied to a function h, it can be expressed as a convolution with the filtering wavelet w ε 0 ,

F ε 0 (h) = R h(x )w ε 0 (x -x ) dx . (12)
F ε 0 is a linear spatial filter designed to mute spatial frequencies larger than (ε 0 λ 0 ) -1 of the space function it is applied to, to zero. The design of the wavelet w ε 0 is not unique. It is built as a compromise

Effective seismic sources 9 between a spatial extent support as small as possible and a spatial frequency cut-off as sharp as possible while keeping low-spatial frequencies intact. In the present work, we use simple cosine taper filters.

Note that we assume here that λ min remains roughly spatially constant through the domain. This is often not the case in practice as wave-speeds strongly vary with space in geological media, at least with depth. In such cases, one can always use the smallest available λ min , but it is not optimal and allowing changes in λ 0 through space might be desirable. A solution toward this possibility valid for layered media is given in Capdeville et al. (2013), appendix B. For 2-D and 3-D smooth variations, a simple variable filtering according to a assumed smooth variation of the dominant minimum wavelength is often good enough. Nevertheless, a complete work giving a general solution to this limitation has yet to be written.

Two-scale homogenization relies on the introduction of two space variables instead of one: x for the large-scale variations and y for small-scale variations. y is defined as

y = x ε 0 . ( 13
)
In the following, x and y are mathematically treated as independent variables. Knowing the y definition above, making x and y independent can be surprising. Nevertheless, the two-scale problem based on these two variables that we are about to describe can be shown to converge to the effective one scale problem, at least in the periodic case [START_REF] Sanchez-Palencia | Non homogeneous media and vibration theory[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF].

Because of the definition of y (13), any partial derivative with respect to the space variable needs to be changed according to:

∂ x → ∂ x + 1 ε 0 ∂ y . (14)
Before going further, we need to introduce F, the same low-pass filter than F ε 0 but applied to the variable y instead of x. It means that, for any function h(x) and h(y) = h(ε 0 y), we have

F ε 0 (h)(x) = F( h)(x/ε 0 ) . ( 15
)
F wavenumber cutoff is always 1/λ min . More details about this technical aspect can be found in [START_REF] Capdeville | Chapter six -an introduction to the two-scale homogenization method for seismology[END_REF], section 2.3.

The last part of the two-scale homogenization is the asymptotic expansion of the wave equation solution as a power series of ε 0 :

u(x, t) = u 0 (x, y, t) + ε 0 u 1 (x, y, t) + ε 2 0 u 2 (x, y, t) + ... (16) σ(x, t) = σ 0 (x, y, t) + ε 0 σ 1 (x, y, t) + ε 2 0 σ 2 (x, y, t) + ... (17)
An important constraint on u i , σ i is that they must contain only small-scale variations on the y axis.

This translates into the fact that if we apply F to any coefficients u i , σ i , we obtain constants in y. This can mathematically expressed as u i , σ i must belong to the space V, defined as

V = {h(x, y), Y -periodic in y such that F(h)(x, y) = h (x)} , (18)
where for any function h(x, y)

h (x) = 1 Y Y 0 h(x, y)dy , ( 19
)
and Y is segment along the y axis, wide enough to fit the spatial support of F. In the following, we will make use of V * , a variant of V defined as

V * = {h(x, y), Y -periodic in y such that F(h)(x, y) = 0} . ( 20
)
The next step is, starting from ρ and E, to build quantities ρ ε 0 (x, y), E ε 0 (x, y) such that the coefficients u i , σ i exist in V. This is the main difficulty of the two-scale homogenization when the periodicity assumption is not made. The complete procedure to do so is described in Capdeville et al. (2010a, Capdeville et al. (2010b).

Introducing ( 14) and (16-17) in (5-6), we find the equations driving the expansion coefficients u i , σ i :

ρ ε 0 (x, y)∂ tt u i (x, y, t) -∂ x σ i (x, y, t) -∂ y σ i+1 (x, y, t) = f δ i,0 , (21)
σ i (x, y, t) = E ε 0 (x, y)(∂ x u i (x, y, t) + ∂ y u i+1 (x, y, t)) . (22)
We have two types of boundary conditions: the boundary condition along the x axis which, for each i, σ i (0, y, t) = σ i (L, y, t) = 0, and periodic boundary conditions for any equation along the y axis.

Solving (21)(22) for consequent indices i up to the order 1, displacement can be found to be written as

u(x, t) = u 0 (x, t) + ε 0 χ ε 0 (x, y)∂ x u 0 (x) + O(ε 2 0 ) , (23)
where y = x/ε 0 (it has only one value here). To the leading order, the strain is

(x, t) = ∂ x u(x, t) = (1 + ∂ y χ ε 0 (x, y)) ∂ x u 0 (x, t) + O(ε 0 ) . ( 24
)
It can be noted that displacement to the leading order, u 0 , does not depend on the small-scale variable y. It can also be shown that σ 0 does not depend upon y. χ ε 0 is the first order corrector, solution of the cell problem

∂ y [E ε 0 (1 + ∂ y χ ε 0 )] = 0 (25)
on Y with periodic boundary conditions. u 0 , σ 0 are solutions of the order 0 effective wave equation:

ρ * ,ε 0 ∂ tt u 0 -∂ x σ 0 = f , (26)
σ 0 = E * ,ε 0 ∂ x u 0 , (27)
where ρ * ,ε 0 = ρ ε 0 and

E * ,ε 0 (x) = E ε 0 (1 + ∂ y χ ε 0 ) (x) . (28)
Luckily, it can be seen that the effective equation is still a classical wave equation. We can therefore still use our favorite wave equation solver to get a solution to this wave equation. Moreover, the effective parameters are free of small-scale, which makes the space discretization simple.

In the 1-D case, there is an analytical solution to the cell problem ( 25), that makes it possible to

show that

E * ,ε 0 (x) = F ε 0 1 E -1 (x) . (29) Moreover, ρ * ,ε 0 (x) = F ε 0 (ρ)(x).
Nothing has been done yet about the source term. It is usually ignored in the homogenization development and introduced back a posteriori, using an energy argument. It can be shown that, to the leading order, for a point single force, the source term can be used without any change and, for a couple-source, it needs to be modified according to

f (x, t) = M 1 + ∂ y χ ε 0 x 0 ε 0 ∂ x δ(x -x 0 )g(t) . (30)
In this 1-D case, due to the analytical solution for the cell equation, the couple-source force correction reduces to

f (x, t) = M E * ,ε 0 E (x 0 )∂ x δ(x -x 0 )g(t) . (31)
In both cases, the source term still involves a Dirac function and therefore is not free of small-scale.

Two-scale homogenization for the source term only

As shown in Capdeville et al. (2010a, Capdeville et al. (2010b), it is often very useful to introduce the concept of homogenization through the periodic case. This helps to better understand the homogenization process and grasp its concepts. Then, the more difficult non-periodic case can be considered.

Unfortunately, while homogenization of the mechanical properties can be derived both in periodic and non-periodic cases, this is not the case for the homogenization of the source term. Indeed, the seismic source is intrinsically non-periodic: there is usually only one or a few earthquakes at the same time.

Moreover, building an artificial problem with finely space sources paving the bar makes no sense.

In this section, we assume that only the source term f contains small-scale that need to be homogenized. This assumption can be made if the mechanical properties (ρ(x), E(x)) are free of small-scale (if they are constant or smooth) but it can also be made if they contain small scales such as discontinuities. In the latter case, this assumption means we choose to homogenize the small-scale of the source while leaving those from the mechanical properties unchanged. This means that the small scales from the mechanical properties will remain and will have to be accounted for with a dense mesh while solving the wave equation with the effective source. This non-intuitive choice can be useful in situations where small-scale need to be kept in an earth model. For example, some medium can be simple to mesh but still have small-scale such as discontinuities. Layered media such as PREM are a good example of such models. In such cases, one might want to avoid having to homogenize the whole model to input an effective source. In those situations, the possibility to homogenize the source while keeping the fine-scale structures of the medium might be desirable. The more general case where both the source and the mechanical properties are homogenized is treated in section Sec. 2.4.

To obtain the effective source and the near-source field corrector, we follow the same procedure as the one used to homogenize the wave equation in the non-periodic heterogeneous mechanical properties case. It consists of three steps:

(i) build the homogenized equations;

(ii) solve the homogenized equations, order by order, assuming it is possible to build the source expansion coefficients such that the expansion coefficients of ( 16) and ( 17) belong to V * ;

(iii) build the two-scale expansion coefficients of the source term based on the conditions derived from the previous step.

We begin with the assumption that the two types of sources can be expanded as

f (x, t) = 1 ε 0 f -1 (x, y)g(t) + f 0 (x)g(t) , (32)
τ (x, t) = 1 ε 0 τ -1 (x, y)g(t) + τ 0 (x)g(t) , (33)
where -1 here means index i = -1 and where y is taken equal to x/ε 0 . This assumption comes from the idea the a Dirac can be decomposed into a smooth and high wavenumber parts:

δ(x) = F ε 0 (δ)(x) + (I -F ε 0 )(δ)(x) , (34)
which, using ( 15), can also be written:

δ(x) = F ε 0 (δ)(x) + (I -F)(δ)(ε 0 y) . (35)
Then, using the property

δ(ε 0 y) = 1 ε 0 δ(y) , (36)
a Dirac function can therefore be decomposed as

δ(x) = 1 ε 0 (I -F)(δ)(y) + F ε 0 (δ)(x) . ( 37
)
Effective seismic sources 13

The Dirac expansion series has therefore two terms, one in (ε 0 ) -1 and one in (ε 0 ) 0 . The fact that the two different types of sources contain a Dirac distribution naturally leads to (32) and ( 33). Note that we have assumed that τ 0 and f 0 depend only on x. This assumption will be proved by construction in Sec. 2.3.3. The expansion (37) of the Dirac distribution is the basis of the homogenized solution presented in this paper.

At this stage, we do not know yet how to build the coefficients f i and τ i . This will be done a posteriori.

The homogenized equations to be solved are therefore

ρ(x)∂ tt u i (x, y, t) -∂ x σ i (x, y, t) -∂ y σ i+1 (x, y, t) = f i (x, y)g(t) , (38)
σ i (x, y, t) = E(x)(∂ x u i (x, y, t) + ∂ y u i+1 (x, y, t)) -τ i (x, y)g(t) . ( 39
)
The boundary conditions are the same as for the homogenization of the mechanical properties.

In the following, for the sake of simplicity, the couple and single-force cases are solved separately.

Resolution for the couple-source case

We assume that the source is a couple-source force, f (x, t) = -∂ x τ (x)g(t). We solve the homogenized equations (38-39) one by one:

• Eqs (38) for i = -2 and (39) for i = -1 give

∂ y σ -1 = 0 , (40)
σ -1 = E∂ y u 0 -τ -1 . ( 41
)
Combining the last two equations, we have

∂ y E∂ y u 0 = ∂ y τ -1 . (42)
We define the leading order corrector θ τ as

u 0 (x, y, t) = u 0 (x, t) + θ τ (y)g(t) . ( 43
)
From the last equation, we see that having θ τ ∈ V * is a necessary and sufficient condition to have u 0 ∈ V. Using ( 43) in ( 42) we find that θ τ is solution of the following cell problem on Y :

E(x)∂ 2 y θ τ (y) = ∂ y τ -1 (x, y) , (44)
with periodic boundary conditions. As it will be seen in Sec. 2.3.3, θ τ is independent of x by construction. Moreover, from (41) and by construction of τ -1 (see Sec. 2.3.3, Eq. 74), we have

σ -1 = 0 . ( 45
)
It is interesting to note that, unlike most homogenization problems, the leading order of displacement u 0 depends on y.

• Eqs (38) for i = -1, (39) for i = 0 and (45) give

∂ y σ 0 = 0 , (46)
σ 0 = E(∂ x u 0 + ∂ y u 1 ) -τ 0 . ( 47
)
Using ( 43) in ( 47), knowing that θ τ only depends on y, we have

σ 0 = E(∂ x u 0 + ∂ y u 1 ) -τ 0 . ( 48
)
Taking the cell average of the last equation and using the fact that

∂ y h = 0 for any function h(x, y) , (49)
we find the order 0 effective constitutive equation:

σ 0 = E∂ x u 0 -τ 0 . (50)
• Eqs (38) for i = 0 gives

ρ(x)∂ tt u 0 -∂ x σ 0 -∂ y σ 1 = 0 , (51)
Taking the cell average of the last equation and using the property (49), we obtain the order 0 dynamic equation:

ρ(x)∂ tt u 0 -∂ x σ 0 = 0 , (52)
At this stage, the order 0 homogenized solution u 0 is found solving (50) and ( 52).

• going further in the order, for i > 0, taking the cell average of (38-39) and using the property (49) we find

ρ∂ tt u i -∂ x σ i = 0 , ( 53
)
σ i = E(x)∂ x u i . ( 54
)
These equations have a unique solution that is u i = σ i = 0. All the expansion coefficients larger than i = 0 are therefore null. This remarkable result is very different from the elastic homogenization for which only the order 1 solution is null. Its practical consequence on the convergence rate from the homogenized solution to the true solution as a function of ε 0 is shown in Sec. 2.5.1, in a homogeneous bar example.

Effective seismic sources 15

To summarize, to obtain the order 0 solution, one need to solve

ρ(x)∂ tt u 0 -∂ x σ 0 = 0 , (55)
σ 0 = E∂ x u 0 -τ 0 , (56)
To obtain the order 1, we have

u(x, t) = u 0 (x) + θ τ (y)g(t) , (57)
where all the expressions are evaluated for y = x/ε 0 . Note that no O(ε 2 0 ) is present in these expressions: this relation is valid for all orders.

For the order 0 strain, we have

(x, t) = 1 ε 0 ∂ y θ τ (y)g(t) + ∂ x u 0 (x, t) . ( 58
)
The order 0 corrector is solution of the cell equation on Y ,

E(x)∂ 2 y θ τ (y) = ∂ y τ -1 (x, y) , (59)
with periodic boundary condition.

We still need to define τ -1 and τ 0 so that the solutions to the above cell problem exist and belong to V * . This is the purpose of section 2.3.3 .

Resolution for the single force case

The resolution of the single force case is similar to the couple case presented in the previous section.

The details are given in Appendix A. We find the following order 0 effective wave equation:

ρ(x)∂ tt u 0 -∂ x σ 0 = f 0 , (60)
σ 0 = E∂ x u 0 . (61)
In contrast to the couple-source case, there is no corrector for the order 0 displacement (it only depends on x):

u 0 (x, y) = u 0 (x) . ( 62
)
Adding the order 1 corrector the order 0 solution, we have:

u(x) = u 0 (x) + ε 0 θ f (y) , (63)
where all the expressions are evaluated for y = x/ε 0 . The corrector is solution of

E(x)∂ 2 y θ f (y) = f -1 (x, y) , (64)
on Y with periodic boundary conditions.

Similarly to the previous section, we can show that all the higher-order solution coefficients are null.

We still need to define f -1 and f 0 so that the solutions to the above cell problem exist and belong to V * . That is the purpose of section 2.3.3 .

Construction of source expansion coefficients

We need to build the source expansion coefficients such that the displacement and stress expansion coefficients u i , σ i belong to V. In order to achieve this, the following sufficient conditions must be met:

(i) θ τ and θ f must be in V * (from Eqs. ( 43) and (A10));

(ii) the expansion ( 32) and ( 33) must be satisfied.

We follow a similar procedure to the one used for the mechanical properties (Capdeville et al. 2010a):

we first assume that all spatial variations are small-scale and then use the filter F to separate the scales on the correctors. Setting M = 1 and F = 1 in ( 7) and ( 8) respectively, omitting the source time dependence for now and transforming δ(x) in the small-scale variable domain using (36), we have

f -1,s (y) = δ(y) , (65)
f 0,s (x) = 0 , (66)
and

τ -1,s (y) = δ(y) , (67)
τ 0,s (x) = 0 . ( 68
)
We then solve a small-scale only version of the cell problems ( 44) and ( 64) on Y with periodic boundary conditions, assuming the variation of the elastic parameters are also on y, to obtain the starting correctors θ s f and θ s τ :

∂ y Ē(y)∂ y θ s σ (y) = ∂ y τ -1,s (y) , (69)
∂ y Ē(y)∂ y θ s f (y) = f -1,s (y) , (70)
where Ē(y) = E(ε 0 y). We then separate the scales to obtain the final correctors

θ i (y) = (I -F)(θ s i )(y) , i ∈ {f, τ } . ( 71
)
At this stage, by construction, both correctors θ i belong to V * , satisfying then condition (i). Moreover, both correctors only depend on y by construction.

Effective seismic sources 17

We also define the order 0 source potential from the smooth part of the starting correctors:

ψ i (x) = F(θ s i )(x/ε 0 ) , i ∈ {f, τ } . ( 72
)
The order -1 source terms are then built as

f -1 (x, y) = E(x)∂ 2 y θ f (y) , (73)
τ -1 (x, y) = E(x)∂ y θ τ (y) . ( 74
)
The last equation justifies (45). Finally, the large-scale variations of the correctors is used to define the order 0 source terms:

f 0 (x) = ∂ x E(x)∂ x ψ f (x) , (75)
τ 0 (x) = E(x)∂ x ψ τ (x) . ( 76
)
The last equations also show that, by construction, f 0 and τ 0 can indeed be built as independent of y.

Next, we can check that the condition (ii) is satisfied. For the case of couple-source, from ( 74) and ( 76), for y = x ε 0 , using ∂ y = ε 0 ∂ x we have

1 ε 0 τ -1 + τ 0 = E(x) ε 0 ∂ y θ τ (y) + E(x)∂ x ψ τ (x) , (77)
= E(x) ε 0 ∂ y θ τ (y) + E(x)∂ x (F)(θ s τ )(x/ε 0 ) , ( 78
) = E(x) ε 0 (∂ y θ τ (y) + ∂ y (F)(θ s τ )(y)) (79) = Ē(y) ε 0 (∂ y (I -F)(θ s σ )(y) + ∂ y (F)(θ s τ )(y)) (80) = Ē(y) ε 0 ∂ y θ s σ (y) . (81)
Using ( 67), ( 69) and the cell problem boundary condition in the last equation, we obtain

1 ε 0 τ -1 + τ 0 = δ(y) ε 0 , (82)
and therefore, using ( 36)

1 ε 0 τ -1 + τ 0 = δ(x) , (83)
which is the wanted result. Verifying (33) can also be done following a similar development.

Solving ( 69) and (70) in practice is not completely trivial. Indeed, the solutions strongly depend on the size of the domain Y and on the boundary conditions, which is unusual for a homogenization problem. This happens because the Dirac sources have a low-wave number content that makes the solution have a non-compact support. This problem can be efficiently fixed by removing the low-

Y. CAPDEVILLE

frequency content on the source and adding it back a posteriori. This solution also makes the cell problem solution independent of the boundary conditions. It is described in Appendix B.

To summarize, to obtain the effective force term and the corrector (here for the couple-source), one needs to follow three steps:

(i) solve the cell problem (44) with the starting source term (67);

(ii) remove the low-wavenumber content of the obtained corrector with (71) to obtain the final corrector;

(iii) build the smooth source term with (76).

As mentioned above, there is a technical extra step which is explained in Appendix B. Examples of these different steps are shown in Fig. A2 and commented in Appendix B.

Combined homogenization of the source and mechanical properties

In many situations, it can be necessary to homogenize both the mechanical properties and the source.

Two ways can be used to do so. The first way is to consider the homogenization problem of the source and the structure as a single problem. The second way is to proceed as two successive homogenization problems: first homogenizing the mechanical properties and then homogenizing the source. The first option is mathematically non-trivial. Furthermore, it can be interesting to have the option to use different values for ε 0 for the elastic properties and for the source, which cannot be done in that case.

Thus, we rely on the successive homogenization approach.

The idea is simple: we first homogenize the elastic properties following the procedure described in Sec. 2.2 to obtain the order 0 solution using ε 0 = ε e . It is the solution of the effective wave equation (26-27) based on the effective elastic parameters ρ * ,εe , E * ,εe . At this point the source is still a point source, but modified according to (31). We then apply the source homogenization procedure described in Sec. 2.3, but to u 0,e and to the effective properties obtained from the elastic homogenization step.

For this second step, we use ε 0 = ε s where ε s can be different from ε e . Limiting our discussion to the couple-source case, using (23), we have, in the first step:

u(x, t) = u 0,e (x, t) + ε e χ εe (x, y)∂ x u 0,e (x) + O(ε 2 e ) . (84)
Following (57), we then apply the source homogenization to u 0,e to obtain u 0,e (x, t) = u 0 (x) + θ τ (y)g(t) (85)

∂ x u 0,e (x, t) = 1 ε s ∂ y θ τ (y)g(t) + ∂ x u 0 (x, t) . (86)
In these last equations, u 0 is the solution obtained using the effective wave equation and the effective Effective seismic sources 19 source, itself obtained in the effective medium. Combining ( 84) and (85-86) together, assuming ε s /ε e behaves like 1, we obtain, to the leading order:

u(x, t) = u 0 (x, t) + θ τ (y)g(t) + ∂ y θ τ (y)χ ε 0 (x, y)g(t) + O(ε 0 ) , (87)
∂ x u(x, t) = ε 0 -1 ∂ y θ τ (y)g(t) + ∂ x u 0 (x, t) (1 + ∂ y χ ε 0 (y)) + O(ε 0 ) (88)
Here, ε 0 is either ε s or ε e : it does not matter because we have assumed their ratio behaves like 1. This, however, is not a necessity as we shall see another choice in 2-D.

So, in practice, to apply successively the two homogenizations, the following steps are required:

(i) compute the effective properties and the corrector χ ε 0 as described in Sec. 2.2;

(ii) compute the effective source and the corrector θ τ or θ f as described in Sec. 2.3 using the effective properties and not the original elastic properties. The point source needs to be corrected according to (31) ;

(iii) solve the effective wave equation using the effective properties obtained in (i) and using the effective source obtained in (ii) to obtain the leading order displacement u 0 ;

(iv) apply the correctors following (87-88).

Implementation of the effective source in SEM

Once the effective source term is known, implementing it in SEM is trivial. We typically have to compute source terms, which are, for example in the couple-source case,

(w, f ) = - Ω w(x)∂ x τ 0 (x) dx , ( 89
)
where w is a test function. The last integral can be computed by parts or not, it has no real influence.

In the heterogeneous case, we have a practical extra term that needs to be taken into account (see Appendix B), but it does not add any further difficulty to the SEM. Because the source term is smooth, we know that the solution u 0 does not present any discontinuity, and therefore, we can be sure that the solution is correctly computed everywhere, including at the source origin. Once u 0 is computed, the discontinuous near-source field is recovered by adding the corrector term θ τ (x/ε 0 )g(t) to it. discontinuity. This implies that the reference solution is not accurate in the element containing the source. Nevertheless, at the scale of the plots, it is barely visible (a small glitch can still be seen on the displacement snapshot for x = 15 m). We compute the effective source with ε 0 = 0.5 and use it in SEM to obtain the effective solution u 0 . Comparing the reference solution with the order 0 effective solution, both for the displacement u 0 and the strain ∂ x u 0 , we can see that the reference solution and the effective solution match perfectly everywhere, except near the source, for distances below λ min from the source location. Near the source, the difference in two solutions is significant.

Nevertheless, once the corrector is added, the effective solution and the reference solution match very well everywhere, including in the near-source field, both for the displacement and the strain. The effective source τ 0 and corrector θ τ used in this example are displayed in Fig. A2, left panels.

Defining the error as

error(x r , ε 0 ) = T 0 u 0 (x r , t) -u ref (x r , t) 2 dt T 0 (u ref ) 2 (x r , t) dt , ( 90
)
where u ref is the reference solution, in Fig. 3, we measure the error as a function of ε 0 , for a receiver position x r in the far-source-field (we use here x r = 3200 m). Starting from large ε 0 , the error is first flat until ε 0 3 and it then decays towards zero quickly until it reaches the numerical error level of the reference solution. This decay of the error as ε 0 becomes smaller is much faster than the decay in ε 2 0 observed when homogenizing the elastic properties (Capdeville et al. 2010a). This observation is consistent with the fact that the only non-zero term in the homogenization asymptotic expansion is the one in order 0: the convergence of the error towards zero as ε 0 becomes smaller is theoretically instantaneous. In practice, because of various numerical errors, the convergence with ε 0 is not exactly instantaneous but still much faster than the ε 2 0 convergence classically observed.

Heterogeneous bar, homogenization of the source only

In Fig. 4, we perform the same test as in the previous section but we use the heterogeneous bar, the one already used in Sec. 3.1. In this example, only the source is homogenized and not the mechanical heterogeneities in the bar. This implies that for computing the effective solution, once the effective source term is obtained, we still need a fine mesh to honor the mechanical discontinuities. This may seem not very useful but we will see in the discussion that this may be necessary for some situations.

We compare the reference solution with the order 0 effective solution once again, both for the displacement u 0 and the strain ∂ x u 0 , obtained with ε 0 = 0.5. We can see that, a few meters away from the source location (about one λ min ), the reference solution and the effective solution match perfectly.

For the strain, the imprint of the heterogeneous model is clear and is well captured in the far-field by the effective solution. Similar to the homogeneous bar test, for distances to the source below one λ min , the difference between the reference solution is significant. Once the zero-order corrector is added, (green lines, left panels) and ∂ x u 0 (green lines, right panels) without the corrector and the effective solution u 0 + θ τ (red lines, left panels) and 1/ε 0 ∂ y θ τ + ∂ x u 0 (red lines, right panels) are shown. The residual difference (magnified ×10) between the reference solution and the effective solution plus corrector is potted in grey line.

u 0 + θ τ g residual×10 residual×10 ∂ x u 0 + 1 ε 0 ∂ y θ τ g ∂ x u 0
ε 0 = 0.5 has been used and only the source has been homogenized, not the mechanical properties (there is no effective media).

the effective solution and the reference solution match very well in the near-source field, both for the displacement and the strain. Interestingly, the fine-scale mechanical heterogeneity effects are also captured in the near-source field thanks to the corrector. The effective source and corrector used in this example are displayed in Fig. A2, right panels.

Heterogeneous bar, successive homogenization for source and mechanical properties

We finally test the successive homogenization in the same numerical experiment framework as in the previous section. This case is probably the most common in realistic situations where reducing the numerical cost is the objective. We first homogenize the heterogeneous bar to obtain the effec- tive mechanical properties ρ * ,ε 0 , E * ,ε 0 and the first order corrector χ ε 0 , using ε 0 = ε e = 0.5. The corresponding effective velocity α * ,ε 0 = E * ,ε 0 /ρ * ,ε 0 is shown in Fig. 5. We then compute the effective source and the source corrector θ τ according to (76), solving (69) but use E * ,ε 0 instead of the original E, once again using ε 0 = ε s = 0.5 (note that we could have used a different value here).

Finally, we solve the wave equation using (ρ * ,ε 0 , E * ,ε 0 ) and the effective source in SEM to obtain the effective solution u 0 . In Fig. 6a and6b, comparisons between snapshots of the reference solution and the effective solution for both the displacement and the strain are shown. It can be seen that the far-field displacement is already well captured. This implies that for most situations, where only the far-field displacement is needed, u 0 is good enough. Nevertheless, the strain is missing the observed small-scale variations and both strain and displacement are missing the near-source field. If we add the source corrector to the displacement (u 0 + θ τ g, Fig. 6c) and to the strain (∂ y θ τ g/ε 0 + ∂ x u 0 , Fig. 6d), we see that the near-source field is recovered. The small scale variations are still missing on the strain.

They are also missing on the displacement, but they have a small amplitude and are not visible at this scale. Finally, we apply the full combined correctors according to (87-88) in Fig. 6e and6f. The small-scale variations are then recovered.

2-D AND 3-D THEORETICAL DEVELOPMENTS

We now consider waves in a 2-D/3-D elastic body Ω with free boundary conditions on ∂Ω. The elastic wave equation is where ρ(x) is the density, c(x) the elastic tensor, u(x, t) the displacement vector and σ(x, t) the stress tensor. The source terms f (x, t) and τ (x, t) are point sources:

ρ∂ tt u -∇ • σ = f , (91)
σ = c : ∇u -τ , (92)
• single force:

f (x, t) = Fδ(x -x 0 )g(t) , (93)
• couple-force or stress-source:

τ (x, t) = Mδ(x -x 0 )g(t) , ( 94
)
where F is a vector and M the moment tensor.

Once again, the stress source term is often inputted in the dynamic equations ( 91) using the following force term: f (x, t) = -∇ • τ (x, t).

In Fig. 7, we use a simple homogeneous domain Ω of size 40 × 40 km 2 in which the mechanical properties are such V P = 1.8 km/s, V S = 1.0 km/s with ρ = 1000 kgm -2 . The source is a point moment tensor whose mechanism is shown in Fig. 7e. As g(t), we use a Ricker with a central frequency f 0 = 0.33 Hz (f max 1Hz), such that λ min = 1 km. The central time of the Ricker is t 0 = 4 s. Two receivers are used, one very close to the source (receiver 1, at a distance of 100 m from the source location) and one in the far-field (receiver 2). We use SEM to solve the wave equation with PML absorbing boundaries around Ω (Festa and Vilotte 2005). We use two meshes: one trivial regular mesh (mesh (a), see Fig. 7a) and one unstructured mesh with smaller elements near the source (mesh (b), see Fig. 7b). For mesh (a), the element edges are of length 1 km and the polynomial degree in each element in each tensorial direction is 5. For such a simple medium, this trivial regular mesh is usually enough to obtain accurate seismograms for receivers not located in the element containing the source, such as receiver 2. Nevertheless, for a receiver close to the source origin, such as receiver 1, we need a denser mesh near the source, such as mesh (b), making sure the receiver is not in the element containing the source. As can be seen in Figs. 7c and7d, the solutions obtained using the two meshes match very well away from the source location, but near the source, they differ significantly. Even though the solution in the dense mesh can be considered as the reference solution near the source, it is still not accurate for the element containing the source. For SEM point source case, 'near' the source approximately means in the element of the source but elements right next to it may also be affected (see Appendix C). Interestingly and as it is well known for SEM, the fact that for the coarse mesh (a), the solution is not accurate in the element containing the source does not prevent the solution to be accurate for all the other elements (see, Figs. 7c and 7d for receiver 2).

As expected and similar to the 1-D case, the wavefield near the source appears to be discontinuous and a cross-section through the source shows a jump in the displacement (Fig. 7c). Of course, the SEM formulation does not allow a real jump, but the dense mesh clearly shows the expected jump at the figure scale (a zoom in the element containing the source would show the field is continuous, but it is an artifact of the formulation which does not make field discontinuities possible).

This simple test leads to the same conclusion as for the 1-D case: obtaining an accurate solution near the source requires a dense mesh. This mesh is potentially difficult to design when mixed with mechanical discontinuities, especially in 3-D, and leads to a high numerical cost. Once again, as we will see in the following sections, homogenization offers a solution to these technical difficulties and an interesting insight into the physical and mathematical processes involved.

Two-scale homogenization of the mechanical properties

The two-scale homogenization of the mechanical properties in 2-D en 3-D can be found in Capdeville, Guillot, and Marigo (2010b) and [START_REF] Cupillard | Non-periodic homogenization of 3-d elastic media for the seismic wave equation[END_REF]. While the 2-D or 3-D homogenization has many similarities with the 1-D case, it also has important differences. In dimensions higher than one, the order 0 stress σ 0 depends on y and, more importantly, there is no analytical solution to the cell problem, preventing an explicit formula such as (29). We therefore always need to use numerical solutions to obtain the effective media and the correctors.

We directly jump to the generalization of Sec. 2.3: the homogenization of the source term in 2-D.

Here, we only solve the case of the double-couple source, but the case of a single force source can be solved similarly without any specific difficulty.

The equations to solve are:

ρ∂ tt u i -∇ x • σ i -∇ y • σ i+1 = 0 , ( 95
)
σ i = c : x u i + y u i+1 -τ i , ( 96
)
where we have assumed that it is possible to find an expansion

τ (x, t) = 1 ε 0 τ -1 (x, y)g(t) + τ 0 (x)g(t) , (97)
such that u i and σ i belong to V, the higher dimension generalization of V (see Capdeville et al. (2010b) for more details). We have once again assumed that τ 0 is independent of y which, as we shall show later, is possible to achieve by construction.

Resolution of the homogenized equations

• Eqs ( 95) for i = -2 and ( 96) for i = -1 give

∇ y • σ -1 = 0 , (98)
σ -1 = c : y (u 0 ) -τ -1 . ( 99
)
Combining the last two equations, we have

∇ y • c : y (u 0 ) = ∇ y • τ -1 . ( 100
)
We define the leading order corrector θ τ as

u 0 (x, y, t) = u 0 (x, t) + θ τ (y)g(t) . ( 101
)
From the last equation, we see that having θ τ ∈ V * is a necessary and sufficient condition to have 101) in ( 100) we find that θ τ is solution of the following cell problem on Y:

u 0 ∈ V. Using (
∇ y • c(x) : y (θ τ )(y) = ∇ y • τ -1 (x, y) , (102)
with periodic boundary conditions on Y. The cell problem ( 102) is a static elasticity loading problem, the same as that found in the classical elastic-property homogenization. Only the source term is different. It can only be solved numerically either with a classical finite element tool (such as the one used in Capdeville et al. 2010b) or with an iterative FFT scheme [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Capdeville | Fast fourier homogenization for elastic wave propagation in complex media[END_REF]. In the following, we use a finite element tool based on triangular elements.

As it is seen in Sec. 3.3.2, θ τ only depends on y by construction.

Using ( 101) and ( 99), we have

σ -1 = c : y (θ y ) -τ -1 , (103)
which, by construction (see Sec. 3.3.2,Eq. 116), leads to

σ -1 = 0 . ( 104
)
• Eqs. ( 95) for i = -1 and ( 96) for i = 0, using (104), give

∇ y • σ 0 = 0 , ( 105
)
σ 0 = c : ( x (u 0 ) + y (u 1 )) -τ 0 . ( 106
)
Taking the cell average of the last equation, using the fact that for any vector v, y (v) = 0, we obtain the order 0 effective constitutive relation:

σ 0 = c : x u 0 -τ 0 . ( 107
)
• Eq. ( 95) for i = 0

ρ∂ tt u 0 -∇ x • σ 0 -∇ y • σ 1 = 0 , (108)
which, taking the cell average leads to effective dynamic equation

ρ∂ tt u 0 -∇ x • σ 0 = 0 . ( 109
)
At this stage, the equation driving u 0 are known and are ( 107) and ( 109). We can access u 0 by adding the corrector in (101). Similarly to the 1-D case, because all the higher order coefficients of the expansion can be shown to be null, we also have u(x, t) = u 0 (x, x/ε 0 , t). Moreover, knowing that

x (θ τ ) = 0 (because θ τ only depends on y by construction) and that u 1 = 0, from (106) we have

σ 0 (x) = c(x) : x u 0 (x) -τ 0 (x) . ( 110
)
We conclude that σ 0 only depends upon x and therefore that σ 0 = σ 0 . This is different from the homogenization of the mechanical properties where σ 0 depends on y. In that case, a small-scale corrector is required to access σ 0 from σ 0 . This leads to a complication for the boundary conditions, typically for the free surface. A boundary layer appears below the free surface and Dirichlet to Neumann operator for the boundary condition is necessary to go beyond the leading order [START_REF] Capdeville | Shallow layer correction for spectral element like methods[END_REF][START_REF] Capdeville | A non-periodic two scale asymptotic method to take account of rough topographies for 2-d elastic wave propagation[END_REF]. Here, for the homogenization of the source, the fact that σ 0 only depends on x implies that we do not expect a boundary layer near the free surface and that dealing with boundary conditions are simple, which is an important simplification compared to the homogenization of the mechanical properties case.

The principle to homogenize the medium first and then the source in 2-D is the same as in 1-D (see Sec. 2.4). Once the medium is homogenized with an ε 0 = ε e to obtain the effective properties (ρ * ,εe , c * ,εe ), the original moment tensor M is modified according to

M * = T G εe : M , (118)
where G εe is the strain corrector (the equivalent of the 1-D '1 + ∂ y χ ε 0 ' term in ( 30)) and T the transpose operator. More details can be found in Capdeville et al. (2010b).

Next, the point source is homogenized with ε 0 = ε s , using the effective medium and the effective moment tensor M * . Then, the leading order solution u 0 solution is obtained solving the wave equation with SEM using the effective medium and the effective source. Finally, the correctors are added to obtain

u i = u 0 i + θ τ,i g + ε e ε s χ εe ijk ∂ y j θ τ,k g + ε e χ εe ijk ∂ x j u 0 k + O(ε e ) , (119)
∂ x i u j = G εe klij ε s -1 ∂ y k θ τ,l g + ∂ k u 0 l + O(ε e ) , (120)
where G εe is the same strain corrector as for ( 118) but taken at the receiver location. Note that ( 119) is in O(ε e ) and not in O(ε 2 e ) as it is often the case for the displacement. This is because only the leading order correction is accounted for the point source in (118), and therefore, the whole process is in O(ε e ).

Extended sources

Once the point source case is solved, it is not difficult to treat the extended source. Starting from (3), discretizing the fault system S with an appropriate quadrature, we have

f (x, t) = i∈S d m(x i , t) • ∇δ(x -x i ) dS i , (121)
where S d is the set of integration points and dS i is the integration weight associated with each integration point. Then, for each integration point in S d , assuming a separation of time and space variables for the moment tensor density (m(x i , t) = m(x i )g i (t)), we solve the cell problem (113) using the moment tensor m(x i )dS i in (111). Doing so, we obtain a set of source potentials ψ i τ and correctors θ i τ associated with each integration points x i . Finally, the effective stress source is

τ 0 (x, t) = i∈S d c(x) : x ψ i τ (x)g i (t) , (122)
and the corrector

θ τ (x) = i∈S d θ i τ (x/ε 0 )g i (t) . (123)
The τ 0 above is the effective source that can be used in the wave equation to obtain the effective wavefield u 0 that accounts for all the source complexity in the far-field. To obtain an accurate nearsource field, we just need to add the corrector to the effective solution, as it is done for the case of the point source.

Finally, let us mention that, if the time-space variable separation assumption of the moment tensor density is not met, the problem remains similar. The difference is that the cell problem has to be solved for each time step of the discretized time history of each m(x i , t).

2-D numerical tests

Homogeneous domain

Here, we use the same setting as for the homogeneous domain Ω presented in Sec. 3.1 and Fig. 7. We follow the procedure presented in Sec. 3.3.2 to obtain the corrector θ τ and the source potential ψ τ .

θ τ and f * = ∇ • c : x (ψ τ ) are shown in Fig. 8. ε 0 = 0.5 has been used. f * is smooth and therefore very simple to implement in SEM as an external source term, making it possible to obtain u 0 using the simple mesh (a) (see Fig. 7a). The corrector θ τ g can be added to u 0 a posteriori to obtain the corrected effective solution. The results are displayed in Fig. 9. Once again, the solution in the dense mesh is used as a reference solution, which is a good assumption except in the element containing the source. At the scale of the plot, this element has the size of the line thickness. One can see that u 0 does not present the typical SEM Lagrange polynomial oscillation in the element containing the source as observed in Fig. 7c for the classical point source implementation. This is because the solution u 0 is smooth, in contrast to u. As expected, u 0 has a very good match with the reference solution except in the vicinity of the source. Once the correction θ τ g is added, the match near the source is also very good. Depending on the element size of the mesh used to compute θ τ near the source, the effective solution can be even more accurate than the solution in the dense mesh at no extra numerical cost for the SEM modelling.

Heterogeneous domain, homogenization of the source only

We now present a numerical test in a heterogeneous domain. In this test, only the source will be homogenized, not the medium. A test where both the medium and the source are homogenized is presented in the next section. The heterogeneous model P wave velocity is presented in Fig. 10, left Table 1. Velocities and densities used in the 2-D heterogenous test defined by their color as shown in Fig. 10 panel. The velocities and densities used for this test are given in Tab. 1. Like in the previous test, two SEM meshes are used: one standard mesh (Fig. 10, middle panel) and one with smaller elements near the source locations (Fig. 10, right panel). In that heterogeneous case, the standard SEM mesh not trivial but unstructured because mechanical discontinuities need to be honored by element boundaries.

The point source has the same moment tensor as in the previous test and is located near the tip of the heterogeneous wedge. Compared to the homogeneous case, this leads to strong waveform distortions as can be seen on the energy density snapshot (Fig. 10, middle and right panels). We follow the procedure presented in Sec. 3.3.2 to obtain the corrector θ τ and the source potential ψ τ using ε 0 = 0.5. θ τ and f * = ∇ • c : x (ψ τ ) are shown in Fig. 11. Compared to the homogeneous case presented in the previous section, f * is not smooth. This is because c is discontinuous. Nevertheless, this is not an issue for a SEM implementation. Indeed, the mesh already honors the mechanical interfaces and this implies that the presence of the discontinuous c is not a difficulty. Moreover, ψ τ is smooth (and

x (ψ τ ) as well), therefore, computing the SEM source term with f * can be done easily and accurately after being integrated by parts. Compared to the homogeneous case, it can be seen that f * is rotated and distorted. Using the homogenization procedure (Capdeville et al. 2010b), we can show the interaction Effective seismic sources 33 In Fig. 12, we compare the SEM solution in the dense mesh with the effective solution without correction u 0 and the effective solution with correction u 0 + gθ. One can see that the near-source field is correctly captured by the the effective solution with correction.

Successive homogenization for source and mechanical properties

Finally, using the same 2-D setting as in the previous section, we test the successive homogenization procedure. We first homogenize the elastic properties to obtain the effective medium (ρ * ,εe , c * ,εe ).

In Fig. 13, some of its components are shown, where ε e = 0.5 has been used. The effective elastic tensor is anisotopic, and we choose to represent V * ,iso P , the P -wave velocity computed with the nearest isotropic elastic tensor projection c * ,iso of c * ,εe , and the 'total anisotropy' measured as the matrix distance between c * ,εe and c * ,iso (see [START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-d numerical illustrations[END_REF] for more details).

The source is also homogenized using another ε 0 = ε s = 1, using the effective medium ρ * ,εe , c * ,εe to solve the cell problem (113). The resulting effective force vector f * is shown in Fig. 14.

Finally, we solve the wave equation using the effective model ρ * ,εe , c * ,εe and the effective source to obtain u 0 . Both the medium and the source are smooth, thus a trivial regular SEM mesh can be used.

A sample of this mesh is shown in Fig. 13, left panel. This simple mesh is finer than what it would be in an homogeneous medium in order to account for the oscillations of the effective medium and effective source (for more details about the mesh element size as a function of ε 0 , see [START_REF] Capdeville | Fast fourier homogenization for elastic wave propagation in complex media[END_REF], Eq. 18). The different correctors are added to u 0 a posteriori following (119). The left panel Left: cross-section in uz for t = 5.5 s, for the heterogeneous test case, used in Sec. 3.6.2 for the dense mesh (reference, black line), the order 0 effective solution u0 z (red line) and the order 0 effective solution plus corrector u0 z + ġθ z (green line). The amplitudes have been normalized by the maximum amplitude of the u0 z cross-section. Right: vertical component traces for receivers 1 and 2 for the reference solution (black line), u 0 z (red line) and the order 0 effective solution plus corrector u 0 z + gθ z (green line). The amplitudes have been normalized by the maximum amplitude for receiver 1.

V * ,iso 15 shows a cross-section in uz for the reference solution, in u0 z + ġθ τ,z using ε e = 0.5 and in u0 z + ġθ z + χ zkl (∂ y j θ τ,k ġ + ∂ x j u0 k ) using ε e = 0.125. The comparison shows that the both the far-field and the near field are accurately recovered with ε e = 0.5 using only the near-source field corrector.

Nevertheless, to model fine details of the near-source field (e.g. near x = 23.9 km) accurately, using a smaller ε e (here ε e = 0.125) and the structure corrector (χ εe ) are necessary. The same observation is made for traces (Fig. 15, right panels): in the near-source field, for the receiver 1, only 100 m away from the source location x 0 , a small ε e (0.25) and the the 3 corrector terms of Eq. 119 are necessary to obtain an accurate result. In the far-field, for receiver 2, ε e = 0.5 (still with ε s = 1) and no corrector is enough to obtain an accurate solution.

In the end, this example shows that we can obtain an accurate solution using a simple regular mesh both in the near-source field and the far-source field and for both a complex source and a point source.

Knowing that a more complex source can be seen as a sum over many point sources, it is easy to see how such a source could be accurately handled on the same trivial regular mesh.

DISCUSSION AND CONCLUSION

We used the non-periodic two-scale homogenization method to homogenize classical seismic point sources. The result is that, for a given frequency band, and therefore a given minimum wavelength, a point source can be replaced by a smooth effective source term. In the far-field, the wavefield generated by the effective source is the same as the true wavefield, up to the desired accuracy. In the near-source field, a corrector needs to be added to the effective wavefield to be accurately recovered.

For a media with fine-scale heterogeneities, we have shown two possibilities. In the first one, we just homogenize the source while keeping the medium unchanged. In this case, one still needs to mesh the medium discontinuities to obtain an accurate solution. In the second solution, we homogenize both the medium and then the source in a successive manner. This option makes it possible to use a simple mesh for both the medium and the source. Although it may seem that the second possibility is the only useful one, in some situations, one might want to keep the fine-scale medium. For example, some medium can be simple to mesh but still have small-scale such as discontinuities. Layered media Effective seismic sources 37 cross-section in uz for t = 5.5 s in the heterogeneous case used in Sec. 3.6.2 for the dense mesh (reference, black line), the order 0 effective solution u0 z + ġθ τ,z using ε e = 0.5 (red line) and the order 0 effective solution plus successive corrector u0

z + ġθ z + χ zkl (∂ yj θ τ,k ġ + ∂ xj u0 k )
for ε e = 0.125 (green line). The amplitudes have been normalized by the maximum amplitude of the u0 z cross-section (not shown). Right: vertical component traces for receivers 1 and 2 for the reference solution (black line), the order 0 effective solution plus corrector u 0 z + gθ z for ε e = 0.5 (red line) and u 0 z + gθ z + χ zkl (∂ yj θ τ,k g + ∂ xj u 0 k ) for ε e = 0.25 (green line). For the effective source, ε s = 1 has been used for each plot. The amplitudes have been normalized by the maximum amplitude for receiver 1. such as PREM are a good example of such models. In such cases, one might want to avoid having to homogenize the whole model to input an effective source. In those situations, the possibility to homogenize the source while keeping the small-scale structures of the medium might be desirable.

Regarding the classical two-scale homogenization technique, the homogenization of the source term in the wave equation presented here has some interesting and unusual aspects. Among them, we note that this case has no periodic equivalent. Indeed, for the homogenization of elastic media, it is always possible to use periodic heterogeneities, for which classical periodic two-scale homogenization is valid, to guide our intuition and test our results. Here, it is not possible to pave the elastic domain with periodically finely spaced sources: no propagating wavefield would result from such a source. The problem is therefore intrinsically non-periodic, which makes it unusual. Another aspect is that, for a couple-source term, the leading order solution has a corrector and therefore depends on the fine-scale y. This is different from most homogenization problems where the leading order solution is usually independent of the small-scale y. Lastly, the homogenization series is very short: depending on the source type, only one (for a single force) or two terms (for a couple) of the series are non-zero. This makes the convergence with ε 0 very fast, which is once again different from many homogenization problems.

The point source homogenization gives an interesting perspective on the fact, that, for SEM, the solution is inaccurate in the source element does not prevent the solution from being accurate anywhere else. On the one hand, the effective source τ 0 is enough to model the far-field, leading to a smooth wavefield, including at the source origin. On the other hand, the corrector θ τ is discontinuous: it cannot be accurately represented on the SEM polynomial basis but it has a zero average and it does not propagate energy. It is just the high wavenumber part of the static response of the medium to the point source. SEM is not able to accurately compute θ τ in the element containing the source, but it does not affect the far-field results because θ τ doesn't propagate energy and weak-form methods explicitly ensure energy conservation.

The interest of our work depends on the context of its application. In a forward modelling context, for weak-form methods such as SEM, the interest is a lower numerical cost only if one is interested in the near-source field. Indeed, for the far-field, classical direct implementation of point sources perform very well and there is no need for an alternate solution. Nevertheless, if one needs to model the nearsource field, especially for complex-fault systems, homogenized sources present a clear interest and make this possible based on a mesh that can ignore the fault complexity. It can also useful if one wishes to take advantage of large very high degree SEM elements, for which the probability to have a receiver in the element of the source is significant [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF]. For strong-form numerical solvers such as FD, this work is interesting even for the far-field: it gives a rigorous solution to obtain a distributed force map for the source, including in the case of a complex heterogeneous medium around the source.

In an inverse problem context, this work has interesting perspectives that need to be studied in future work. Indeed, while inverting for a moment tensor is a simple linear problem, inverting for its position is a strongly non-linear problem. Moreover, if the point source approximation is not valid anymore, there is no simple way to go beyond it. With that perspective, inverting for a smooth source potential ψ τ is appealing because knowing an approximate location of the source, the associated inverse problem is linear for everything, including the source position or its complexity. However, this idea remains to be studied and checked for feasibility. and, taking the cell average of the last equation, we obtain the order 1 constitutive relation:

σ 0 = E∂ x u 0 . (A13)
Taking the average of ( 38) for i = 0, we obtain the order 0 dynamic equation:

ρ(x)∂ tt u 0 -∂ x σ 0 = f 0 , (A14)
At this stage, we have found the order 0 homogenized solution. Going further, we find that all the higher-order terms of the homogenization expansion are zero.

APPENDIX B: PRACTICAL MATTERS TO SOLVE THE SOURCE CELL PROBLEM

As mentioned in Sec. 2.3.3 and Sec. 3.3.2, solving the source cell problem ( 44) and ( 64) in 1-D and (102) in 2-D leads to some difficulties. Indeed, the solution to the cell equations is not localized near the source and extends to the boundaries of the domain, leading to difficulties in filtering and separating the scales (see Fig. A1, left panel). In particular, it can lead to an effective source that extends to the whole domain, which is not optimal once introduced in the wave equation solver. Here, we propose a practical solution that localizes the solution around the source location and eliminates its dependency on the boundary conditions.

B1 1-D case, couple-source

The main idea is to have a smooth source term which, when used in the cell problem, gives a solution that behaves similar to the true solution (with the true source) source far away from its origin x 0 and then to remove it from the true source term. We first define a smooth space wavelet w εr (x) = F εr (δ(xx 0 ))(x) , (B1)

where ε r can be different from ε 0 . We also define wεr (y) = ε 0 w εr (ε 0 y), the same wavelet defined on the y axis. We then replace the Diract source term in ( 65) and ( 67) by f -1,s (y) = δ(y)α f wεr (y) , (B2)

τ -1,s (y) = δ(y)α τ wεr (y) , (B3)

where α f and α τ are coefficients designed to make sure that θ f and θ τ , solutions of ( 69) and ( 70) respectively, quickly go to zero away from the source origin. These two coefficients can be determined numerically, but an educated guess also performs very well. It E(x) is smooth with respect to λ min ,

α f = α τ = 1 is sufficient. If E(x)
is not smooth, α f = 1 is still good enough, but for the couplesource, we can use

α τ = E * ,εr E (x 0 ) , (B4)
and obtain a very good result. E * ,εr is the effective coefficient defined in (29). One can see the relation between the last equation and the effective source term in (31). One example of the effect of w εr on θ s τ is given in A1: as it can be seen, it removes the large-scales of the solution, allowing θ s τ to have a more compact support and to gently tends toward zero away from the source location.

We then follow the procedure described in Sec. 2.3.3, for f 0 (x) = ∂ x E(x)∂ x ψ f (x) + α f w εr (x) , (B5) τ 0 (x) = E(x)∂ x ψ τ (x) + α τ w εr (x) .

(B6)

The choice of the ε r value is up to the user and has a little impact on the final results. Choosing an ε r smaller than ε 0 leads to an effective source narrower than what is wanted with ε 0 . Thus, it is a good idea to choose ε r somewhat larger than ε 0 . Interestingly, if E(x) is not smooth, the choice of ε r affects τ 0 (or f 0 ) and the correctors, but once combined into the final result u 0 , it does not affect the solution.

Two examples of the different steps of the procedure are given in Fig. A2 for the couple-source case: one in a homogeneous bar and one in a heterogeneous bar (the same one as in Fig. 1). For this example, we use ε 0 = 0.5 and ε r = 2 (ε r = 4ε 0 ). The corresponding wavelet w εr function is shown in Fig. A1. We use (B3) with α f = 1 in the homogeneous case and α τ computed following (B4) in the heterogeneous case. Even if the cell problems can be solved analytically in this simple 1-D case, we use a finite element method to solve (69). The starting corrector θ s τ (Fig. A2, top panel) is, as expected, a step function with missing low-wavenumbers. One can note that the finite element methods introduce some error (oscillations) in the element to which belongs the source position x 0 and that the step function is approximate. Nevertheless, the element size can be as small as needed and this is not a problem as it has no impact on the wave propagation. The final corrector θ τ (Fig. A2, middle panels) shows only small-scale variations. The order 0 source expansion term τ 0 (Fig. A2, bottom panels) has no more singularities. In the heterogeneous case (right panel), it has discontinuities because of the choice to keep the elastic model discontinuity in the modelling. It would be completely smooth if the model was also homogenized or smooth from the start.

B2 2-D case, double-couple source

For the higher dimension case, the recipe is similar. We only develop the double-couple source case.

We first define a smooth space wavelet w εr (x) = F εr (δ(xx 0 ))(x) . . Left: positions of the source (blue star) and the receivers (red diamonds) relative to the mesh element boundaries (black lines). Right: cross-sections in uz along the dash line (right plot) for the solution obtained using the regular mesh (red line) and using the dense mesh (gray line, see Fig. 7b). The element boundary positions along the cross-section are marked by a vertical dash-line and the position of the source by a vertical dotted-line. A3, computed using the regular mesh (red line, see Fig. 7a) and the dense mesh (black line, see Fig. 7b). The dense mesh solution can be considered as the reference solution.

Figure 1 .

 1 Figure 1. Solution of the wave equation (5-6) computed in a heterogeneous bar, for the couple (8) and single force (7) point sources located in x 0 = 15 m. The source time function g is a Ricker of maximum frequency 1.5 kHz corresponding to a minimum wavelength λ min of approximately 3.3 m. Panel a: wave velocity in the elastic bar around the source position x 0 . b: normalized displacement snapshot around the source position at t = 1.5 ms after the Ricker central time for the couple-source. c: corresponding strain. d: normalized displacement snapshot around the source position at t = 1.5 ms after the Ricker central time for the single force source. e: corresponding strain. The absolute value of the maximum amplitudes of each displacement snapshots have been normalized to 1.

Figure 2 .

 2 Figure2. Homogeneous bar tests. Displacement (left panels) and strain (rigth panels) snapshots for t = 2 ms (top panels) and t = 4 ms (bottom panels) for the reference solution (black lines), the effective solution u 0 (green lines, left panels) and ∂ x u 0 (green lines, right panels) without the corrector and the effective solution u 0 + θ τ (red lines, left panels) and 1/ε 0 ∂ y θ τ + ∂ x u 0 (red lines, right panels). The residual difference (magnified ×10) between the reference solution and the effective solution plus corrector is potted in grey line. ε 0 = 0.5 has been used. The displacement amplitudes are normalized with respect to the maximum displacement amplitude for t = 4 ms.

Figure 3 .

 3 Figure3. Log-log plot of the error, defined in (90), for x r = 3200 m, as a function of ε 0 (solid line) for the homogeneous 1-D case defined in Sec. 2.5.1. y = ε 2 0 is also plotted for reference (dash line).

Figure 4 .

 4 Figure 4. Heterogeneous bar test. The displacement (left panels) and strain (rigth panels) snapshots for t = 2 ms (top panels) and t = 4 ms (bottom panels) for the reference solution (black lines), the effective solution u 0

Figure 5 .

 5 Figure 5. Original velocity α(x) (black line) and the corresponding effective velocity α * ,ε0 (x) (red line) for a sample of the bar used for the test in Sec. 2.5.3. ε 0 = 0.5 for a λ min 3.3 m has been used.

Figure 6 .

 6 Figure 6. Displacement (left column panels) and strain (rigth column panels) snapshots for t = 4 ms obtained for a couple-source and in the heterogeneous bar (see Fig. 3.1) for the reference solution (black line) and the successive homogenization schemes (red lines). a: u 0 , b:∂ x u 0 , c: u 0 + θ τ g, d: ∂ y θ τ g/ε 0 + ∂ x u 0 , e: u 0 + θ τ g + θ τ χg, f: (∂ y θ τ g/ε 0 + ∂ x u 0 )(1 + ∂ y χ). For e and f, the difference, magnified by a factor 10, between the reference solution and the effective solution plus corrector is also plotted (grey line).

Figure 7 .

 7 Figure 7. SEM wavefield modelling for the seismic point source moment tensor (e) propagating in a homogeneous medium Ω using a regular mesh (a) and a mesh with smaller elements near the source (b). In (a) and (b), a snapshot of the energy density of the propagating wavefield for t = 10 s (background color), the mesh (gray lines), the source (blue star) and the receiver (red diamonds) positions are plotted. (c): cross-sections in the velocity wavefield vertical component uz for t = 5.5 s along the dashed line shown in (a) and (b) for the regular mesh (black line) and the dense mesh (red line) modelling. (d): vertical velocity trace for receivers 1 and 2, for the regular mesh (black line) and the dense mesh (red line) modeling. The amplitudes have been normalized by the maximum amplitude for receiver 1.

Figure 8 .

 8 Figure 8. (a) and (b): effective source force f * = ∇ • c : x (ψ τ ) in the homogeneous case used in Sec. 3.6.1. (c) and (d): x and z components of the corrector θ τ , respectively. Each quantity amplitude has been normalized to 1.

Figure 9 .

 9 Figure9. Left: cross-section in uz for t = 5.5 s in the homogeneous case used in Sec. 3.6.1 for the dense mesh (reference, red line), the order 0 effective solution u0 z (black line) and the order 0 effective solution plus corrector u0z + ġθ z (green line). The amplitudes have been normalized by the maximum amplitude of the u0 z cross-section. Right: vertical component displacement traces for receivers 1 and 2 for u 0 z (black line) and the order 0 effective solution plus corrector u 0 z + gθ z (green line). The amplitudes have been normalized by the maximum amplitude for receiver 1.

Figure 10 .

 10 Figure 10. Left: domain Ω for the 2-D heterogeneous tests. The background colors represent V P and the dash line box the zoomed area used for the two right panels. Middle: zoom on the standard mesh used to solve the wave equation, designed to match the mechanical discontinuities. The background color represents an energydensity snapshot for the source located at the blue star and t = 5.5 s. Note that the seismic phase seen for this time is mostly the main S wave. The receiver positions are indicated with red diamonds. Rigth: same as the middle plot, but for the mesh with fine elements near the source (the "dense" mesh)

Figure 11 .

 11 Figure 11. (a) and (b): effective source force f * = ∇ • c : x (ψ τ ), for the heterogeneous test case, used in Sec. 3.6.2, with an amplitude normalized to 1. The green star show the original point source position. (c) true moment tensor. (d) apparent moment tensor due to the interaction of the true moment tensor and the heterogeneous elastic properties.

  Figure12. Left: cross-section in uz for t = 5.5 s, for the heterogeneous test case, used in Sec. 3.6.2 for the

Figure 13 .

 13 Figure 13. Left: nearest isotropic P-wave velocity V * ,isoP

Figure 14 .

 14 Figure 14. (a) and (b): effective source force f * = ∇ • c * ,εe : x (ψ τ ), for the heterogeneous test case, computed with the successive homogenization (see Sec. 3.6.3). ε e = 0.5 and ε s = 1 have been used. Each quantity amplitude has been normalized to 1.

Figure 15 .

 15 Figure 15. Successive homogenization results, where both the source and the medium are homogenized. Left:

Figure A1 .

 A1 Figure A1. Left panel: corrector θ s τ computed with using (69) with w εr = 0 (red line) and with w εr = 0 using (B3) (black line), both computed in a constant elastic model for a couple-source located in x 0 = 15 m. Rigth panel: w εr function used to obtain the black line corrector in the left panel.

Figure A2 .

 A2 Figure A2. Different steps used to build the effective source and corrector in two cases, in an homogeneous elastic bar (left column) and in the same heterogeneous bar as the one used in Sec. 3.1 (right column). The source is a couple-source with same properties as for Sec. 3.1. In the first row are shown the starting correctors θ sτ obtained solving (70) using (B3) with ε r = 4ε 0 (the corresponding w εr function is shown in Fig.A1), right panel). The middle plot row shows the final corrector θ s , obtained in (71) with ε 0 = 0.5. In the bottom plot raw are show the final effective source τ 0 , obtained in (B6). The four top graph amplitudes are normalized with respect to the maximum amplitude of the homogeneous case θ s τ .

  Figure A3. Left: positions of the source (blue star) and the receivers (red diamonds) relative to the mesh element

Figure A4 .

 A4 Figure A4. Vertical displacement u z for the four receiver positions shown in Fig. A3, computed using the

ACKNOWLEDGMENTS

The author would like to thank an anonymous reviewer, Gaetano Festa and the editor Carl Taper for their useful and constructive comments. I also thank Paul Cupillard, Zoé Renat, and Sneha Singh for fruitful discussions and practical suggestions for improving the manuscript. Numerical computations were performed using the "Centre de Calcul Intensif de Pays de la Loire (CCIPL)" resources. This Effective seismic sources 39 work was funded by ANR (Agence Nationale pour la Recherche) ANR-16-CE31-0022-01 "HIWAI" project. The programs associated to this work are available upon request.

Data availability

No new data were generated or analysed in support of this research.

We need to build the source expansion such that the displacement and stress expansion coefficients u i , σ i belong to V. For that, the following sufficient conditions must be met:

(i) θ τ must be in V * (from Eqs. ( 101) and (A10));

(ii) the expansion (97) must be satisfied.

We follow the same procedure as for the 1-D case: we first assume that all spatial variations are smallscale and then use the filter F to separate the scales on the correctors:

We then solve the y scale version of the cell problem (102) on Y with periodic boundary conditions, assuming the variation of the elastic tensor are also on y, and obtain the starting corrector θ s τ :

where c(y) = c(ε 0 y). We then separate the scales to obtain the final correctors

At this stage, θ τ belongs to V * by construction, and satisfy condition (i). Moreover, it only depends on y, also by construction.

We define the order 0 source potential from the smooth part of the starting corrector,

The order -1 source term is then built as

The last equation also justifies (104). Finally, the large-scale variations of the correctors is used to define the order 0 source terms:

We finally could check that the condition (ii) is satisfied following the same procedure as for the 1-D case.

Similar to the 1-D case, solving (113) in practice is not completely trivial for the same reason. The same solution, but generalized to the higher dimensions, can be used and is described in Appendix B.

APPENDIX A: RESOLUTION FOR THE SINGLE FORCE 1-D CASE

We assume in this appendix that

• Eqs (38) for i = -2 and (39) for i = -1 give

which imply

Knowing that E is positive, the last equation implies that ∂ y u 0 = 0 (see Capdeville et al. 2010a) and therefore,

This result is, unlike for the couple-source case, similar to the homogenization of mechanical parameters: the order 0 solution does not depend on the small-scale variable y.

• Eqs (38) for i = -1 and (39) for i = 0 give

Combining the two last equations, we have

We introduce the order 1 corrector

In the 1-D case, we can show that u 1 = 0 (Capdeville et al. 2010a). In the end,

We also define wεr (y) = ε 0 w εr (ε 0 y), the same wavelet defined on the y domain. We then replace the Dirac source term in (111) by

where α τ is a tensor coefficient designed to make sure that θ τ components quickly go to zero away from the source origin. This tensor can be determined numerically, but an educated guess also performs

where G εr is the strain corrector obtained for ε 0 = ε r and T the transpose operator.

We then follow the procedure described in Sec. 2.3.3,modifying (117) for