Parabolic equations are among the most popular numerical techniques in many fields of physics. This article considers extra-wide-angle parabolic equations, wide-angle parabolic equations, and narrow-angle parabolic equations (EWAPEs, WAPEs, and NAPEs, respectively) for sound propagation in moving inhomogeneous media with arbitrarily large variations in the sound speed and Mach number of the (subsonic) wind speed. Within their ranges of applicability, these parabolic equations exactly describe the phase of the sound waves and are, thus, termed the phase-preserving EWAPE, WAPE, and NAPE. Although variations in the sound speed and Mach number are often relatively small, omitting the second-order terms pertinent to these quantities can result in large cumulative phase errors for long propagation ranges. Therefore, the phase-preserving EWAPE, WAPE, and NAPE can be preferable in applications. Numerical implementation of the latter two equations can be performed with minimal modifications to existing codes and is computationally efficient. Numerical results demonstrate that the phase-preserving WAPE and NAPE provide more accurate results than the WAPE and NAPE based on the effective sound speed approximation.

The narrow-angle parabolic equation (NAPE) and wide-angle parabolic equation (WAPE) are among the most popular computational techniques in atmospheric acoustics,1–5 ocean acoustics,6–8 nonlinear acoustics,9,10 biomedical acoustics,11 and electromagnetic propagation.12 The NAPE and WAPE are well suited to relatively small computers, large domains, and high frequencies and can handle many phenomena such as atmospheric and ocean stratification and refraction; scattering by turbulence, internal waves, and other inhomogeneities; ground impedance and ocean bottom interactions; and propagation over slowly varying terrain and ocean bathymetry.

In atmospheric acoustics, the NAPE and WAPE should account for the effect of the wind velocity on the sound field. The corresponding equations have been formulated in a number of references, e.g., Ref. 2 and references therein. The WAPEs used in atmospheric acoustics are rather involved because they include the derivatives of the sound, density, and wind velocity, e.g., Refs. 13 and 14. These NAPEs and WAPEs are usually formulated assuming that the variations in the sound speed are small and the wind-speed Mach numbers are low.

A recent article15 suggests considering EWAPEs (extra-wide-angle parabolic equations or one-way equations), NAPEs, and WAPEs in the high-frequency approximation, in which the derivatives of the ambient quantities can be omitted. (Also, see Sec. 11.2 in Ref. 2.) The resulting PEs become much simpler than the previous formulations and still can be used in many applications. Specifically, Ref. 15 considers, in detail, the EWAPE given by Eq. (39) and the NAPE and WAPE derived from that equation. The WAPE from Ref. 15 has already been used in the literature, e.g., Refs. 16–19.

Reference 15 also suggests the other EWAPE, which is given by Eq. (B1), but considers the latter equation very briefly. One of the goals of the present article is to study this EWAPE in detail and determine when it is preferable compared to the EWAPE given by Eq. (39). To this end, starting with the EWAPE given by Eq. (B1) in Ref. 15, the NAPE and WAPE are derived. The newly formulated NAPE is remarkably simple and valid for arbitrary variations in the sound speed and Mach numbers. The new NAPE and the WAPE in the Padé (1,1) approximation can be implemented numerically with minimal modifications to existing codes and are computationally efficient.

In the literature, EWAPEs, NAPEs, and WAPEs are derived by omitting some terms in the Helmholtz-type equations. The ranges of applicability of the resulting equations are often difficult to assess. Hence, another goal of the present article is to study the ranges of applicability of these equations, which is examined by two approaches. In the first approach, starting with these equations, the dispersion relations (which are equivalent to the eikonal equations) are derived and compared to geometrical acoustics. When these dispersion relations agree with geometrical acoustics, the corresponding parabolic equations exactly describe the phase of a sound wave. The present article formulates the EWAPE, NAPE, and WAPE exactly describing the phase of a sound wave within the ranges of their applicability for arbitrary variations in the sound speed and Mach numbers. Therefore, these equations are termed the phase-preserving parabolic equations. On the other hand, many NAPEs and WAPEs from the literature preserve the phase of a sound wave only for low Mach numbers and/or small variations in the sound speed when the second-order terms can be omitted. Although the variations in the sound speed and Mach numbers are often relatively small, omitting the second-order terms might result in large cumulative phase errors for long ranges. For example, Ref. 15 argues that the second-order small terms cannot be omitted for the wind speed greater than or about 15 m/s; in this case, equations for arbitrary Mach numbers should be used. Thus, in many applications (such as sound propagation in the near-ground atmosphere and infrasound propagation in the upper atmosphere), it is desirable to use the phase-preserving EWAPE, NAPE, and WAPE.

The second approach for studying the ranges of applicability of parabolic equations is to compare them with the exact equation for sound propagation in stratified moving media. This article shows that the phase-preserving EWAPE, NAPE, and WAPE correctly describe sound propagation in stratified moving media within the ranges of their applicability.

This article is organized as follows. Section II briefly outlines results from Ref. 15 and presents new results such as the phase-preserving NAPE. In Sec. III, the ranges of applicability of EWAPEs, NAPEs, and WAPEs are studied by comparing them to the geometrical acoustics equations. In Sec. IV, these parabolic equations are compared to the exact equation in a stratified moving medium. Section V provides numerical algorithms for the phase-preserving NAPE and WAPE. Numerical results are presented in Sec. VI. The results are summarized in Sec. VII.

This section considers parabolic equations for sound propagation in moving media.

First, we briefly overview some results from Ref. 15 and provide additional analysis of these results. Reference 15 expresses the sound pressure p of a monochromatic sound wave in a moving medium in terms of the auxiliary function ϕ such that
(1)
Here, R = ( x , y , z ) is the vector of Cartesian coordinates, v ( R ) is the medium velocity, ϱ ( R ) and ϱ 0 are the density and its reference value, respectively, and ω is the angular frequency of the sound wave. The auxiliary function satisfies the convective Helmholtz equation,
(2)
where c ( R ) is the sound speed. Equations (1) and (2) are derived by applying the high-frequency approximation (for which the sound wavelength is smaller than the characteristic scale of medium inhomogeneities) to the linearized equations of fluid dynamics. This approach was originally suggested in Ref. 20. Additionally, Eq. (2) is derived assuming that v = 0. This equality is valid in a stratified moving medium. It is also approximately valid if the deviation from a stratified medium is relatively small or the derivatives of the medium velocity can be omitted.
Starting with Eq. (2), the following one-way equation is obtained [Eq. (38) in Ref. 15]:
(3)
Equation (3) describes sound propagation in the positive direction of the x axis and, following Refs. 15 and 21, is termed the EWAPE (or the parabolic equation for extra-wide angles). In Eq. (3), σ is a factor specified below and
(4)
Here, k0 and c0 are the reference wavenumber and sound speed, respectively; kc is the wavenumber in a motionless medium; vx and v = ( v y , v z ) are the velocity components in the direction of the x axis and the transverse plane, respectively; Mx and γx are the Mach number and Lorentz factor, respectively, pertinent to vx; = ( / y , / z ); and the operators m ̂ and τ ̂ act on the transverse coordinates r = ( y , z ). Furthermore, in Eq. (3),
(5)
where Δ = ( 2 / y 2 , 2 / z 2 ). Note that the factor σ appears in the notations given by Eq. (5) but not in Eq. (4).
Equation (3) is valid if the angle θ between the direction of sound propagation and the x axis is less than 90°. In deriving this equation, the derivatives of c and v are omitted; these derivatives are beyond the measurement accuracy in many applications. It is important to emphasize that Eq. (3) is valid for arbitrary subsonic Mach numbers, M = v / c, and arbitrary variations in the sound speed (provided that the derivatives of c and v can be omitted), which are characterized by the function ε:
(6)
The function ε is the deviation of the refractive index squared in a motionless medium from one. In contrast, many parabolic equations from the literature are obtained assuming that ε and/or M are small parameters, keeping terms of order ε and M, and omitting higher-order terms such as ε 2 and M2. However, omitting the second- and higher-order terms may result in significant phase errors for relatively long propagation ranges. For example, Ref. 15 argues that for the wind speed greater than 15 m/s ( M 0.05), the second- and higher-order terms in M should not be neglected such that parabolic equations valid for arbitrary Mach numbers should be used.
The EWAPE given by Eq. (3) does not actually depend on the factor σ. (This can be readily shown by moving σ premultiplying the square root into this square root.) However, when a NAPE or WAPE is obtained from Eq. (3) by approximating the square-root operator in this equation, it does depend on the choice of σ. We desire that the resulting NAPE and WAPE be valid for arbitrary variations in the sound speed and Mach numbers, respectively, i.e., arbitrary ε and M. First, consider the NAPE, which can be obtained from Eq. (3) by approximating the square-root operator with the first two terms of the Taylor series. This approximation can be used only if the terms containing the operators μ ̂ σ and χ ̂ σ in Eq. (3) are small compared to one. To estimate these operators, we use a zeroth-order approximation for the auxiliary function: ϕ exp [ i k 0 ( x cos θ + e r sin θ ) ], where e is a unit vector characterizing the propagation direction of the sound wave in the (y,z)-plane. For relatively small angles θ, for which the NAPE is applicable, we have | ϕ | k 0 θ ϕ and Δ ϕ k 0 2 θ 2 ϕ. Then, in Eq. (3), the operators m ̂ , μ ̂ σ, and χ ̂ σ can be estimated as
(7)
It follows from these relationships that for small θ, the terms containing the operator μ ̂ σ are small compared to one. However, the terms containing the operator χ ̂ σ are small only if σ 2 k 0 2 = k c 2 γ x 2 or
(8)
This is the desired value of σ for which the square-root operator in Eq. (3) can be accurately approximated with the first two terms of the Taylor series for arbitrary ε and M.

