{'fraction_non_alphanumeric': 0.037402569764039266, 'fraction_numerical': 0.012476390625415657, 'mean_word_length': 3.952180213410618, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
True amplitude one-way propagation in heterogeneous media

This paper deals with the numerical analysis of two one-way systems derived from the general complete modeling proposed by M.V. De Hoop. The main goal of this work is to compare two different formulations in which a correcting term allows to improve the amplitude of the numerical solution. It comes out that even if the two systems are equivalent from a theoretical point of view, nothing of the kind is as far as the numerical simulation is concerned. Herein a numerical analysis is performed to show that as long as the propagation medium is smooth, both the models are equivalent but it is no more the case when the medium is associated to a quite strongly discontinuous velocity.


Introduction
The numerical solution of seismic acoustic waves propagation in heterogeneous media is generally based on the solution of the second-order full wave equation. Then the second-order wave equation can be solved completely by using a finite difference scheme and it is well-known that such an approach results in a high computational burden, especially in the case of three-dimensional problems. This surely explains that a lot of people prefer to solve an approximate problem which involves either a truncate expansion of the solution [10] or an approximate wave equation arising from the factoring of the exact one [18,13]. In the simplest case of a homogeneous medium, the solution can be obtained from the inversion of a system of one-way wave equations [1]. It is a very interesting way of solving the wave equation because it is based on the decomposition of the wave into a down-going part and an up-going one which reproduces the physical phenomenon very faithfully. During several years, people (see [7] and its references) tried to extend this approach by introducing correcting terms in the model to account for heterogeneities into the propagation medium. In 1996, M.V. De Hoop [8] derived a new formulation based on the micro-local analysis which allowed to derive a complete system of one-way wave equations coupling by exact correcting terms. The use of the theory of pseudo-differential operators makes the derivation of the system easy. Nevertheless, its plain writing is complicated because it involves the composition of pseudo-differential operators which means that each term is defined from an asymptotic expansion. Thus in practice the complete system is approximated by truncating the asymptotic expansions. Such an approach may look as if it complicates the solution as compared to the now well-controlled solution of the full wave equation. But its formulation allows to unpack the multiples from the primary reflections which is of outstanding importance for the geological interpretations and allows to reduce the computation time. Paper [8] has been followed by numerous publications and among them, we refer to as [16] in which one can find a complete bibliography on the topic. As far as the numerical solution is concerned, it is associated to the inversion of an approximate system generally based on only keeping the main term in each asymptotic development. J. Le Rousseau was the first to obtain accurate snapshots (see [15] and [9]). However if using his numerical method for the computation of arrival times, one gets erroneous results on the amplitudes level. More recently, Zhang et al. [18] have proposed a corrected one-way wave equation which allows to compute the correct amplitude of the acoustic pressure. This new one-way wave equation is obtained from the factorization of the full wave problem. The equation is factorized by using a WKBJ solution in which the amplitude is taken into account as well as its phase. Their new formulation of the one-way system includes a correction term. In this work, we intend to show that the amplitude of the numerical solution can be corrected by adding a transmission term in the system, proposed in [15]. The heterogeneities of the medium are modeled as discontinuities of the velocity which is supposed to vary in all the directions. The correcting term can be included in the system by two ways and we show that one of them is optimal. In fact we show that the best improvement of the amplitude is obtained when including the transmission operator into the right-hand-side of the system. This might come upon the reader since the model which is the nearest of the one of Zhang and al. [18] is obtained by including the transmission operator into the one-way equations, i.e. into the left-hand-side of the system. But some numerical tests indicate that the two approaches are close in the case of smooth media. The paper is divided into 7 sections, plus this one and a concluding part. The next one deals with the initial model whose unknown is the acoustic pressure and its transformation into a first-order system by introducing the velocity as an unknown also. The third part is devoted to the reduced system whose derivation is based on selecting the depth variable as the leading direction and by plugging the other terms into the frequencies domain after using a Fourier-Laplace transform. The fourth part concerns the first-order approximation of the reduced system. By accounting for the complete coupling terms of order 0, we get two equivalent systems which are interesting to consider since their numerical solution can be obtained by two different ways. The fourth next parts deal with some numerical aspects which are essential for the method. We have chosen to neglect the description of the propagation because we intend to focus on the transmission operator. The numerical tests are developed in the 2D case but we mention that some 3D test have been performed in [16].
In the following, we use standard notations for the micro-local analysis of classical pseudo-differential operators and we refer to [17] for their definitions. We only precise that the symbol of an operator P is denoted by σ(P ) and its principal symbol is σ P (P ).

Initial model
The analysis of the waves propagation is an efficient tool for imaging the soil. Assume the region of interest Ω is located between the surface of the earth {z = 0} and a given depth {z = z max }. The phenomenon of propagation is supposed to be generated by a source located at the top of Ω. Then the discontinuities of the propagation velocity can be defined by computing the time arrivals which are recorded by a set of receivers located at a given depth. By assuming that Ω is surrounded by two homogeneous regions Ω sup and Ω inf (see Fig.1), the receivers do not record any wave propagating below or above Ω. Hence only Ω is under study.
The phenomenon is governed by the wave equation set in IR 3 : whose solution is the acoustic pressure p = p(t, x). The nonnegative variable t denotes the time while x = t (x, y, z) stands for the position vector. In the following, x ′ = t (x, y) designates the transverse variable. The propagation velocity v := v(x) is constant-valued, equal to c sup in Ω sup and c inf in Ω inf .
In Ω, v := c(x) varies in all the directions. The source S is a Ricker function modelling an impulsion, where δ xs denotes the Dirac distribution at x s = t (x s , y s , 0), α = πν * and t * = 1/ν * . The constant ν * is a given frequency. We use standard notations for the Differential Calculus, such as: ∂ t denotes the partial derivative in the time variable, if ∂ i stands for the partial derivative with respect to i = x, y, z, divv = ∂ x v x + ∂ y v y + ∂ z v z represents the divergence of a vector field v with components v x , v y , v z and the operator ∇ is the gradient defined as ∇ϕ = t (∂ x ϕ, ∂ y ϕ, ∂ z ϕ).
Model (1) can be derived from two laws of the mechanic of continuous media (see for instance [11]) which involve the velocity in the fluid u = u(t, x) = t (u x , u y , u z ). The first concerns the conservation of the displacement quantity and reads as: for i = x, y, z, where f = t (f x , f y , f z ) designates the density of the exterior forces and ρ represents the density of the medium supposed to be constant herein. The stress tensor σ describes the behavior laws of the fluid whose lines are denoted by σ i with i = x, y, z. We use the notation: In the particular case of a perfect fluid, the tensor σ can be written as : where I 3 denotes the identity matrix of order 3. The law (2) modifies then into : In the framework of seismic prospection, the exterior forces are often negligible. The source is an impulsion which is taken into account in the second law describing the conservation of the mass: where q = q(t, x) represents an injection rate per mass unity. By applying the divergence to Eq.(4) and by deriving Eq.(5) with respect to the time, one can eliminate the velocity u and the acoustic pressure p is solution to the secondorder wave equation with right-hand-side (rhs) given by S = ∂ t q. Thus, we have: The wave equation (1) can then be rewritten as a first-order hyperbolic system of the form: For the numerical simulations, the direction of the depth is selected as the one governing the sense of propagation. This approach is classical in case of homogeneous media where the physical parameters do not vary. By using the formalism of pseudo-differential operators, it can be generalized to inhomogeneous media. In this work, we are interested in system (6) to which we associate a reduced system which is the subject of the following section.

Reduced system
The initial system is given by (6). We recall that the initial data are supposed to be null. We propose to describe the propagation of waves along the direction Oz. This idea was formerly exploited in [6] next in [7,1] where the problem was to compute solutions propagating into plane and homogeneous layers. In such a way, the problem turned into the solution of a system with constant coefficients and the solution was written as a linear combination of eigenvectors of the problem. This approach was successful for the study of the propagation of plane waves into 1D homogeneous layers [6] and the theory developed next in [7] provides an extension to the 2D case. Nevertheless, in case of more complex media, the constitutive parameters ρ and c of the domain Ω vary in all the directions. This is why we need the formalism of pseudo-differential operators to generalize this method of modelling, just as it was formerly suggested by M.V. de Hoop in [8].
In the following, we limit our attention to the solution in Ω. We eliminate the time variable by applying a Laplace transform to the system and the related variable to t is denoted by ω. System (6) then transforms into a stationary system, which reads in Ω as: where ∇ ⊥ = t (∂ x , ∂ y ), u ′ = t ( u x , u y ) and div ⊥ u ′ = ∂ x u x + ∂ y u y . In the above equations, ϕ represents the Laplace transform of ϕ. The time variable being suppressed, the transverse unknowns u ′ are eliminated by plugging the first equation into the third. Then we get a system with a pair of unknowns of the form: with U = t ( p, v z ), D z = I 2 ∂ z and I 2 represents the 2 × 2 identity. The operator L is defined by: and the source F is given by F = t (0, q) . By factoring by iω, L reads as L = iωL ♯ where According to the theory developed by Hörmander [12], Operator L ♯ is a pseudodifferential operator of OPS 0 depending on the parameter ω, and if L ♯ denotes its symbol, we have the following representation: for any test-function ϕ, where s ′ = t (s x , s y ) ∈ IR 2 and k ′ = t (k x , k y ) is the dual variable of x ′ such as the symbol of ∇ ⊥ is given by ik ′ . Symbol of L ♯ is defined by: Each term of symbol L ♯ is an homogeneous function of k ′ with degree 0. Hence it belongs to a class of matrices whose terms are in S 0 . Then it admits an asymptotic development with respect to the parameter |k ′ | ω which is small at high frequency, that is when ω is large. The operator L ♯ can then be seen as a semi-classical operator. Moreover the symbol of L = iωL ♯ admits an asymptotic development according to the small parameter ω −1 . Hence we are in the framework of the ω-calculus developed by O. Lafitte in [14] and L can be considered as a ω-pseudo-differential. This property is interesting for the construction of high-order models and is the subject of a current work [5]. To solve (7) turns into describing the propagation of the field U along the depth z. The formalism of pseudo-differential operators is well-adapted to extend the classical way for solving differential systems with constant coefficients and based on the diagonalization of matrix L. In the case of pseudo-differential operators (or differential operators with variable coefficients ), Taylor [17] has developed a method of factoring strictly hyperbolic systems. This method allows to replace the initial system by a system of two one-way equations. One will describe the down-going propagation and the other will model the up-going displacement. Both of the equations are coupled by a off-diagonal matrix which are assimilated to reflection terms. The idea of Taylor has been developed in [8,9], later in [15] for acoustic waves, in the simplest case where ρ is constant. Notice also that it has been applied in [2] for the construction of radiation conditions for electromagnetism.

First-order Formulation
The first-order formulation reads as a transport equation whose derivation is based on the factoring of the symbol L ♯ . The eigenvalues of L ♯ are given by the symbols γ 0 and −γ 0 with In order to define the square-root involved in (8), the plane IR 2 is divided into three regions. The first is defined by .
As soon as k ′ ∈ H, matrix L ♯ 0 admits two eigenvalues which are single and real. This region corresponds to the one in which System (6) is strictly hyperbolic and then H is called the hyperbolic region. It is exactly the region in which propagation occurs and the sign of each eigenvalue indicates the sense of propagation. The eigenvalue γ 0 is associated to the downgoing part of the wave i.e. the part propagating in the direction of increasing z. In the same way the eigenvalue −γ 0 is related to the upgoing part of the wave. The second region is defined by If k ′ ∈ E, the eigenvalues of L ♯ 0 are purely imaginary. System (6) is then elliptic in this region which is called the elliptic region. It defines a subset of frequencies which are not linked to the propagation. At last the third region is defined in which (6) degenerates. Indeed the eigenvalue 0 is double and the problem is neither hyperbolic nor elliptic. Region G is called the glancing region and is related to grazing rays.
Then we have the following result: Let m ∈ Z and U be a solution to the reduced system. There exists at least an operator where Λ 0 ∈ OP S 0 is a unique diagonal operator and R 0 ∈ OP S 0 depends on P m .
Proof. Let M 0 be the matrix defined as We aim at constructing a matrix P such that M 0 = P −1 L ♯ P. If w = t (w + , w − ) denotes one of the eigenvectors, associated with the eigenvalue ±γ 0 , its coordinates satisfy the equation: ρw − = ±γ 0 w + which admits an infinite number of solutions. We propose to solve this equation in such a way that the coordinates of the eigenvectors are in the same class of symbols. Hence if S m , m ∈ Z denotes this class [17], we define the conjugated matrix P m = [w + , w − ] where w + and w − have been chosen in S m . Then P m defines an operator of OPS m denoted by P m . Matrix P m is invertible and its inverse P −1 m is the principal symbol of the inverse of P m denoted by Q −m ∈ OPS −m . As far as the symbols are concerned, we have the relation Then we define the operator Λ 0 via its symbol by the relation: According to the pseudo-differential theory [17], we deduce that there exists a regularizing operator R −1 ∈ OPS −1 such that Operator R −1 is uniquely defined by its symbol whose expression is known thanks to the composition rule of pseudo-differential operators [17]. Next let us consider U as a solution to (7). We set V = Q −m U . Then V satisfies: and by composing on the left by Q −m , we get : where ∂ z P m is the pseudo-differential operator with symbol ∂ z P m , the system modifies into But by definition of L, we have : which implies, by taking (10) into account, Proposition 2.1 is then proved by setting : Hence the model is defined from a transport equation which involves two one-way equations of order +1. In the rhs, R 0 takes the coupling terms into account and they are modelled by the off-diagonal terms of the rhs. Our aim is to account for the lateral variations of the velocity. Hence the density mass can be constant or variable in all the directions. In this work, we limit our attention to the case ρ is constant.
The coupling operator can be split into the sum of two operators, the first being defined by the diagonal part of R 0 and the second by the off-diagonal part. Introduce By considering this decomposition, one can construct a new model in which the transport equation is diagonal and involves the sum of two operators in OPS 1 and OPS 0 respectively. Then we get: There exists at least an operator P m ∈ OP S m with inverse Q −m ∈ OP S −m such that if U denotes a solution to the reduced system, V = Q −m U is solution to: The model in 4.2 consists of two one-way equations which are coupled by the operator R ad 0 . According to its definition, R 0 is not easy to define exactly because it is given by an asymptotic expansion expressing its symbol. It is obvious that we cannot use the exact symbol of R 0 and that we must be content with approximating the system by a truncated one. The truncation order must then be defined. Since the initial reduced system is defined from an operator L in OPS 1 with a rhs in OPS 0 , we propose to keep this property in the one-way model. Hence we address the question of putting the principal part R 0,P in place of R 0 . Operator R 0,P is defined as the one whose symbol is σ P (R 0 ). Then we get the following result:

Proposition 4.3. Under the same assumptions than in Proposition 4.1, an approximate one-way model is given by
where the symbol of R 0,P is given by: Proof. We aim at computing the principal symbol of R 0 . We use the definition: We introduce the following notation. Let A and B be two matrices whose terms are symbols. Let By applying the composition rule of pseudo-differential operators [17], we get: Moreover we have where M −2 is the symbol of an operator in OPS −2 . Then, since we can deduce that By plugging (15) and (14) into (13), we complete the proof of Proposition 4.3.

Setting of the numerical method
The numerical method is now based on the solution to the approximate one-way model: where Λ 0 is the diagonal operator in OPS 0 whose symbol is the diagonal matrix and its opposite −γ 0 are the eigenvalues of L ♯ . The auxiliary unknown V whose components are the down-going field V + and the up-going one V − is linked to U by the relation U = P 0 V where P 0 is the operator in OPS 0 whose symbol is the matrix P 0 . We choose P 0 such as . Since R d 0,P does not play a part in coupling the down-going and up-going fields, it can be seen as a transmission term. Introducing a parameter ǫ ∈ {0, 1}, we write Then if ǫ = 0, R d 0,P acts in the right-hand-side of the model (like in the model proposed by De Hoop), and if ǫ = 1, R d 0,P is included in the left-hand-side (like in the one-way equations proposed by Zhang et al.). Only R ad 0,P accounts for the coupling and it can be seen as the reflection operator. In the following, we use the notations: R d 0,P = T 0 I 2 , R ad 0,P = R 0 J 2 with the letters T to indicate the Transmission and R corresponds to the Reflection. We introduce also the operator: To compute the solution to (16) requires to invert D z + M ǫ and the inverse operator is denoted by G ǫ , the so-called propagator. It governs the propagation along the depth, in the two senses. Like M ǫ , G ǫ is diagonal and Assume that G ǫ is given. Then (16) can be transformed into: Hence the solution is computed once (I 2 − K ǫ ) has been inverted. By construction, K ǫ ∈ OPS −1 because K ǫ arises from the composition of G ǫ ∈ OPS −1 and of ((1 − ǫ)T 0 I 2 + R 0 J 2 ) ∈ OPS 0 . Hence we have, according to [17]:

INRIA
This representation corresponds to a Neumann series for the inverse of I 2 − K ǫ , under the assumption K ǫ < 1.
Then we use (20) to expand V in the following form: This representation of V is called the Bremmer series and refers to as precursory Bremmer's works [6] in the simple case of a 1D model. Each term V (j) has two components denoted by V − where + corresponds to the down-going part propagating in the sense of increasing z and − is associated to the up-going part. By definition of the source F = t (0, q), the auxiliary source S has a priori two components S + et S − and the first term in the series is given by: Then the second term V (1) is obtained by the formulas: where K ǫ lm , 1 ≤ l, m ≤ 2 represent the terms of K ǫ . However the numerical solution is computed into a region limited by z = 0 and z = z max and which is surrounded by two homogeneous regions. By definition V − is possible and the case of a source located into Ω can be considered also. If K ǫ is represented from its principal symbol: we have σ P (K ǫ ) = 0 in Ω sup . Then since K ǫ is a pseudo-differential operator depending continuously on the parameter z, K ǫ l2 V (0) − = 0, 1 ≤ l ≤ 2. We then deduce that in that case, it is not useful to compute V while the next terms are linear combinations of the down-going and up-going fields: V (j) with j ≥ 2. If ǫ = 1, each term is simpler because K ǫ 11 = K ǫ 22 = 0, which implies that V (1) + = 0 and next, using the chain rule, we get:

RR n 0123456789
Hence each term is defined from either a down-going field or a up-going one which shows off the uncoupling of the computations. This property makes think out models based on paraxial equations and a comparison between oneway models and paraxial systems is a current work [4].
To summarize, the solution of (16) involves the following steps: 1. decomposition of the source: S = P −1 0 F 2. computation of the Bremmer terms from the propagator G ǫ , the reflection operator R 0 and the transmission term T 0 .
The center of the algorithm is given by the item 2. In order to illustrate its principle, we consider the case of a velocity model consisting of three homogeneous layers.

Bremmer series for a 2D stratified medium
The numerical solution is based on the expansion of the fields as a Bremmer series. We assume that Ω is defined as in Fig. 2. The first layer is homogeneous with thickness z 1 − δ 2 where δ is a small parameter and the corresponding constant velocity is denoted by c 1 . Next on a thin layer with thickness δ, the velocity continuously varies from c 1 to the constant value c 2 which is the propagation velocity in the layer with thickness (z 2 − z 1 ) − δ. Then on a thin layer with thickness δ, the velocity continuously varies from c 2 to the constant value c 3 .
The interfaces between each medium are flat and they are linked to variations of the velocity which are so important as δ is small. As long as the wave propagates into a homogeneous layer, no reflection phenomenon occurs. This is well-reproduced by the symbol R 0 since it is equal to 0 as soon as the z-derivative of the velocity vanishes, as in the case of a homogeneous medium. Then we have:

Series terms before any discretization
We propose to write the terms of the Bremmer series with item 0, 1 and 2. The next terms can be straightforwardly deduced from the term with item 2. Let us assume that if ϕ = t (ϕ + , ϕ − ) is a test-function depending on x ′ and z, we have: and where any G ǫ ± (h, z) is a pseudo-differential operator depending continuously on h and z. It has been proven in [16] that this assumption is available. Then the first term in the Bremmer series is obtained from the following arguments. First, just as was formerly noted, it is not useful to compute V + reads as: As far as the second term of the series is concerned, we have: where: Firstly consider the down-going term V + . Then, Property (22) implies that: and we deduce from (23) that: Next using (22) again, we get: and combining this relation with (23) yields: for any z into the homogeneous layer [z 1 + δ/2, z 2 − δ/2]. Now, if z belongs to [z 2 − δ/2 ≤ z ≤ z 2 + δ/2], we have: and in the homogeneous layer {z ≥ z 2 + δ/2} V (1) Moreover according to (22), we have: By applying the same reasoning than previously, we get: In the case where ǫ = 0, G 0 represents the propagation without accounting for the discontinuity of the medium. Hence V To construct V (2) we follow the same scheme than for V (1) : The component V + reads as follows: In the first layer, we have seen that σ(T 0 ) = σ(R 0 ) = 0. Thus we can deduce that + (z) = 0. Next developing the same ideas than for the construction of V (1) + , we get the following relations:

Formulas describing V
(2) + into the different layers are more intricate than the ones given V which can be split into each layer as: In the case where ǫ = 1, V − vanishes since V + vanishes too.

INRIA
The next terms V (j) − with j ≥ 3 read exactly in the same way than in the case j = 2 but by replacing the subscript (2) by (j) and the subscript (1) by (j − 1) into the formulas.
In the case where ǫ = 1, the method requires two times less operations than in the case where ǫ = 0. Hence one could already suppose that it will be more judicious to choose ǫ = 1. This is why we propose to investigate this question in the following section devoted to numerical experiments.

Numerical Bremmer terms
Let ∆z be the depth step. In practise, ∆z = δ where δ is the thickness of the layer in which the velocity varies continuously between two constant values. Each of the integrals involved into the Bremmer terms is approximated by the rectangle method and to choose ∆z = δ implies that the integrals in the variable z is replaced by the multiplication by ∆z. For instance, when the receivers are located at the surface z = 0, they are supposed to record V − (0) which is obtained step by step by computing the terms of the Bremmer series and we propose to truncate the series up to the third term. The approximate field is then given by: The terms are obtained by following always the same scheme. We focus our discussion on V − computed at a given depth z. Nevertheless the scheme depicted at Fig. 3 gives a general survey of the contribution of each term of the series. Assuming the propagator G ǫ ± (., .) satisfies the Chasles relation: the approximate term can be written as a function of the source: ǫ))T 0 G ǫ + (0, z 1 )S + Reflection at the second interface

