© The Author(s), Published by EDP Sciences, 2025

1 Introduction

Although the negative impact of noise pollution on human health is established [1], the effects of anthropogenic noise on biodiversity are under study [2, 3]. Researching the impact of noise on wildlife communications (mainly birds and mammals) brings together various scientific fields and leads to cooperations. The study of sound propagation in forests is part of this effort because it can support in-situ recordings of sound signals made by eco acousticians [4]. To better address this topic, several long-range outdoor sound propagation models are available and provide quantitative noise predictions. Usually, noise maps are produced using physics-based engineering methods which present limitations in specific cases, e.g. for porous ground and/or diffuse field [5].

More generally, outdoor sound environment modeling is frequently performed with a large panel of numerical methods, presenting various advantages and limitations, either wave-based methods in the time domain (e.g. FDTD) or wave-based methods in the frequency domain (e.g. BEM, PE) or analytical methods (e.g. Weyl–Van der Pol formulation) or ray-tracing methods [6–8]. Depending on the specificities (symmetry, frequency range...) of the scenario to be modeled, frequency domain, time domain or modal decomposition methods can be chosen for their assets [8, 9]. Among these methods, the Boundary Element Method (BEM) stands out for its efficiency in boundary-dominated scenarios, as it avoids numerical dispersion errors within the volume and requires discretization of only the boundaries, reducing the dimensionality of the problems by one. However, the Transmission Line Matrix (TLM) method is considered to model complex sound environments with a high number of reflections (such as with tree trunks in forests) using impedance boundary conditions. Its computational cost scales with the size of the domain, but it is better suited to handle broadband scenarios with in-volume complexities. The low order characteristic of this scheme and its implementation makes it a potential candidate to solve the wave equation on such scenarios taking into account the backscattered components of the sound field [10, 11]. It has already been used to model realistic outdoor propagation with meteorological effects or multiple reflections by cylinder arrays normal to an impedance plane [12, 13] with source-receiver distances on the order of hundreds of times the minimum modeled wavelength [14].

The main issue of using a low-order spatial integration scheme, such as TLM or centered finite differences in the time domain (FDTD) on a Cartesian mesh is the apparition of numerical dispersion [15]. This pure discretization-induced artifact leads to spurious waves whose phase and group speeds depend on frequency [16]. To solve this problem, there are techniques to avoid or reduce numerical dispersion but all of them have a cost in either efficiency or complexity. For instance, increasing the number of points in a fixed-order FDTD scheme and optimizing the numerical dispersion equation is commonly done in aeroacoustics [17, 18]. For the TLM method, attempts have been made to modify the geometry of the mesh, but it increases the complexity of the implementation [19, 20]. In most cases, the method remains second-order, with the computational complexity similar to that of a matrix-vector product [11]. The effect of the numerical dispersion on TLM results has been highlighted from the free-field case [11]. This paper investigates the performance with respect to a singular reflection to assess the impact of numerical dispersion on the modeling of boundaries. This is a step toward realistic simulations of sound propagation in forests.

The study presented below focuses on the TLM model and its limitations in specific application domains. Instead of increasing the order of the model (i.e. numerical complexity at the boundaries), this paper thoroughly analyzes the effect of numerical dispersion on long-range modeling with specular reflection included. The approach aims to understand how much numerical dispersion affects the results so that the low-order scheme can be used knowingly.

First, a brief overview, complementary to the literature, is given on the theory and the implementation process of a boundary condition in the TLM method. Then the analytical solution of the specific case considered is given, together with some implementation details. Next, the setup of the numerical experiment is presented, and the results are analyzed for three different sound source characteristics. Finally, a theoretical comparison between the TLM model and a common FDTD method is given in the last part.

2 Theoretical background

The TLM method is based on the Huygens-Fresnel wavefront decomposition to model the propagation of sound waves within a fluid as an iterative process of spreading pressure pulses through a Cartesian mesh [21]. The simulated domain is decomposed into transmission lines connecting nodes located by the vector of indices: $\mathbf{r}=(j_1,\cdots ,j_d)$ where d is the dimension of the problem. The spatial step Δℓ between the lines is equal in every Cartesian directions and the time is decomposed into steps Δt such as $t_n=n\mathrm{\Delta }t,n\in \mathbb{N}$. This nomenclature allows writing the incident pressure pulse to the node r along the line m ∈ {1,…,2d} at the time step n as ${\displaystyle {}_n{I}_{\mathbf{r}}^m}$. Similarly, ${\displaystyle {}_n{S}_{\mathbf{r}}^m}$ represents the scattered pulse reflected instantaneously. From their modelling, it is possible to calculate _n P _r , an approximated value of the exact acoustic pressure p(x ₁,…,x _d ,t _n ) taken at the point (x ₁, …, x _d ) at time t _n .