Next, consider the WAPE obtainable from Eq. (3) by approximating the square-root operator with the Padé series expansion. The Padé series converges for all values of the operators μ ̂ σ and χ ̂ σ and, hence, for all values of σ.22 However, if μ ̂ σ and χ ̂ σ are small compared to one, fewer terms in this series are needed for a good approximation of the square root than in the case of large μ ̂ σ and χ ̂ σ. (This result is illustrated in Table I of Ref. 23.) This, again, dictates the choice of σ given by Eq. (8), at least for small θ.

TABLE I.

Comparison between the NAPEs and WAPEs. The first two columns specify the parabolic equation and the function which satisfies this equation. The coefficients h n , m pertinent to these equations are in the last four columns.

Equation Function h 1 , 0 h 1 , 2 h 2 , 0 h 2 , 2
NAPE  ϕ ̂  n eff 1  1 2 n 
WAPE  ϕ ̂  1 4 n 2 γ x 2  n eff 1  1 2 n + n eff 1 4 n 2 γ x 2 
WAPEM  ϕ ̂  1 + n 2 1 4  1 4 γ x 2  ( n 2 1 ) γ x 2 2 ( 1 + n 2 1 4 ) τ  1 2 τ 4 γ x 2 
NAPEeff  p ̂  n eff 2 1 2  1 2 
WAPEeff  p ̂  1 + n eff 2 1 4  1 4  n eff 2 1 2  1 2 
Equation Function h 1 , 0 h 1 , 2 h 2 , 0 h 2 , 2
NAPE  ϕ ̂  n eff 1  1 2 n 
WAPE  ϕ ̂  1 4 n 2 γ x 2  n eff 1  1 2 n + n eff 1 4 n 2 γ x 2 
WAPEM  ϕ ̂  1 + n 2 1 4  1 4 γ x 2  ( n 2 1 ) γ x 2 2 ( 1 + n 2 1 4 ) τ  1 2 τ 4 γ x 2 
NAPEeff  p ̂  n eff 2 1 2  1 2 
WAPEeff  p ̂  1 + n eff 2 1 4  1 4  n eff 2 1 2  1 2 
Thus, the factor σ given by Eq. (8) is preferable for approximating the square-root operator in the EWAPE [Eq. (3)]. Substituting Eq. (8) into Eqs. (3) and (4), we obtain
(9)
(10)
Equation (9) is the desired EWAPE with the operators μ ̂ and χ ̂ given by Eq. (10). In these operators, the subscript σ is omitted because the factor σ is specified with Eq. (8). With σ given by Eq. (8) and for small θ, the operators μ ̂ and χ ̂ are estimated as
(11)
The terms containing these operators are small compared to one for small θ. This leads to the NAPE that is valid for arbitrary ε and M and requires fewer terms in the Padé series expansion for the WAPE.

Equations (9) and (10) coincide with Eqs. (B1) and (B2) in Ref. 15. However, Ref. 15 considers the EWAPE given by Eq. (9) only briefly in Appendix B. One of the main goals of this article is to study in detail this EWAPE and the NAPE and WAPE derived from this equation.

In the literature, the factor σ differs from that given by Eq. (8). In a motionless medium, σ is often set to one, e.g., Ref. 7. In this case, χ ̂ σ = ε and the Taylor series expansion of the square-root operator in Eq. (3) is valid only for small ε. In a moving medium, the bulk of Ref. 15 considers the EWAPE, NAPE, and WAPE pertinent to σ = γ x. With this choice of σ, the NAPE is valid only for small ε, whereas the WAPE requires more terms in the Padé series.

Thus, the EWAPE given by Eq. (9) might be preferable compared to other equations used in the literature because the terms containing the operators μ ̂ and χ ̂ are small compared to one for arbitrary ε and M (for small θ). This feature of Eq. (9) is achieved by the factor k c γ x 2 premultiplying the square-root operator, whereas other EWAPEs do not have this factor. The EWAPE given by Eq. (9) does not contain the reference wavenumber k0. It is worthwhile to point out that Eqs. (9) and (10) are simpler (i.e., have fewer terms) than the corresponding Eqs. (39) and (40) in Ref. 15, which pertain to σ = γ x. As expected, these equations coincide if c is constant.

In the literature, EWAPEs, NAPEs, and WAPEs are derived by omitting some terms in starting equations such as Eq. (2). The ranges of applicability of the resulting equations are often difficult to assess. Another goal of this article is to study the ranges of applicability of these equations in detail. To this end, two approaches are used. First, starting with the EWAPE, NAPE, and WAPE, the dispersion relations (which are equivalent to the eikonal equations) are derived and compared to geometrical acoustics. Section III shows that within the ranges of applicability of the EWAPEs given by Eqs. (3) and (9) ( θ < 90 °), the dispersion relations following from these equations coincide with geometrical acoustics. Therefore, Eqs. (3) and (9) are termed the phase-preserving EWAPEs.

In the second approach, EWAPEs, NAPEs, and WAPEs are compared with the exact equation for sound propagation in stratified moving media. Section IV shows that for θ < 90 °, the EWAPEs given by Eqs. (3) and (9) coincide with this exact equation.

Equations presented in this section are the starting point of the subsequent analysis.