INRIA
In Geophysics, the two first terms are the most important. But it is also of great outstanding to compute the multiples also because in the case where the velocity varies strongly, they can have a large amplitude and a simultaneous calculus of the two first events with the multiples can give rise to polluted results which are very difficult to interpret correctly. Hence it is very interesting to dispose of a numerical method which allows to compute the multiples and to uncouple them from the first events. Then the multiples can be removed from the primary fields which facilitates the interpretation of the results. To compute the multiples, it is necessary to combine the reflection operator at least three times. Precisely the m-th multiple is the result of 2m + 1 reflections. In order to explain the formula giving V −,app (0), we restrict ourselves to write V According to the Chasles (25) relation, we have: Thus we get the relation: for any z * ∈ [z, z max ], In practise, z * = z + ∆z. Hence if V The case of a stratified medium is very useful to show clearly why the cases ǫ = 0 and ǫ = 1 are different. Indeed, since the symbols do not depend on the variable x ′ , we can write the propagator in the wave number domain and the symbol t 0 of T 0 is non zero only at each interface. Hence if we consider a layer {z 1 < z < z 2 }, the propagators G ǫ ± are represented in the Fourier variables as [8,15,16]: g ǫ + = e −i(z2−z1)ωγ0 e ǫ∆zt0(z1) and g ǫ − = e −i(z2−z1)ωγ0 e −ǫ∆zt0(z1) By injecting these relations into the definition of V −,app (0), we get for the primary reflections: when ǫ = 0 V P −,app (0) = e −iz1ωγ0 (−∆zr 0 (z 1 ))e −iz1ωγ0 S + +e −iz1ωγ0 (Id − (∆zt 0 (z 1 )))e −i(z2−z1)ωγ0 (−∆zr 0 (z 2 ))e −i(z2−z1)ωγ0 (Id + (∆zt 0 (z 1 )))e −iz1ωγ0 S + and when ǫ = 1 We can then observe that the term (Id ∓ (∆zt 0 (z 1 ))) acting in the case ǫ = 0 is replaced by the exponential e ∓∆zt0(z1) in the case ǫ = 1. These two terms are of the same order when ∆zt 0 (z 1 ) is small, which is the case when the medium is smooth, i.e. when the vertical variations of the velocity are small.