2.1 Relation to the wave equation

Deriving from the TLM model pressure pulses definition, the equivalent numerical pressure scheme can be formulated as follows [11]:

$$\begin{array}{cc}\hfill {}_{n+1}{P}_{\mathbf{r}}{+}_{n-1}{P}_{\mathbf{r}}=& \frac{1}{d}\sum _{m=1}^d{[}_n{P}_{(j_1+{\delta }_{m1},\cdots ,j_d+{\delta }_{\mathit{md}})}\hfill \\ \hfill & {+}_n{P}_{(j_1-{\delta }_{m1},\cdots ,j_d-{\delta }_{\mathit{md}})}],\hfill \end{array}$$(1)

where δ denotes the Kronecker delta. To analyze the order of approximation of this scheme, and to derive the solved wave equation, Taylor expansions are applied to all the terms of equation (1) [14, 21], leading to:

$$\begin{array}{c}\hfill \frac{{\partial }^{2}p}{\partial t^2}-c_{{}_{\text{TLM}}}^2{\mathrm{\nabla }}^{2}p=\mathcal{O}\left(\mathrm{\Delta }{t}^2\right)+\mathcal{O}\left(\frac{\mathrm{\Delta }{\ell }^4}{\mathrm{\Delta }{t}^2}\right),\end{array}$$(2)

where $c_{\text{TLM}}=\frac{\mathrm{\Delta }\ell }{\sqrt{d}\mathrm{\Delta }t}$. This reveals the wave equation and shows that the TLM model is a second-order approximation method in both time and space. It also highlights that the TLM solves the wave equation only if the condition c ² _TLM = c ₀ ² is satisfied. Expanding this condition yields:

$$\begin{array}{c}\hfill \frac{c_0\mathrm{\Delta }t}{\mathrm{\Delta }\ell }=\frac{1}{\sqrt{d}},\end{array}$$(3)

which corresponds with the Courant–Friedrichs–Lewy (CFL) stability condition of the so-called finite difference Leap-Frog scheme [7].

2.2 Numerical dispersion

To assess the model stability, the numerical dispersion relation of the method can be expressed as:

$$\begin{array}{c}\hfill cos(\omega \mathrm{\Delta }t)=\frac{1}{d}\sum _{m=1}^{d}cos\left(k_m\mathrm{\Delta }\ell \right),\forall \mathrm{\Delta }t,\forall \mathrm{\Delta }\ell .\end{array}$$(4)

This equation indicates that the TLM method is unconditionally stable in a homogeneous, non-dissipative medium. However, it also reveals that the model exhibits numerical dispersion along the primary grid directions (e.g., in the 2D case, $\theta =\beta \frac{\pi }{2},\beta \in \mathbb{Z}$, where θ is the angle between the plane wave vector and the horizontal direction of the grid). Equation (4) is equivalent to the dispersion relation of a second-order FDTD scheme centered in space and time. More theoretical aspects of the homogeneous and inhomogeneous TLM model are thoroughly described for the free field case in [11]. The following will focus on the modeling of a reflection boundary condition within the method.

2.3 Reflection condition

To model a simple specular interaction on a reflective boundary, the TLM model uses a pressure reflection coefficient [22]. Compared to the free-field connection laws [11], the formula governing the transmission line between a boundary and the nearest node (Fig. 1) is modified as described in equation (5). Figure 1 shows additionally that the considered node is located at a distance Δℓ/2 from the boundary. This specificity ensures synchronization between the arrival of pulses reflected from a boundary and the regular pulses traveling along transmission lines of length Δℓ.

Figure 1. Representation of the 2D simple boundary reflection at the local node r = (j1, j2) close to a boundary of reflection coefficient R1. On the line m = 1, scattered pulse “S” (orange) at time step n becomes a reflected incident pulse “I” (blue) at next step n+1.