To derive the NAPE from Eq. (9), the operator χ ̂ = 2 m ̂ + m ̂ 2 is substituted into this equation. Then, we take into account that in the square-root operator in the resulting equation, m ̂ M θ while μ ̂ k 0 2 θ 2 / ( k c 2 γ x 2 ); see Eqs. (7) and (11). Approximating the square-root operator with the Taylor series and keeping terms of order θ and θ 2 yields
(12)
Note that in this formula, terms of order m ̂ 2 cancel out. Substituting Eq. (12) into Eq. (9), we obtain
(13)
Here, c eff = c + v x is the effective sound speed, which appears in this equation as a convenient notation rather than an approximation. Equation (13), which is a new result, is the desired NAPE that is valid for the arbitrary variations in the sound speed and arbitrary Mach numbers, i.e., arbitrary ε and M. Equation (13) does not contain the reference wavenumber k0. After the auxiliary function ϕ is calculated with Eq. (13), the sound pressure is determined with Eq. (1). In a motionless medium, Eq. (13) coincides with Eq. (3.38) from Ref. 6 and Eq. (4.17) from Ref. 24 if in these equations, the derivatives of the sound speed are set to zero.

The NAPE [Eq. (13)] is valid for relatively small propagation angles θ. By studying the phase errors of a sound wave in motionless and moving media, Refs. 7 and 15 conclude that NAPEs are valid for the propagation angles θ 20 °. However, this conclusion depends on an acceptable phase error and the problem considered. For near-ground sound propagation, there is generally a consensus2,13 that a NAPE [including Eq. (13)] can be used if θ 20 °.

Section III shows that the dispersion relation derived from Eq. (13) coincides with geometrical acoustics within the ranges of applicability of Eq. (13); therefore, this equation is termed the phase-preserving NAPE. Section IV shows that within the ranges of its applicability, Eq. (13) coincides with the exact equation in a stratified moving medium. Given the wide range of applicability of the phase-preserving NAPE [Eq. (13)], it is remarkably simple. Numerical implementation of Eq. (13) is provided in Sec. V A.

Other NAPEs appearing in the literature are valid either for low Mach numbers and/or small variations in the sound speed. They preserve the phase of a sound wave only approximately and do not exactly coincide with the exact equation in a stratified moving medium. For example, the NAPE given by Eq. (46) in Ref. 15 is valid only for | ε | 1; that equation is also much more involved than the phase-preserving NAPE [Eq. (13)].

The NAPE in a moving medium is usually formulated for the sound pressure p rather than for the auxiliary function ϕ. It is often written in the form
(14)
[Equation (14) coincides with Eq. (2.111) in Ref. 2 if p is replaced with p ̂ exp ( i k 0 x ).] Equation (14) describes sound propagation only for small variations in the sound speed ε and low Mach numbers M. With the last term omitted, Eq. (14) can be obtained from the NAPE in a motionless medium using the effective sound speed approximation; that is, by replacing c with c eff.1 Because Eq. (14) is valid only for small ε and M the deviation of c eff from c0 is also relatively small, i.e., | c eff c 0 | / c 0 1. Therefore, keeping c eff 2 in Eq. (14) is inconsistent with its ranges of applicability. It can be shown that within these ranges, 1 + c 0 2 / c eff 2 in Eq. (14) can be replaced with 2 c 0 / c eff. In this case, the third terms in Eqs. (13) and (14) coincide.
WAPEs are usually derived by approximating the square-root operator in Eq. (9) with the Padé (m,m) series expansion such that
(15)
In this series, the coefficients a j , m and b j , m are given by25 
(16)
Substituting Eq. (15) into Eq. (9), we obtain
(17)
Equation (17) is the desired WAPE valid for the arbitrary variations in the sound speed and Mach numbers. This characteristic is attained because the Padé series converges for all values of μ ̂ and χ ̂. Note that Eq. (17) is formulated without the reference wavenumber k0. This equation is not explicitly presented in Ref. 15. In the Padé (1,1) approximation, only the first term in the series in Eq. (17) is retained. Numerical implementation of this WAPE is considered in Sec. V B. Keeping more terms in the Padé series enables larger propagation angles θ to be considered.

The Padé (1,1) approximation is a more accurate approximation of the square-root operator in Eq. (9) than the first two terms in the Taylor series. (For example, see Table I in Ref. 23.) Therefore, a WAPE in the Padé (1,1) approximation is valid for larger propagation angles than a NAPE. By studying the phase errors pertinent to WAPEs in the Padé (1,1) approximation, Refs. 7 and 15 conclude that these equations can be used if the propagation angles θ 35 °. As with NAPEs, this conclusion depends on an acceptable accuracy and the problem at hand.

Section III shows that the dispersion relation obtained from Eq. (17) coincides with geometrical acoustics within the ranges of applicability of Eq. (17). Therefore, Eq. (17) is termed the phase-preserving WAPE. In a stratified moving medium, Eq. (17) coincides with the exact equation, again, within the ranges of its applicability (See Sec. IV.).

As the Padé series converges for all values of the operators μ ̂ σ and χ ̂ σ in Eq. (3), the WAPE given by Eq. (44) in Ref. 15, which is obtained from Eq. (3) with σ = 1, is also phase preserving. However, as explained in Sec. II, more terms in the Padé series might be needed in Eq. (44) from Ref. 15 than in Eq. (17). Furthermore, Eq. (17) is simpler (has fewer terms) than Eq. (44) from Ref. 15. (See also Sec. V B.) Other WAPEs in moving media, such as those in Refs. 2, 13, and 14, are valid for small ε and M.

In this section, the ranges of applicability of EWAPEs, NAPEs, and WAPEs are studied by comparing them with the geometrical acoustics equations.

In moving media, the dispersion relation for geometrical acoustics [e.g., Eq. (3.35) in Ref. 2] can be written as
(18)
Here, k is the wave vector and k = | k |. For sound propagation in the positive direction of the x axis, it is convenient to write k = ( k x , k ), where k = ( k y , k z ). Squaring both sides of Eq. (18) and using notations from Eq. (4) yields
(19)
where k = | k |. Multiplying both sides of Eq. (19) by γ x 2, it can be recast in the equivalent form
(20)
Equation (20) is an exact consequence of the dispersion relation [Eq. (18)] and describes sound propagation in the positive and negative directions of the x axis. Equation (20) is recast as a 2 b 2 = ( a b ) ( a + b ) = 0, where a2 and b2 are the first and second terms in this equation, respectively. For sound propagation in the positive direction of the x axis, the latter equation reduces to a b = 0. Thus,
(21)
This is the desired dispersion relation for sound propagation in the positive direction of the x axis and is valid for θ < 90 °. In Eq. (21), k θ for small θ. Therefore, the square root has the form 1 + O ( θ ) + O ( θ 2 ) and can be effectively approximated using the Padé or Taylor series for arbitrary ε and M.
The wave vector k of a sound wave and its phase Φ are related by k = Φ. Substitution of this formula into Eq. (18) yields the eikonal equation for the phase Φ. For sound propagation in the direction of the x axis, it is convenient to write
(22)
Substituting kx and k from Eq. (22) into Eq. (21), we can obtain the eikonal equation for the phase Φ of a sound wave propagating in the positive direction of the x axis.
In this section, starting from Eqs. (1) and (3), we derive the dispersion relation and compare it with Eq. (21). The sound pressure p and auxiliary function ϕ are expressed in the form used in geometrical acoustics:
(23)
(24)
Here, Θ is the eikonal and k0 is the same reference wavenumber as in Eq. (4). In geometrical acoustics, 1 / k 0 plays the role of a small parameter and facilitates formulations of the dispersion relation and transport equation.
Equations (23) and (24) are substituted into Eq. (1). Keeping terms of order k 0 ( 0 ) and omitting higher-order terms, such as k 0 ( 1 ), yields
(25)
The geometrical acoustics equations can be written in forms that do not include the reference wavenumber k0. To this end, recall that the wave vector of a sound wave can be expressed in terms of the eikonal
(26)
Using the first of these relationships, Eq. (25) can be written without k0 such that
(27)
Next, Eq. (24) is substituted into the EWAPE [Eq. (3)]. To derive the dispersion relation (or the eikonal equation), we keep terms of order k0 and omit higher-order terms such as k 0 ( 0 ) , k 0 ( 1 ), etc. The first term in Eq. (3) is transformed as follows:
(28)
Here, following Eq. (26), k 0 Θ / x is replaced with kx. The last equality is obtained by taking into account that k x ϕ 0 is of order k0 while ϕ 0 / x is of order k 0 ( 0 ). Thus, when deriving the dispersion relation, the operator / x in Eq. (3) is replaced by ikx, i.e., / x i k x. The operator, m ̂ , in this equation can be assessed similarly:
(29)
Here, we used Eq. (26) for k and took into account that k ϕ 0 is of order k0 while ϕ 0 is of order k 0 ( 0 ). Thus, m ̂ m .
Equations (28) and (29) exemplify the rule that when one or more differential operators are applied to the ansatz e i k 0 Θ ϕ 0 F (where F is a function that may appear in transformations), they should act only on the exponential function e i k 0 Θ yielding terms proportional to k0, whereas the derivatives of F and ϕ 0 can be omitted to order k0. As a result of this rule, m ̂ n m n while the operators μ ̂ σ and χ ̂ σ appearing in Eq. (3) transform to
(30)
(31)
To complete the derivation, the square-root operator in Eq. (3) is expressed as the Taylor series and assessed as follows:
(32)
Here, we took into account that the operators μ ̂ σ and χ ̂ σ commute to order k0.
Using all these transformations of the operators in Eq. (3), the dispersion relation following from this equation takes the form
(33)
Taking into account the definitions of m and μ σ, it can be shown that this dispersion relation does not depend on σ and coincides with Eq. (21). Thus, the EWAPE given by Eq. (3) exactly describes the phase of a sound wave within the ranges of its applicability, i.e., θ < 90 °, for any σ.