Description of the software
We use the language Fortran 90. The code has been developed from the C one written by Jérôme Le Rousseau. It runs in the 2D and 3D cases. Its architecture has been modified as compared to the first version. Each routine is now independent of the other ones and can be replaced by anyone else. This is an interesting property for the analysis of the numerical method which can now be coupled with other methods like the ones involving finite difference schemes for instance. The computation of the Bremmer series has been optimized in the case where ǫ = 0 [16]. Hence, the computational cost is almost the same for ǫ = 0 and ǫ = 1.

Data
The source F is obtained from a FFT (Fast Fourier Transform) of the source q. One gets a table as a function of the frequency ω. The source q is a Ricker function which implies that its FFT decreases quickly to zero. This is why we define a window of calculus into the interval [0, ω max ]. Next one sets the number of terms in the Bremmer series by pointing at the number N of multiples which must be modelled. Then the code will compute V

Composition
For each value of m, we have computed the field recorded at z = 0 and this has been done for each value of the frequency from 0 to ω max . Then applying P 0 to this field, we get the unknown U . Next by applying an inverse FFT to its first component, we get the value of the acoustic pressure as a function of time. By depicting the acoustic pressure at the surface z = 0 we get a seismic section which represents the time value of the pressure as a function of the transverse variable x ′ . In all the pictures we state the variable x ′ is the abscissa and the ordinate gives the time oriented to the bottom. In this way we can represent the kinematic of the phenomenon. To depict the dynamic, we use a grey scale for the amplitude of the acoustic pressure at point (x ′ , t) and z = 0. It is also possible to represent snapshots (value of the wave field at (x ′ , z, t)). Then we use an auxiliary table T EM P := T EM P (x ′ , l, ω) in which we store the values of the field at (x ′ , z, ω) and by applying an inverse FFT, we get the field evaluated at (x ′ , z, t).