The energetic absorption coefficient α for the boundary is described by:

$$\begin{array}{c}\hfill \alpha =1-{\left|R_m\right|}^2,\end{array}$$(6)

where R_m corresponds to the plane wave reflection coefficient assigned to the boundary perpendicular to the branch m. The condition R_m = 1 describes an infinitely rigid boundary (zero particle velocity), R_m = −1 a free boundary (zero acoustic pressure), and R_m = 0 a completely absorbing boundary.

3 Analytical solution – image source theory

The analytical solution of the pressure field with specular reflection on a perfectly reflective ground surface for an omnidirectional point source is used as a reference to estimate the reliability of the numerical results. The time-domain formulation of the pressure at a point r is relatively straightforward and can be written as [11, 23, 24]:

$$\begin{array}{c}\hfill p(\mathbf{r},t)={\mathrm{TF}}^{-1}\{-i\omega {\rho }_0{\widehat{Q}}_0(w){\widehat{G}}_d(\mathbf{r},\omega )\},\end{array}$$(7)

with TF⁻¹ {} the inverse Fourier transform, ρ₀ the medium density, ${\widehat{Q}}_0(w)$ the mass flow and ${\displaystyle {\widehat{G}}_d(\mathbf{r},\omega )}$ the sum of two Green functions of the free-field given by:

$$\begin{array}{c}\hfill {\widehat{G}}_d(\mathbf{r},\omega )=\frac{i}{4}(H_0^{(1)}\left(k_{\omega }\parallel \mathbf{r}-{\mathbf{r}}_s\parallel \right)+H_0^{(1)}\left(k_{\omega }\parallel \mathbf{r}-{\mathbf{r}}_s^{\prime }\parallel \right))\end{array}$$(8)

$H_0^{(1)}$ being the Hankel function of the first kind, k_w the wavenumber and r_s′ is the image source location.

This formulation of the analytical solution is based on the image source theory. This formulation is valid only for purely reflective boundaries, i.e. with pressure reflection coefficients R_m = 1 (reflective boundary), R_m = −1 (soft boundary) and R_m = 0 (absorbing boundary) [25]. Thus, the values needed to calculate the pressure field p(r,t) are the distances ‖r − r_s‖ and $\Vert \mathbf{r}-{\mathbf{r}}_s^{\prime }\Vert$ between the source S and the receiver and the image source S’ and the receiver, respectively (Fig. 2). In equation (8), the second term H₀⁽¹⁾(k_ω∥r − r_s′∥) accounts for the reflection on the surface, modeling a contribution from an image source in free field.

Figure 2. Numerical experiment setup with the image source theory. The dotted lines represent free boundaries and the orange line a perfectly reflective boundary.

For a simulation setup as in Figure 2, with a boundary condition located at y = 0, r = ‖r − r_s‖ and $r^{\prime }=\Vert \mathbf{r}-{\mathbf{r}}_s^{\prime }\Vert$ become:

$$\begin{array}{c}\hfill \parallel \mathbf{r}-{\mathbf{r}}_s\parallel =\sqrt{{(x-x_{}\text{s})}^2+{(y-y_{}\text{s})}^2},\end{array}$$(9a)

$$\begin{array}{c}\hfill \parallel \mathbf{r}-{\mathbf{r}}_s^{\prime }\parallel =\sqrt{{(x-x_{}\text{s})}^2+{(y+y_{}\text{s})}^2}.\end{array}$$(9b)

However, when these distances are discretized throughout the TLM model, special care must be taken because of the spatial discretization of the domain. Indeed, as shown in Figure 3, the distance to the boundary surface is not exactly Δx_i = j_iΔℓ because the meshing process introduces two artificial nodes to interface with the solver code [26]. The presence of the marker nodes (“outside” and “boundary”) leads to:

Figure 3. Modeling of the boundaries in the TLM implementation. Orange dots: “outside” nodes; black dots: “boundary” nodes; blue dots: “fluid” nodes. The red dotted line represents the boundary as modeled in the concept representation of Figure 1.

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}_i=(j_i-\frac{3}{2})\mathrm{\Delta }\ell .\end{array}$$(10)

Thus, the distances of interest on the grid (noted r and r^′ for clarity) are computed as:

$$\begin{array}{cc}& r_{j_1,j_2}=\sqrt{{\left((j_1-j_{}\text{src})\mathrm{\Delta }\ell \right)}^2+{\left((j_2-j_{}\text{src})\mathrm{\Delta }\ell \right)}^2},\hfill \end{array}$$(11a)

$$\begin{array}{cc}& r_{j_1,j_2}^{\prime }=\sqrt{{\left((j_1-j_{}\text{src})\mathrm{\Delta }\ell \right)}^2+{\left((j_2+j_{}\text{src}-3)\mathrm{\Delta }\ell \right)}^2}.\hfill \end{array}$$(11b)

Reintroducing these distances in equation (8) gives the solution for the analytical pressure on the grid ${}_n{p}_{j_1,j_2}$.

4 Numerical experiment

A numerical experiment is now carried out with the introduction of a perfectly reflective surface. The numerical results from the TLM model are compared to the image source formulation of the analytical solution described in Section 3.

4.1 Simulation setup

The numerical experiment setup displayed in Figure 4 is implemented for two main geometric parameters: the height h between the ground and the source and the maximum distance r_max between the last receiver point and the source. The number of receivers in the array is also determinant since a too-coarse matrix can lead to a lack of precision in the detection of the interference patterns. For the following experiments, a matrix of 100 ×100 receivers is used. Thus, the most important modeling parameter for the analysis is the number of points per wavelength N_ppw = λ/Δℓ which is set to 10 at the maximal frequency f_max of the source emission.

Figure 4. Numerical experiment setup.

4.2 Data processing

Various approaches are available to evaluate the numerical dispersion error. In order to introduce the method, two third-order tensors ${}_n{\mathfrak{n}}_{j_1,j_2}$ and ${}_n{\mathfrak{a}}_{j_1,j_2}$ are introduced to represent numerical (TLM simulations) and analytical results respectively. They are consistent with the formalism used before (Sect. 3): ${}_n{\mathfrak{n}}_{j_1,j_2}$ is the pressure signal at the receiver r = (j₁,j₂) and at the time step n. The tensors are filled with real values that are the analytical and modeled pressure, their size is the number of time steps and the number of microphones along the two axes of the mesh.

4.2.1 Attenuation relative to a reference receiver

The SPL attenuation relative to a reference receiver is used as a metric to assess the discrepancies between simulated pressure signals and their corresponding analytical solutions. Errors are evaluated using sound attenuation levels rather than acoustic pressure to reduce sensitivity to small pressure variations that can cause disproportionate local errors, following standard practice in environmental acoustics. For every line of receivers, the reference receiver is the closest to the source. The SPL attenuation is calculated as:

$$\begin{array}{c}\hfill A_{j_1,j_2}\left(\mathfrak{n}\right)=10{log}_{10}\left(\frac{{\sum }_{n=0}^{N-1}{\left({}_n{\mathfrak{n}}_{j_1,j_2}\right)}^2}{{\sum }_{n=0}^{N-1}{\left({}_n{\mathfrak{n}}_{0,j_2}\right)}^2}\right).\end{array}$$(12)

Equation (12) can also be written for ${}_n{\mathfrak{a}}_{j_1,j_2}$. To compare the numerical results to the analytical solutions, the following SPL absolute error is used:

$$\begin{array}{c}\hfill {\epsilon }_{j_1,j_2}=|A_{j_1,j_2}(\mathfrak{n})-A_{j_1,j_2}(\mathfrak{a})|[\text{dB}].\end{array}$$(13)

For a better normalization process, a polar array of receivers is chosen to ensure that all the receivers of an arc are at the same distance from the source.

5 Results

Since the numerical dispersion affects more the high frequencies during the sound propagation [11] the analysis is performed by comparing different sources, regarding distinct spectral distributions. For a better understanding of the results, the number of points per wavelength N_ppw is presented for each source. An overview of one simulation is displayed in Figure 5 for a Gaussian pulse.

Figure 5. Sound pressure field for a Gaussian source at the end of the TLM simulation. The whole calculation domain is represented. At the last iteration of the calculation, reflected and incident waves have passed through all the receivers (blue dots) in the network.