The EWAPE given by Eq. (3) is derived from the convective Helmholtz equation [Eq. (2)]. The dispersion relation can also be formulated starting with Eq. (2). To this end, Eq. (24) is substituted into the convective Helmholtz equation. Similar to the analysis above, it can be shown that the operator transforms to i k and the operator Δ transforms to k 2. With these transformations, the dispersion relation following from Eq. (2) exactly coincides with that used in geometrical acoustics [Eq. (18)].

This section considers dispersion relations pertinent to NAPEs. To this end, Eq. (24) is substituted into Eq. (13). In the resulting equation, the operators v and Δ are transformed similarly to Eqs. (29) and (30), respectively. As a result, we obtain the following dispersion relation:
(34)
The same result can be derived from the dispersion relation in geometrical acoustics [Eq. (21)]. To this end, the square root in Eq. (21) should be approximated similarly to Eq. (12). Thus, within the ranges of its applicability, the NAPE [Eq. (13)] correctly describes the phase of a sound wave for arbitrary variations in the sound speed and Mach numbers.

Other NAPEs in moving media correctly describe the phase of a sound wave only approximately when M and/or ε are small. For example, starting with Eq. (14), the dispersion relation can be derived. It can be shown that it coincides with the dispersion relation given by Eq. (34) only for small ε and M and the NAPE given by Eq. (46) from Ref. 15 preserves the phase of a sound wave only for small ε.

Next, we consider the dispersion relation obtained from the WAPE given by Eq. (17), where the auxiliary function ϕ is specified by Eq. (24). The operators μ ̂ and χ ̂ are estimated with Eqs. (30) and (31), respectively, in which σ is given by Eq. (8) such that
(35)
(36)
Using these results, Eq. (17) becomes
(37)
The same equation can be obtained from the dispersion relation in geometrical acoustics [Eq. (21)]. To this end, the square root in Eq. (21) is expressed as the Padé series [Eq. (15)]. Therefore, in the Padé approximation, the phase-preserving WAPE [Eq. (17)] correctly describes the phase of a sound wave for arbitrary variations in the sound speed and Mach numbers, i.e., arbitrary ε and M.

The WAPE obtained from Eq. (3) in the Padé approximation also preserves the phase of a sound wave within the ranges of its applicability but might require more terms in the corresponding Padé series. (See Sec. II.) Other WAPEs appearing in the literature2,13,14 preserve the phase of a sound wave only for small ε and M.

In this section, the ranges of applicability of EWAPEs, NAPEs, and WAPEs are studied by comparing them to an exact result for sound propagation in a stratified moving medium.

As is commonly done in atmospheric acoustics, here, we assume a stratified moving medium. Thus, the ambient quantities depend only on the vertical coordinate z, i.e., c = c ( z ) , ϱ = ϱ ( z ), vz = 0, and v ˇ ( z ) = ( v x ( z ) , v y ( z ) ). Equations (2.63) and (2.64) in Ref. 2 exactly describe sound propagation in a stratified moving medium. For a monochromatic sound wave, by omitting the time dependence e i ω t, these equations can be recast as
(38)
(39)
Here, r ˇ = ( x , y ) , κ = ( κ x , κ y ), and the function ϕ ̂ ( κ , z ) is proportional to the two-dimensional (2D) spatial spectrum of the sound pressure. In this section, the check marks above r ˇ and v ˇ are used to distinguish these quantities from those in Secs. II and III. In Eq. (39), the derivatives of c, ϱ, and v ˇ are omitted, as consistent with the high-frequency approximation considered in this article.
Equation (39) can be recast in the following equivalent form:
(40)
Multiplying both sides of this equation by γ x 2, it can be written as
(41)
Because we omit the derivatives of c and v ˇ , Eq. (41) can be recast as a 2 b 2 = ( a b ) ( a + b ) = 0. For sound propagation in the positive direction of the x axis, the latter equation reduces to a b = 0. Thus,
(42)
Equation (42) exactly describes sound propagation in the positive direction of the x axis (i.e., θ < 90 °) provided that the derivatives of c and v ˇ can be omitted.
Next, we consider the EWAPE [Eq. (3)] in a stratified moving medium. The auxiliary function ϕ is expressed in the form
(43)
Substitution of Eq. (43) into Eq. (1) yields Eq. (38). Thus, the relationship between p and ϕ ̂ is the same as that in Ref. 2 and the formulations pertinent to the EWAPE given by Eq. (3). When Eq. (43) is substituted into Eq. (3), the operators are transformed as follows:
(44)
(45)
(46)
(47)
Here, the subscript s indicates that the corresponding quantity is pertinent to a stratified medium. When deriving these relationships, we also took into account that the operators m ̂ , μ ̂ σ, and χ ̂ σ commute if the derivatives of c and v ˇ are omitted. Because of this property, in Eq. (3),
(48)
With these results, Eq. (3) for ϕ transforms to the following equation for ϕ ̂:
(49)
This is the EWAPE for propagation in a stratified moving medium, which is derived from Eq. (3). Substituting with m , s, it can be shown that Eq. (49) does not depend on σ and coincides with Eq. (42). Thus, the EWAPE given by Eq. (3) exactly describes propagation in a stratified moving medium within the ranges of its applicability (for θ < 90 ° and when the derivatives of c and v ˇ can be omitted) for any σ.
Formulations for the NAPE proceed similarly to those for the EWAPE. To this end, Eq. (43) is substituted into Eq. (13). The operators / x and v ˇ are assessed with Eqs. (44) and (45). The result is
(50)
This is the NAPE for a stratified moving medium, which follows from Eq. (13). The same equation can be obtained from Eq. (42) by approximating the square root with Eq. (12).

Thus, within the ranges of its applicability, the phase-preserving NAPE [Eq. (13)] correctly describes sound propagation in a stratified moving medium with arbitrary variations in the sound speed ε and Mach number M. Other NAPEs [such as Eq. (14)] are valid only if M and/or ε are small.