Setting of the velocity model
We consider a simple velocity model (See Fig.4) which consists of two layers. The topographic data of the model are the following: depth 5 km thickness 21.59 km The shot is set at the middle of the domain i.e. at x s = (10800, 0) and the receivers are located at the depth z = 3000m. By this way, we only account for the transmission effects. The transverse dimension of the model has been chosen such that the cardinal of the discrete set of points x j is matched for using FFTs. Each layer corresponds to a constant value of the velocity and their interface is located at z = 2500m. The mesh dimensions are collected in the following table:  In the following we refer to Velocity Model 1 (VM 1) when the velocity is equal to 1600 ms −1 in the first layer and 2400 ms −1 in the second one. VM 2 corresponds to the case where the velocity is equal to 1600 ms −1 in the first layer and 5000 ms −1 in the second one. The last model VM 3 is defined by a velocity equal to 2400 ms −1 in the first layer and 1600 ms −1 in the second one. VMs 1 and 3 are examples of velocities with a quite low contrast while VM 2 gives an example of high velocity contrast. On Fig.5, we have collected four seismograms. The first one depicts the reference solution obtained from the second-order full wave equation. On the right of this picture, we have set the solution of the one-way equations system as it was computed initially in [15]. One can observe that the kinematics is wellreproduced. Nevertheless the extrema of the grey scale shows that the amplitude of the one-way solution is erroneous. This motivates the two following pictures underneath. The left seismogram depicts the best results: both the kinematics and the amplitude of the wave field are correct. The right seismogram is not so bad even if the amplitude is more erroneous than in the previous case. Hence this first collection of numerical tests shows that to include the transmission term into the model really improves the amplitude of the wave field. Moreover if we limit our analysis to consider seismograms, we can conclude that the transmission term can be numerically handled as a term of the right-hand-side of the one-way system or as a proper term of the one-way equations. But the seismograms give a global view of the results and they are not precise enough to estimate the accuracy of the amplitude. This is why we have performed a series of numerical tests and we have chosen to represent them from the value of Q which is defined by: On Fig. 6, the lower curve depicts the value of Q for the one-way solution with ǫ = 0. Its values are very close to 1 which confirms what was observed on the corresponding seismogram. The top curve represents the values of Q for the one-way solution without the transmission term. This picture strengthens the previous conclusion since the minimum value of the error is 20%. Hence it is essential to include the transmission into the model to get accurate amplitudes. Now, being convinced that the transmission must be included, the theory shows that it can be equivalently introduced whether into the one-way equations or into the rhs of the system. The first way allows to write a system of one-way equations which are quite close to the equations considered by Zhang and al. [18] as an approximation of the full wave equation. The second one is the most natural in the formalism of Bremmer series and was introduced by Le Rousseau and De Hoop [15,8]. On Fig. 7, we have collected the values of Q for ǫ = 0 and ǫ = 1. When ǫ = 1, the results are correct in a neighbourhood of the shot and then they spoil quickly. Hence this numerical test indicates that for VM 1, the best approach  consists in including the transmission into the rhs. Thus the result for ǫ = 1 may deteriorate more and more as fast as the velocity contrast is high. That motivates the next test on Fig. 8 where we can observe that actually the results spoil more and more: the minimum value of the error is now 8% versus 1% for VM 1.
On Fig. 9, we have collected three curves obtained from the one-way solution with ǫ = 1. Each of them are distinguished from the value of the velocity contrast. The pattern shows that the lower the contrast is the larger the neigbourhood of the shot in which the results are satisfactory is. This figure seems to indicate that the model with ǫ = 1 is correct when the propagation medium is smooth.
On Fig. 10 we have depicted the results for VM 3. It is just to show that the same conclusion holds even if the upper velocity is larger than the lower one. Let us mention that each curve shows some instabilities (depicted as local maxima) which are due to the periodicities created by the FFTs.