In the following subsections, the absolute errors εj₁,j₂ on the SPL attenuation at each receiver are displayed (cf. Eq. (13)). The objective is to determine whether dispersion-induced delayed interferences occur in the receivers’ matrix. Indeed, if there is an interference pattern for one of the delayed frequencies at a receiver, a lack of energy compared to the analytical solution will appear in the simulated signals. The following results are presented for f_max = 500 Hz, h = 2m and the sound speed c₀ = 344.24 m.s⁻¹. The simulation time is t_sim = 4.19e–01 s with Δℓ = 6.88e–02 m and Δt = 4.14e–04 s. The height of the source, h is equal for all the simulations since it has been shown to affect mostly the interference patterns between the direct and reflected wavefront, leading to similar order of magnitude in the errors [14].

Figure 6. (a) Time signal according to time (top) and spectral distribution (bottom) of a Gaussian pulse compared to Nppw (orange dashed line) according to normalized frequency. (b) Polar map of the SPL absolute error εj1,j2 (Eq. (13)).

5.1 Gaussian source

In this case, the source signal is implemented according to:

$$\begin{array}{c}\hfill S({ṟ}_s,n)=exp(-{\pi }^2{(\frac{f_{}\text{max}}{2}n\mathrm{\Delta }t-1)}^2).\end{array}$$(14)

As it can be observed in Figure 6, most of the spectral components of the Gaussian source are over-discretized and those corresponding to f ≥ f_max and N_ppw ≤ 10 represent less than 2% of the signal energy. On the vertical and horizontal axes, negligible absolute errors appear between the analytical and the numerical results. In addition, an interference pattern of the error is observable with a “ray” of errors in the areas where the dispersed end of the incident wavefront interacts with its reflection.

Figure 7. Statistical distribution of the absolute errors εj1,j2 for the Gaussian pulse case.

In Figure 7, the statistical distribution of the SPL error is represented without considering the spatial locations. It confirms that for an over-discretized source signal, the error related to numerical dispersion is negligible (95% of the errors are below 0.03 dB).

5.2 Shifted Gaussian source

In this case, the source signal is implemented according to:

$$\begin{array}{c}\hfill S({ṟ}_s,n)=exp(-{\pi }^2{(1.2\times f_{}\text{max}n\mathrm{\Delta }t-1)}^2).\end{array}$$(15)

The source has a broader spectrum than in the previous case. The excitation signal is therefore less over-sampled and the impact on the results is visible on the Figure 8b. The ray of interference also appears, but the error magnitudes are higher than for the regular Gaussian case because of the lower discretization of the source spectral components.

Figure 8. (a) Time signal according to time (top) and spectral distribution (bottom) of a shifted Gaussian pulse compared to Nppw (orange dashed line) according to normalized frequency. (b) Map of the SPL absolute error εj1,j2.

Indeed, the maximum error is of 0.44 dB and 95% of the errors in the considered quarter circle are below 0.1 dB Figure 9. This shift of approximately a factor ten in the errors confirms that the discretization of the modeled wave packet influences the numerical dispersion effect on long-range results (r_max ≈ 150 × λ_min here).

Figure 9. Statistical distribution of the absolute errors εj1,j2 for the shifted Gaussian pulse case.

5.3 Windowed sinusoidal source

In this case, the source signal is implemented according to:

$$\begin{array}{c}\hfill {sin}_{}\text{w}[n]={\text{w}}_{}\text{hann}\times sin\left(2\pi f_{}\text{max}n\mathrm{\Delta }t\right)\end{array}$$(16)

with W_hann a Hann window of duration 0.15 × t_sim. The two previous cases implied broadband sources. Now, a windowed sine is used as it is pseudo-harmonic. As shown in Figure 10a, most of the source energy is focused around the frequency f_max and has a 10-point spatial discretization. Thus, the direct wave is more likely to interfere with the reflected one. With this spectral distribution, significant interference-induced errors are generated by the numerical dispersion (up to 11.06 dB as shown in Fig. 10b). Indeed, at some receivers, the energy of the received signal is either increased or lost due to artificial interference patterns induced by the TLM model. Comparing Figures 10b, 8b and 6b also leads to the conclusion that the spectral distribution of the source has an impact on the spatial distribution of the error: the wider the spectral distribution of the source is, the wider the spatial distribution of the error.