Next, consider the WAPE [Eq. (17)] in a stratified moving medium, where ϕ is specified by Eq. (43). To derive the equation for ϕ ̂, the operators / x and v ˇ are replaced with Eqs. (44) and (45). With σ given by Eq. (8), Eqs. (46) and (47) for the operators μ ̂ σ and χ ̂ σ, respectively, take the form
(51)
Here, the subscript σ is omitted because the factor σ is specified with Eq. (8). Also, similar to Eq. (48), the following transformation holds:
(52)
With these results, the WAPE for ϕ ̂ takes the form
(53)
This is the WAPE in a stratified moving medium, derived from Eq. (17). The same equation can be obtained from Eq. (42) by approximating the square root with the Padé series [Eq. (15)]. Thus, within the ranges of its applicability, the phase-preserving WAPE [Eq. (17)] coincides with the exact equation for sound propagation in a stratified moving medium for arbitrary variations in the sound speed ε and Mach numbers M.

It can be shown that the WAPE given by Eq. (44) in Ref. 15 also coincides with the exact equation for sound propagation in a stratified medium within the ranges of its applicability. However, this equation might require more terms in the Padé series than in Eq. (53), as discussed in Sec. II. The WAPEs in Refs. 2, 13, and 14 coincide with the exact equation in a stratified moving medium only approximately for small ε and M.

In this section, we consider numerical implementation of the phase-preserving NAPE [Eq. (13)] and WAPE [Eq. (17)] for 2D sound propagation in the vertical (x,z) plane. As is typical for sound propagation modeling, we also assume that vz = 0.

First, consider the phase-preserving NAPE. For 2D sound propagation, the operators in Eq. (13) simplify such that v = 0 and Δ = 2 / z 2. In this equation, the auxiliary function ϕ is expressed in the form
(54)
Here, k0 is the reference wavenumber and ϕ ̂ is the complex amplitude of the auxiliary function. With these results, the phase-preserving NAPE [Eq. (13)] takes the form
(55)
where n is the refractive index in a motionless medium and n eff is the effective refractive index,
(56)
Note that n eff appears in Eq. (55) as a convenient notation rather than the result of the effective sound speed approximation.
For further analysis, it is convenient to express Eq. (55) in the equivalent form
(57)
where the coefficients h n , m are given by
(58)
The initial condition for the function ϕ ̂ at x = 0 is chosen as that adopted in the literature for a point source and NAPE,1,6
(59)
Here, ϕ 0 characterizes the amplitude of the point source and zs is the source height above the ground. In numerical examples to follow, ϕ 0 = 1 Pa. Numerical implementation of Eq. (57) is outlined in Sec. VI A of Ref. 15. It involves only minimal modifications to existing Crank-Nicholson NAPE algorithms and is just as computationally efficient.
After the auxiliary function ϕ ̂ is found as a solution of Eq. (57), the sound pressure is calculated as
(60)
This result follows from Eqs. (1) and (54). Numerical implementation of Eq. (60) is also considered in Sec. VI A of Ref. 15.
In the literature, the sound pressure in a moving medium, such as the atmosphere, is often calculated with the NAPE based on the effective sound speed approximation or the NAPEeff.1 This equation is given by Eq. (14) with the last term omitted. In the NAPEeff, the sound pressure is expressed in the form
(61)
where p ̂ is the complex amplitude of the sound pressure. For 2D sound propagation calculations using Eqs. (14) and (61), p ̂ also satisfies Eq. (57) but with different coefficients h n , m, namely,
(62)
The initial condition for the NAPEeff is formulated for p ̂ ( 0 , z ) and is given by the right-hand side of Eq. (59).
Next, consider the phase-preserving WAPE [Eq. (17)]. For 2D sound propagation in the (x,z)-plane, the operators χ ̂ and μ ̂ in this equation simplify to
(63)
In Eq. (17), ϕ is replaced with Eq. (54). With these transformations and in the Padé (1,1) approximation, the phase-preserving WAPE takes the form
(64)
It follows from Eq. (16) that in the Padé (1,1) approximation, a 1 , 1 = 1 / 2 and b 1 , 1 = 1 / 4. Multiplying both sides of Eq. (64) by 1 + b 1 , 1 μ ̂, this equation can be recast as Eq. (57), in which the coefficients h n , m are given by
(65)
Although written in a somewhat different form, Eq. (65) for the coefficients h n , m is consistent with Eq. (B5) in Ref. 15. The initial condition for the WAPE in the Padé (1,1) approximation is chosen as that adopted in the literature [see Eq. (G.75) in Ref. 1]:
(66)
In the literature, the sound pressure p in a moving medium is often calculated with the WAPE based on the effective sound speed approximation or the WAPEeff, e.g., see Ref. 1 and Eq. (56) in Ref. 15. In the Padé (1,1) approximation and substituting with Eq. (61), the WAPEeff can also be written as Eq. (57) with the following coefficients:
(67)
The initial condition for the WAPEeff coincides with the right-hand side of Eq. (66).

It follows from Secs. V A and V B that the phase-preserving NAPE and WAPE and the corresponding equations based on the effective sound speed approximation, NAPEeff and WAPEeff, reduce to the same equation [Eq. (57)] but with different coefficients h n , m. These coefficients are summarized in Table I. The rows “NAPE” and “WAPE” correspond to the phase-preserving NAPE and WAPE [Eqs. (13) and (17)], and the rows “NAPEeff” and “WAPEeff” correspond to the equations based on the effective sound speed approximation. The row “WAPEM” corresponds to Eq. (78) in Ref. 15, which is the 2D WAPE in the Padé (1,1) approximation and can be obtained from Eq. (3) in this article by setting σ = 1 and written in the form given by Eq. (57). As mentioned in the Introduction, Eq. (78) from Ref. 15 has already been used in the literature16–19 and is valid for arbitrary Mach numbers but small variations in the sound speed. In the WAPEM row, τ = ( n M x ) M x γ x 2.

Interestingly, although the parabolic equations in Table I describe similar sound propagation phenomena, the coefficients h n , m in these equations differ significantly. Also, these coefficients are more complex for the WAPEM than for the phase-preserving WAPE. This result elucidates the statements above that the EWAPE, NAPE, and WAPE, which are obtained from Eq. (3) by setting σ equal to Eq. (8), are simpler than those pertinent to σ = 1.

Another difference between these parabolic equations is that the phase-preserving NAPE and WAPE and the WAPEM are formulated for the complex amplitude of the auxiliary function ϕ ̂, while the NAPEeff and WAPEeff are formulated for the complex amplitude of the sound pressure p ̂.

In this section, numerical results obtained with the phase-preserving NAPE and WAPE and the NAPEeff and WAPEeff are compared to analytical and reference solutions. We will also briefly consider results obtained with the WAPEM.

In all numerical examples to follow, we consider a point source with the frequency f = ω / ( 2 π ) = 200 Hz located above a perfectly reflecting, flat surface. The reference sound speed is c 0 = 340 m/s. The computational domain is 500 m long and 125 m high. At the lower boundary, the vertical derivative of ϕ for NAPE and WAPE and the vertical derivative of p for NAPEeff and WAPEeff are zero due to a perfectly reflecting surface. For the upper boundary, a perfectly matched layer (PML) is used.26 The thickness of the PML is 30 grid points. It is implemented as described in Ref. 19. In all numerical examples, for x greater than about 100 m but smaller than 500 m, the propagation angles θ 35 °.