conclusion
In this paper we consider the numerical analysis of two one-way systems derived from the general modeling of M. De Hoop [8]. Such a formulation is used to replace the full wave equation by a system of one-way wave equations whose computational burden is lower than the one associated to the finite difference solution of the wave equation. Moreover it permits to unpack and identify the multiples from the primary reflections. We have include a transmission term in the one-way model. Then numerical tests have been performed in the 2D case and they indicate that accounting for the transmission improves significantly the amplitude of the solution. The computational algorithm has been optimized in such a way that its complexity is now of the same order than the one of the one-way solver without transmission. Hence to add the correcting term does not penalize the computational method. Now the solution is based on a method of propagators whose numerical approximation accuracy is very sensitive to the position of this correcting term. To put it into the one-way system, which is the case ε = 1, seems to deteriorate the results. But this approach should be interesting since it corresponds to the same idea than the one proposed by Zhang et al. in [18] where time-arrivals are computed from the solution of a second-order wave equation obtained by factoring the full wave equation.
In this paper, we have performed a numerical test which shows that the accuracy of the method with ǫ = 1 is similar to the one of the method with ǫ = 0 when the velocity contrasts decrease. In the proposed methods, the transmission operator is approximated by a zero order pseudodifferential operator which is exact for stratified media. A higher order operator, that should account for media with lateral velocity variations, is the topic of a current research [3].