Figure 10. (a) Time signal according to time (top) and spectral distribution (bottom) of a windowed sine compared to Nppw (orange dashed line) according to normalized frequency. (b) Polar map of the SPL absolute error εj1,j2 (Eq. (13)).

In this case, the histogram in Figure 11 presents a broader occurrence distribution of the errors compared to Figures 7 and 9. The most critical value is that 5% of the errors are in the interval [5–11 dB]. This value suggests that depending on the source definition, more than ten points per wavelength are needed to model the propagation of wave packets with the TLM model. To make the concept of misplaced interference patterns clearer, the next Section 5.4 shows in more detail the attenuation along receiver lines of interest.

Figure 11. Statistical distribution of the absolute errors εj1,j2 for the windowed sine case.

5.4 Attenuation along line

To better explain the maxima of errors presented in the windowed sine case, Figure 12 displays the SPL attenuation relative to the reference receiver for three lines of receivers: θ = [4°, 10°, 15°]. The plots in Figure 12 help visualize better the lines of receivers on the previous color maps.

Figure 12. Numerical SPL attenuation relative to a reference receiver according to the propagation distance r (blue line), compared to the analytical solution (orange dashed line) for three source-receiver angles: (a) θ = 4°, (b) θ = 10° and (c) θ = 15°. Source located at hs = 2m and emitting a windowed sinusoidal signal at fmax = 500 Hz (Nppw = 10). On the error maps on the left, red lines indicate the receiver lines corresponding to the SPL attenuation graphs on the right.

The “Sawtooth” effect on the attenuation profiles comes from the Cartesian discretization of the numerical domain. It does not affect the results as it is present for both analytical and numerical solutions. In Figure 12a and 12c related to θ = 4° and θ = 15°, it is observed that maxima of errors are due to misplacement of attenuation peaks along the source-receiver distance. The θ = 15° line is the one that presents the largest gap because the third attenuation peak is simulated closer to the source than the analytical solution.

6 Comparison to FDTD schemes

Since the FDTD is the closest numerical method to the TLM in terms of discretization of partial differential equations, it seems relevant to compare the error on the modeled axial sound speeds from the two models. As highlighted by the stability analysis of the TLM model, the numerical dispersion is maximal along the main directions of the mesh [11]. Thus, considering θ the angle between the plane wave vector and the horizontal direction of the mesh, the numerical phase speed in a 2D Cartesian grid c_ph can be expressed as a function of N_ppw on the main directions of the mesh ($\theta =p\frac{\pi }{2}$, $p\in \mathbb{Z}$), as

$$\begin{array}{c}\hfill c_{}\text{ph}=c\frac{N}{\pi \sqrt{2}}arccos\left[\frac{1}{2}(cos\left(\frac{2\pi }{N_{}\text{ppw}}\right)+1)\right],\end{array}$$(17)

with c the modeled sound speed in the propagation medium.

Similar expressions can be formulated for the group speed, for 3D cases or other “Taylor-expansions series based” numerical methods [16]. In Figure 13, the phase and group speeds associated with centered finite difference schemes are compared to the ones related to the TLM method. From these figures, it is straightforward that increasing the method order reduces the dispersion errors and that the phenomenon is stronger on the group speed. For instance, more than 30 points per wavelength would be needed for the TLM to get the accuracy that the fourth-order FDTD scheme has at N = 10.

Figure 13. (a) Effective speeds compared to the physical constant; (b) Relative errors on the phase and group velocities [%] related to the model anisotropy depending on the number Nppw. “DFOD2” and “DFOD4” stand for the second and fourth-order centered finite difference schemes, respectively.

As mentioned in Section 2, it is straightforward on Figure 13b that the TLM and the second-order FDTD numerical phase speeds are completely equivalents. Hence, the same conclusion is applied to the numerical group speed. Additionally, Figure 13b shows that the same number of points per wavelength would result in the same dispersion error for the two methods, so the numerical cost of using the TLM or a second order FDTD is the same. The comparison results with fourth-order FDTD confirm that for in-volume modeling, the TLM method is not the most performant method for universal applications. However, this statement can be moderated by the fact that a lower-order scheme implies less implementation complexity on the boundary conditions. Thus, when numerous scatterers are considered, second-order schemes could be attractive for their computational costs.

7 Conclusion