This section considers the point source located 50 m above the surface in a uniformly moving medium with Mach number M = 0.5 and constant sound speed c = c0. The relatively large Mach number is chosen to elucidate the difference between the results obtained with the phase-preserving NAPE and WAPE and the NAPEeff and WAPEeff. Note that for the considered case c = c0, the WAPEM coincides with the phase-preserving WAPE, i.e., the coefficients hnm in Table I for these parabolic equations are the same. Using these parabolic equations, the sound pressure level (SPL) and relative level are calculated. The latter is defined as follows
(68)
where p0 is the sound pressure produced by the point source in free space. In these calculations, the steps along the x and z axes are 0.05 m, which corresponds to λ / 34. Here, λ is the sound wavelength. These relatively small steps are needed to correctly predict the sound pressure at interference minima. In practical applications, the grid steps in the WAPE and NAPE can be set to those typically used in the literature, e.g., λ / 10.1 Therefore, the computational cost of these equations is about the same as for the NAPEeff and WAPEeff.

The SPL calculated with the parabolic equations is depicted in Fig. 1, which also presents the SPL corresponding to the exact analytical solution of the problem considered, presented in Sec. VI B of Ref. 15. The SPL has several minima and maxima due to interference between the direct wave from the source to receiver and that reflected from the surface. For x 100 m, the SPL calculated with the phase-preserving WAPE is indistinguishable from the exact analytical solution. The phase-preserving NAPE results are also close to the analytical solution but at larger propagation ranges. The results obtained with the WAPEeff and NAPEeff noticeably deviate from the analytical solution.

FIG. 1.

(Color online) SPL in a uniformly moving medium with M = 0.5 produced by a point source located 50 m above a perfectly reflecting, flat surface. The sound frequency is 200 Hz. The subplots correspond to the exact analytical solution and numerical results obtained with the parabolic equations as labeled.

FIG. 1.

(Color online) SPL in a uniformly moving medium with M = 0.5 produced by a point source located 50 m above a perfectly reflecting, flat surface. The sound frequency is 200 Hz. The subplots correspond to the exact analytical solution and numerical results obtained with the parabolic equations as labeled.

Close modal

These conclusions are reinforced by Fig. 2, in which the relative level Δ L obtained with the analytical solution and parabolic equations, is plotted versus the propagation range x for four heights z above the surface. For x 100 m, propagation angles θ can be relatively large such that the ranges of applicability of the NAPEs and WAPEs might not be fulfilled. Therefore, the predictions based on the parabolic equations can significantly deviate from the analytical solution. For x 100 m, the results obtained with the phase-preserving WAPE are very close to the analytical solution. The results obtained with the phase-preserving NAPE slightly deviate from the analytical solution in the range 100 m x 400 m because the NAPE is valid for smaller θ than the WAPE. However, for x 400 m (where θ becomes small), the NAPE is close to the analytical solution. On the other hand, Fig. 2 clearly shows that the results obtained with the WAPEeff and NAPEeff significantly deviate from the analytical solution. For the problem considered (but with different source height and frequency), inaccurate predictions of WAPEeff are also reported in Ref. 15.

FIG. 2.

(Color online) Relative level Δ L versus the propagation range x for four heights z above the surface. The geometry of the problem is the same as that for the SPL in Fig. 1. Different curves correspond to the analytical solution and parabolic equations.

FIG. 2.

(Color online) Relative level Δ L versus the propagation range x for four heights z above the surface. The geometry of the problem is the same as that for the SPL in Fig. 1. Different curves correspond to the analytical solution and parabolic equations.

Close modal
The difference between the results obtained with the parabolic equations and analytical solution can be quantified by the normalized errors,
(69)
where p PE is the sound pressure calculated with one of the parabolic equations from Table I, and p an is the sound pressure obtained with the exact analytical solution. (Note that the errors ϵ are normalized by p0 rather than by p an because p an can be very small at the interference minima; see Fig. 2.) The normalized errors ϵ versus the propagation range x are plotted in Fig. 3 for four heights z above the ground. These errors correspond to the results shown in Fig. 2. In Fig. 3, the x axis starts with x min = 100 m because at shorter ranges, the parabolic equations might not applicable. It follows from Fig. 3 that the normalized errors ϵ pertinent to the phase-preserving WAPE are almost zero. The normalized errors ϵ for the phase-preserving NAPE are small and tend to zero as the propagation range x increases. The normalized errors pertinent to the WAPEeff and NAPEeff are large and, generally, do not decrease with range.
FIG. 3.

(Color online) Normalized errors ϵ between the analytical solution and results obtained with the parabolic equations from Fig. 2.

FIG. 3.

(Color online) Normalized errors ϵ between the analytical solution and results obtained with the parabolic equations from Fig. 2.

Close modal
Next, we consider the cumulative normalized errors, which are obtained by integrating the normalized errors ϵ along the x axis:
(70)
where x max = 500 m. The cumulative errors pertinent to the normalized errors in Fig. 3 are shown in Table II. It follows from the table that the phase-preserving WAPE has the smallest cumulative errors, followed by those for the phase-preserving NAPE and the NAPEeff and WAPEeff.
TABLE II.

Cumulative errors ϵ ¯ pertinent to the normalized errors ϵ shown in Fig. 3.

Equation z = 1 m z = 2 m z = 5 m z = 10 m
NAPE  0.03  0.03  0.09  0.15 
WAPE  0.00  0.01  0.01  0.02 
NAPEeff  0.44  0.32  0.45  0.43 
WAPEeff  0.46  0.51  0.59  0.76 
Equation z = 1 m z = 2 m z = 5 m z = 10 m
NAPE  0.03  0.03  0.09  0.15 
WAPE  0.00  0.01  0.01  0.02 
NAPEeff  0.44  0.32  0.45  0.43 
WAPEeff  0.46  0.51  0.59  0.76 
This section considers sound propagation in a motionless medium. The source is located 10 m above the surface. The sound speed linearly increases with the height z such that
(71)
where α = 1 s−1. As in Sec. VI A, the relatively large value of the coefficient α is chosen to elucidate the difference between the results obtained with the phase-preserving WAPE and NAPE and the WAPEeff and NAPEeff. It follows from Table I that in a motionless medium, the WAPEM coincides with the WAPEeff.

The SPL calculated with these parabolic equations is depicted in Fig. 4, which also shows the SPL calculated with the finite-difference time-domain (FDTD) solution of linearized equations of fluid dynamics. (For details, see the  Appendix.) For the linear sound speed profile, the FDTD approach can be considered as a reference solution.

FIG. 4.

(Color online) SPL for a linear profile of the sound speed produced by a point source located 10 m above a perfectly reflecting, flat surface. The sound frequency is 200 Hz. The subplots correspond to the results obtained with the FDTD solution and parabolic equations as labeled.

FIG. 4.

(Color online) SPL for a linear profile of the sound speed produced by a point source located 10 m above a perfectly reflecting, flat surface. The sound frequency is 200 Hz. The subplots correspond to the results obtained with the FDTD solution and parabolic equations as labeled.

Close modal

It can be observed from Fig. 4 that the phase-preserving WAPE provides the best agreement with the FDTD solution, followed by the WAPEeff, NAPE, and NAPEeff. This conclusion is further illustrated in Fig. 5, where the relative level Δ L calculated with the FDTD approach and parabolic equations is plotted versus the propagation range x for four heights z above the ground. Figure 5 clearly shows that for x 100 m, the results obtained with the phase-preserving WAPE are in a very good agreement with those obtained with the FDTD solution.

FIG. 5.

(Color online) Relative level Δ L versus the propagation range x for a linear profile of the sound speed and four heights z above the surface. The geometry of the problem is the same as that for the SPL in Fig. 4 Different curves correspond to the results obtained with the FDTD solution and parabolic equations.

FIG. 5.

(Color online) Relative level Δ L versus the propagation range x for a linear profile of the sound speed and four heights z above the surface. The geometry of the problem is the same as that for the SPL in Fig. 4 Different curves correspond to the results obtained with the FDTD solution and parabolic equations.

Close modal

The normalized errors ϵ pertinent to the results displayed in Fig. 5 are depicted in Fig. 6. It follows from Fig. 6 that the normalized errors for the phase-preserving WAPE are smaller compared to those for the WAPEeff. Similarly, the phase-preserving NAPE has smaller normalized errors compared to the NAPEeff. The WAPEeff outperforms the NAPE because this section considers sound propagation in a motionless medium for which the effective sound speed approximation is not relevant.

FIG. 6.

(Color online) Normalized errors ϵ between the FDTD solution and results obtained with the parabolic equations from Fig. 5.

FIG. 6.

(Color online) Normalized errors ϵ between the FDTD solution and results obtained with the parabolic equations from Fig. 5.

Close modal

Table III shows the cumulative errors ϵ ¯ pertinent to the normalized errors ϵ in Fig. 6. It follows from Table III that the phase-preserving WAPE has smaller cumulative errors compared to the WAPEeff, as does the phase-preserving NAPE compared to the NAPEeff.

TABLE III.

Cumulative errors ϵ ¯ pertinent to the normalized errors ϵ depicted in Fig. 6.

Equation z = 1 m z = 2 m z = 5 m z = 10 m
NAPE  0.34  0.51  0.68  1.01 
WAPE  0.08  0.09  0.11  0.17 
NAPEeff  0.52  1.09  1.15  1.55 
WAPEeff  0.27  0.41  0.50  0.69 
Equation z = 1 m z = 2 m z = 5 m z = 10 m
NAPE  0.34  0.51  0.68  1.01 
WAPE  0.08  0.09  0.11  0.17 
NAPEeff  0.52  1.09  1.15  1.55 
WAPEeff  0.27  0.41  0.50  0.69 
This section considers a stratification for which the sound speed and medium velocity depend on z. Again, the sound speed profile is given by Eq. (71) but with a much smaller coefficient α = 0.1 s−1. The medium velocity logarithmically increases with height,
(72)
where v 0 = 9 m/s and β = 1 m−1. As in Sec. VI B, the source is located at zs = 10 m above the surface.

Figure 7 shows the relative level Δ L versus the propagation range x for four heights z above the surface. The results correspond to the FDTD solution, the phase-preserving NAPE and WAPE, and the NAPEeff and WAPEeff. For the considered case of a stratified moving medium, the FDTD solution is valid for arbitrary Mach numbers (see Sec. IV D in Ref. 27). It follows from Fig. 7 that the results obtained with the phase-preserving WAPE are very close to the FDTD solution and the phase-preserving NAPE outperforms the NAPEeff and WAPEeff. The results obtained with the latter two parabolic equations often deviate from the reference solution not only at the interference minima and maxima but also at the ranges between them.

FIG. 7.

(Color online) Relative level Δ L versus the propagation range x for a stratified moving medium (Sec. V C) and four heights z above a perfectly reflecting, flat surface. The source is located 10 m above the surface and the sound frequency is 200 Hz. Different curves correspond to the results obtained with the FDTD solution and parabolic equations.

FIG. 7.

(Color online) Relative level Δ L versus the propagation range x for a stratified moving medium (Sec. V C) and four heights z above a perfectly reflecting, flat surface. The source is located 10 m above the surface and the sound frequency is 200 Hz. Different curves correspond to the results obtained with the FDTD solution and parabolic equations.

Close modal

Table IV shows the cumulative errors ϵ ¯ pertinent to the results in Fig. 7. It follows from the table that the phase-preserving WAPE has the smallest errors, followed by the phase-preserving NAPE, WAPEeff, and NAPEeff. The WAPEM was also used to calculate the relative level for the stratification considered in this section. The results obtained are very close to those for the phase-preserving WAPE with the same cumulative errors. This can be explained by the fact that for the stratification considered, the variation in the sound speed is relatively small.

TABLE IV.

Cumulative errors ϵ ¯ pertinent to the results depicted in Fig. 7.

Equation z = 1 m z = 2 m z = 5 m z = 10 m
NAPE  0.09  0.08  0.13  0.13 
WAPE  0.04  0.02  0.02  0.05 
NAPEeff  0.18  0.80  0.56  0.92 
WAPEeff  0.05  0.14  0.17  0.21 
Equation z = 1 m z = 2 m z = 5 m z = 10 m
NAPE  0.09  0.08  0.13  0.13 
WAPE  0.04  0.02  0.02  0.05 
NAPEeff  0.18  0.80  0.56  0.92 
WAPEeff  0.05  0.14  0.17  0.21 

This article considered the phase-preserving EWAPEs, WAPEs, and NAPEs for sound propagation in moving media. These equations [Eqs. (9), (17), and (13), respectively] were derived in the high-frequency approximation while omitting the derivatives of the sound speed and medium velocity. Within the ranges of their applicability, the EWAPE, WAPE, and NAPE exactly describe the phase of sound waves and are valid for arbitrary variations in the sound speed and arbitrary (subsonic) Mach numbers. These equations also correctly describe sound propagation in a stratified moving medium within the ranges of their applicability.

Although variations in sound speed and wind Mach number are often relatively small, omitting the corresponding second-order terms can result in significant phase errors for long propagation ranges. Therefore, it is preferable to use the phase-preserving WAPE and NAPE, at least in some practical applications.

This article further showed that the WAPE given by Eq. (44) in Ref. 15 can also correctly describe the phase of a sound wave but might require more terms in the Padé series expansion. Moreover, that equation involves more terms than the phase-preserving WAPE given by Eq. (17) in this article. The WAPEs from Refs. 2, 13, and 14 preserve the phase of sound waves only approximately for small variations in sound speed and Mach numbers. Previous NAPEs from the literature (e.g., Refs. 2 and 15) correctly describe the phase of sound waves only for small variations in the sound speed and/or low Mach numbers.

Numerical simulations showed that the results obtained with the phase-preserving WAPE and NAPE agree better with the analytical and reference solutions than those with the parabolic equations based on the effective sound speed approximation (the WAPEeff and NAPEeff), which are often used in the literature. Limitations of this approximation have been reported elsewhere, e.g., Refs. 2, 15, and 28. There are, of course, problems for which the effective sound speed approximation might be suitable. However, because the complexity of numerical implementation of the phase-preserving WAPE and NAPE is about the same as that for the WAPEeff and NAPEeff, it makes sense to use the phase-preserving equations.

This research was partially funded by the U.S. Army Engineer Research and Development Center (ERDC) basic research program. Permission to publish was granted by Director, Cold Regions Research and Engineering Laboratory. Any opinions expressed in this article are those of the authors and are not to be construed as official positions of the funding agency or the U.S. Department of the Army unless so designated by other authorized documents. This work was also performed within the framework of the LABEX CeLya (Grant No. ANR-10-LABX-0060) of Université de Lyon within the program Investissements d'Avenir (Grant No. ANR-16-IDEX-0005), operated by the French National Research Agency (ANR).

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

This appendix briefly outlines the FDTD approach employed in Secs. VI B and VI C. The approach is based on solution of Eqs. (17) and (18) in Ref. 27, which, for 2D sound propagation in a stratified moving medium, become
(A1)
where t is time and (wx,wz) are the components of the acoustic particle velocity. The source term is given by
(A2)
where q 0 = 1 m2 s−1 is the source flow rate and B = 0.24 m is its radius.

The FDTD solution uses the optimized fourth-order finite-difference technique (see Ref. 29 for details). The mesh grid is uniform with a grid step of 0.1 m, which corresponds to 17 points per wavelength for f = 200 Hz. The Courant-Friedrichs-Lewy number is set to 0.5, yielding a time step 1.47 × 10 4 s. The FDTD simulation runs for 1.54 s.