The limitations of the TLM model for long-range sound propagation simulations over a rigid ground were investigated thanks to numerical experiments and a theoretical comparison. A polar array of virtual receivers has been used to study the interactions between simulated direct sound waves and their reflection. The results from three distinct source signals emphasize the importance of the number of points per wavelength as a simulation parameter. The presence of a perfect specular reflection in the model shows that numerical dispersion can affect the simulated sound pressure levels. In specific cases of reflection, misplaced simulated interferences appear due to the wrongly approximated group speed of some spectral components of the wave packet. Relatedly, the link between the spatial distribution of the errors and the spectral distribution of the source has been highlighted. Finally, a comparison between the TLM model and a higher-order centered FDTD method stated that for general in-volume modeling, the TLM is not the more efficient method against numerical dispersion as it requires a finer discretization of the mesh.

However, since a higher-order spatial integration requires a higher-order integration scheme for the boundary conditions, the idea of using a low-order scheme to model a complex environment needs further analysis. In addition, the use of a non-dissipative method such as TLM in multiple scatterers and source scenarios where the dispersion-induced error could be mitigated seems to be a promising and computationally efficient approach. A comparison with in-situ measurements under real outdoor conditions should be made to characterize the reliability of the simulated sound pressure levels in environmental acoustics applied cases. In addition, further research to increase the accuracy of the TLM model from second to fourth order could help in the effort to compare its performance with finite difference schemes.

Acknowledgments

The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d’Avenir) and the CPER Alsacalcul/Big Data.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability statement

Data are available on request from the authors.

References

All Figures

	Figure 1. Representation of the 2D simple boundary reflection at the local node r = (j1, j2) close to a boundary of reflection coefficient R1. On the line m = 1, scattered pulse “S” (orange) at time step n becomes a reflected incident pulse “I” (blue) at next step n+1.
In the text

	Figure 2. Numerical experiment setup with the image source theory. The dotted lines represent free boundaries and the orange line a perfectly reflective boundary.
In the text

	Figure 3. Modeling of the boundaries in the TLM implementation. Orange dots: “outside” nodes; black dots: “boundary” nodes; blue dots: “fluid” nodes. The red dotted line represents the boundary as modeled in the concept representation of Figure 1.
In the text

	Figure 4. Numerical experiment setup.
In the text

	Figure 5. Sound pressure field for a Gaussian source at the end of the TLM simulation. The whole calculation domain is represented. At the last iteration of the calculation, reflected and incident waves have passed through all the receivers (blue dots) in the network.
In the text

	Figure 6. (a) Time signal according to time (top) and spectral distribution (bottom) of a Gaussian pulse compared to Nppw (orange dashed line) according to normalized frequency. (b) Polar map of the SPL absolute error εj1,j2 (Eq. (13)).
In the text

	Figure 7. Statistical distribution of the absolute errors εj1,j2 for the Gaussian pulse case.
In the text

	Figure 8. (a) Time signal according to time (top) and spectral distribution (bottom) of a shifted Gaussian pulse compared to Nppw (orange dashed line) according to normalized frequency. (b) Map of the SPL absolute error εj1,j2.
In the text

	Figure 9. Statistical distribution of the absolute errors εj1,j2 for the shifted Gaussian pulse case.
In the text

	Figure 10. (a) Time signal according to time (top) and spectral distribution (bottom) of a windowed sine compared to Nppw (orange dashed line) according to normalized frequency. (b) Polar map of the SPL absolute error εj1,j2 (Eq. (13)).
In the text

	Figure 11. Statistical distribution of the absolute errors εj1,j2 for the windowed sine case.
In the text

Figure 12. Numerical SPL attenuation relative to a reference receiver according to the propagation distance r (blue line), compared to the analytical solution (orange dashed line) for three source-receiver angles: (a) θ = 4°, (b) θ = 10° and (c) θ = 15°. Source located at hs = 2m and emitting a windowed sinusoidal signal at fmax = 500 Hz (Nppw = 10). On the error maps on the left, red lines indicate the receiver lines corresponding to the SPL attenuation graphs on the right.

In the text

	Figure 13. (a) Effective speeds compared to the physical constant; (b) Relative errors on the phase and group velocities [%] related to the model anisotropy depending on the number Nppw. “DFOD2” and “DFOD4” stand for the second and fourth-order centered finite difference schemes, respectively.
In the text