The root mean square (rms) sound pressure is determined from the FDTD solution using the last three periods of the sound signal. For comparison with the parabolic equations, the rms pressure is normalized using the factor q 0 ω ϱ exp { ω 2 B 2 / [ 4 c 2 ( z s ) ] }.

1.
E.
Salomons
,
Computational Atmospheric Acoustics
(
Kluwer Academic
,
Dordrecht
,
2001
), p.
335
.
2.
V. E.
Ostashev
and
D. K.
Wilson
,
Acoustics in Moving Inhomogeneous Media
, 2nd ed. (
CRC Press
,
Boca Raton, FL
,
2015
),
521
pp.
3.
J. F.
Lingevitch
,
M. D.
Collins
,
D. K.
Dacol
,
D. P.
Drob
,
J. C. W.
Rogers
, and
W. L.
Siegmann
, “
A wide angle and high Mach number parabolic equation
,”
J. Acoust. Soc. Am.
111
,
729
734
(
2002
).
4.
C.
Khodr
,
M.
Azarpeyvand
, and
D. N.
Green
, “
An iterative three-dimensional parabolic equation solver for propagation above irregular boundaries
,”
J. Acoust. Soc. Am.
148
(
2
),
1089
1100
(
2020
).
5.
C. M.
Nyborg
,
K.
Bolin
,
I.
Karasalo
, and
A.
Fischer
, “
An inter-model comparison of parabolic equation methods for sound propagation from wind turbines
,”
J. Acoust. Soc. Am.
154
,
1299
1314
(
2023
).
6.
F. D.
Tappert
, “
The parabolic approximation method
,” in
Wave Propagation in Underwater Acoustics.
Lecture Notes in Physics No. 70, edited by
J. B.
Keller
and
J. S.
Papadakis
(
Springer
,
New York
,
1977
).
7.
F. B.
Jensen
,
W. A.
Kuperman
,
M. B.
Porter
, and
H.
Schmidt
,
Computational Ocean Acoustics
, 2nd ed. (
Springer
,
New York
,
2011
), pp.
457
529
.
8.
M. D.
Collins
, “
A split-step Padé solution for the parabolic equation method
,”
J. Acoust. Soc. Am.
93
,
1736
1742
(
1993
).
9.
F.
Dagrau
,
M.
Rénier
,
R.
Marchiano
, and
F.
Coulouvrat
, “
Acoustic shock wave propagation in a heterogeneous medium: A numerical simulation beyond the parabolic approximation
,”
J. Acoust. Soc. Am.
130
,
20
32
(
2011
).
10.
L.-J.
Gallin
,
M.
Rénier
,
É.
Gaudard
,
T.
Farges
,
R.
Marchiano
, and
F.
Coulouvrat
, “
One-way approximation for the simulation of weak shock wave propagation in atmospheric flows
,”
J. Acoust. Soc. Am.
135
,
2559
2570
(
2014
).
11.
P. V.
Yuldashev
,
I. S.
Mezdrokhin
, and
V. A.
Khokhlova
, “
Wide-angle parabolic approximation for modeling high-intensity fields from strongly focused ultrasound transducers
,”
Acoust. Phys.
64
(
3
),
309
319
(
2018
).
12.
R.
Janaswamy
, “
A curvilinear coordinate-based split-step parabolic equation method for propagation predictions over terrain
,”
IEEE Trans. Antennas Propagat.
46
(
7
),
1089
1097
(
1998
).
13.
V. E.
Ostashev
,
D.
Juvé
, and
P.
Blanc-Benon
, “
Derivation of a wide-angle parabolic equation for sound waves in inhomogeneous moving media
,”
Acta Acust. Acust.
83
(
3
),
455
460
(
1997
).
14.
P.
Blanc-Benon
,
L.
Dallois
, and
D.
Juvé
, “
Long range sound propagation in a turbulent atmosphere within the parabolic approximation
,”
Acta Acust. Acust.
87
(
6
),
659
669
(
2001
).
15.
V. E.
Ostashev
,
D. K.
Wilson
, and
M. B.
Muhlestein
, “
Wave and extra-wide-angle parabolic equations for sound propagation in a moving atmosphere
,”
J. Acoust. Soc. Am.
147
(
6
),
3969
3984
(
2020
).
16.
D. K.
Wilson
,
M. J.
Shaw
,
V. E.
Ostashev
,
M. B.
Muhlestein
,
M. E.
Swearingen
, and
S. L.
McComas
, “
Numerical modeling of mesoscale infrasound propagation in the Arctic
,”
J. Acoust. Soc. Am.
151
(
1
),
138
157
(
2022
).
17.
B.
Kayser
,
D.
Mascarenhas
,
B.
Cotté
,
D.
Ecotiere
, and
B.
Gauvreau
, “
Validity of the effective sound speed approximation in parabolic equation models for wind turbine noise propagation
,”
J. Acoust. Soc. Am.
153
(
3
),
1846
1854
(
2023
).
18.
D.
Mascarenhas
,
B.
Cotté
, and
O.
Doaré
, “
Propagation effects in the synthesis of wind turbine aerodynamic noise
,”
Acta Acust.
7
(23),
1
16
(
2023
).
19.
J.
Colas
,
A.
Emmanuelli
,
D.
Dragna
,
P.
Blanc-Benon
,
B.
Cotte
, and
R. J. A. M.
Stevens
, “
Wind turbine sound propagation: Comparison of a linearized Euler equations model with parabolic equation methods
,”
J. Acoust. Soc. Am.
154
(
3
),
1413
1426
(
2023
).
20.
A. D.
Pierce
, “
Wave equation for sound in fluids with unsteady inhomogeneous flow
,”
J. Acoust. Soc. Am.
87
(
6
),
2292
2299
(
1990
).
21.
V. E.
Ostashev
,
M. B.
Muhlestein
, and
D. K.
Wilson
, “
Extra-wide-angle parabolic equations in motionless and moving media
,”
J. Acoust. Soc. Am.
145
(
2
),
1031
1047
(
2019
).
22.
H. S.
Wall
,
Analytic Theory of Continued Fractions
(
Chelsea Publishing
,
New York
,
1973
).
23.
M. D.
Collins
, “
Applications and time-domain solution of higher-order parabolic equations in underwater acoustics
,”
J. Acoust. Soc. Am.
86
(
3
),
1097
1102
(
1989
).
24.
L.
Fishman
and
J. J.
McCoy
, “
Derivation and application of extended parabolic wave theories. I. The factorized Helmholtz equation
,”
J. Math. Phys.
25
,
285
296
(
1984
).
25.
A.
Bamberger
,
B.
Engquist
,
L.
Halpern
, and
P.
Joly
, “
Higher order paraxial wave equation approximations in heterogeneous media
,”
SIAM J. Appl. Math.
48
,
129
154
(
1988
).
26.
F.
Collino
, “
Perfectly matched absorbing layers for the paraxial equations
,”
J. Comput. Phys.
131
(
1
),
164
180
(
1997
).
27.
V. E.
Ostashev
,
D. K.
Wilson
,
L.
Liu
,
D. F.
Aldridge
,
N. P.
Symons
, and
D. H.
Marlin
, “
Equations for finite-difference, time-domain simulation of sound propagation in moving inhomogeneous media and numerical implementation
,”
J. Acoust. Soc. Am.
117
(
2
),
503
517
(
2005
).
28.
J.
Assink
,
R.
Waxler
, and
D.
Velea
, “
A wide-angle high Mach number modal expansion for infrasound propagation
,”
J. Acoust. Soc. Am.
141
(
3
),
1781
1792
(
2017
).
29.
D.
Dragna
and
P.
Blanc-Benon
, “
Towards realistic simulations of sound radiation by moving sources in outdoor environments
,”
Int. J. Aeroacoust.
13
(
5-6
),
405
426
(
2014
).