The influence of topography on sonic boom propagation is investigated. The full two-dimensional Euler equations in curvilinear coordinates are solved using high-order finite-difference time-domain techniques. Simple ground profiles, corresponding to a terrain depression, a hill, and a sinusoidal terrain, are examined for two sonic boom waves: a classical N-wave and a low-boom. Ground reflection of the sonic boom is affected by elevation variations: a concave ground profile induces compression, which tends to increase the peak pressure in particular, while the opposite is true for convex elevation variations, which lead to expansion and a reduction in peak pressure. The reflected boom is then strongly altered. Furthermore, a sufficiently concave topography can cause focal zones, which generate extra contributions at ground level in the form of U-waves in addition to the reflected wave. This mechanism has the largest effect on waveforms at ground level. The variations of standard metrics are of a few dBs compared to a flat ground for both sonic boom waves, and they are notably greater for the terrain depression than for the hill. Finally, in the case of a sinusoidal terrain, the pressure waveforms are composed of multiple arrivals due to successive focal zones.

## I. INTRODUCTION

Several supersonic business jet projects are under development (Sun and Smith, 2017), with entry-into-service planned in the coming years. This should mark the beginning of a second era of commercial supersonic flight. However, due to the sonic boom, overland commercial supersonic flight remains banned in most countries. This has serious consequences on the possible demand for supersonic transport and on its commercial success (Liebhardt *et al.*, 2020). For public acceptance, it is thus important to reduce sonic boom annoyance by designing low-boom aircraft, but also to be able to accurately predict the sonic boom that will be generated at ground level by this new generation of supersonic aircraft.

To this end, meteorological and ground effects have to be taken into account. While several studies are devoted to the former, ground effects have not received much attention. It is thus common in prediction models to account for ground reflection by multiplying the sonic boom incident to the ground by a constant factor, usually equal to 1.9 [for instance, in recent papers Ueno *et al.* (2017) and Phillips and West (2019)]. However, for specific configurations such as focused booms and booms near the carpet cut-off, studies (Coulouvrat, 2002; Rendón and Coulouvrat, 2005) have shown that ground absorption significantly affects the sonic boom signature with a reduction of the peak overpressure and an increase in rise time. Similar effects have been observed for sea roughness (Boulanger and Attenborough, 2005). The influence of irregular terrain on sonic boom reflection has been studied little, even though the dependence of sonic boom at the ground level on the surrounding topography can be expected. In their comprehensive review, Maglieri *et al.* (2014) only refer to the work of Dini and Lazzeretti (1969) on this subject. Therefore, there is a need to investigate how a sonic boom is affected by the topography at ground level using modern prediction methods.

The prediction of the sonic boom remains mainly based on the analytical or numerical calculation of the flow field near the aircraft, coupled with efficient geometric acoustic methods to calculate the propagation through the atmosphere to the ground [see, e.g., Cleveland *et al.* (1996), Rallabhandi (2011), and Kanamori *et al.* (2018)]. Ray-tracing was thus the only method employed for far-field prediction in the second and third AIAA Sonic Boom Prediction Workshops (AIAA, 2017, 2020; Rallabhandi and Loubeau, 2019). One-way wave equation methods are favored when propagation in a turbulent medium is examined (Blanc-Benon *et al.*, 2002; Yuldashev *et al.*, 2017). State-of-the-art simulations (Luquet *et al.*, 2019; Stout and Sparrow, 2019) are now able to consider a three-dimensional geometry.

Recently, the first direct simulations of the entire flow field around a supersonic body were performed by Yamashita and Suzuki (2016) with the aim of computing the generation of the sonic boom and its propagation in the atmosphere in a single simulation. While appealing, this approach is computationally demanding, especially to capture the high-frequency characteristics of the sonic boom wave. A grid refinement strategy is thus required near the shocks, which is relevant for a steady sonic boom but seems to be hardly applicable as soon as atmospheric turbulence or irregular terrain are accounted for. Such simulations have been limited to a stratified atmosphere until now, and ground reflection was neglected.

Numerical simulations based on the fluid mechanics equations are, however, also attractive for sonic boom propagation in the planetary boundary layer. Compared to ray-tracing methods, they readily handle diffraction and do not suffer from the angular limitations of one-way wave equations. This is of particular importance for propagation over irregular terrain because shadow zones, caustics, and focal zones and backward diffraction are significant in this case.

The objectives of this study are to investigate the influence of irregular terrain on sonic boom reflection and to identify the physical mechanisms at play. With this aim, numerical simulations based on the Euler equations are carried out. For ease of analysis, simple ground profiles are considered. Two sonic boom signatures are examined: a classical N-wave, as generated by the former generation of civil supersonic aircraft, and a low-boom wave, which is expected for the next generation.

The paper is organized as follows. In Sec. II, the propagation model is detailed, including the equations, the boundary conditions, and the injection of the sonic boom wave. The configurations investigated are then described in terms of both input and computational parameters and a mesh convergence study is performed. The results for isolated terrain irregularities, namely, a terrain depression and a hill, are shown in Sec. III. The change in the sonic boom signature along the ground and in the corresponding metrics are discussed. Section IV focuses on a sinusoidal terrain and cumulative effects are examined. Concluding remarks are given in Sec. V.

## II. SETUP OF NUMERICAL SIMULATIONS

### A. Presentation of numerical methods

Sonic boom propagation above non-flat terrain is considered in a two-dimensional geometry, corresponding to the vertical plane through the axis of the aircraft (see Fig. 1). The Euler equations (Anderson, 1995) are solved using numerical simulation, hence neglecting viscous and thermal diffusion terms. A curvilinear coordinate transformation is applied to allow for variations of the ground profile, as illustrated in the computational domain representation in Fig. 2. In what follows, *ρ* denotes the density, *p* the pressure, *ρe* the energy density, *u* and *v* the horizontal and vertical components of velocity. Their mean values are indicated by the subscript 0. After transformation from the Cartesian to curvilinear coordinate systems $(e\u2192x,e\u2192z)$ and $(e\u2192\xi ,e\u2192\eta )$, the Euler equations in the conservative form are written as

where *t* is time, $U=[\rho \u2009\rho u\u2009\rho v\u2009\rho e]T$ is the vector of unknowns, and **E** and **F** are the following vectors:

Finally, in Eq. (1), the Jacobian of the transformation is denoted by $J=(\xi x\eta z\u2212\xi z\eta x)$, where the notation $ij=\u2202i/\u2202j$ is used for the Jacobian terms. The equations are closed using the definition of energy density for a perfect gas,

with *γ* the heat capacity ratio.

The terrain-following transformation proposed by Gal-Chen and Sommerville (1975) is employed,

where *h* is the ground profile and $zmax$ is the maximal altitude of the computational domain. This formulation assumes the terrain is sufficiently smooth, without any slope discontinuity, so that topographies such as cliffs or buildings cannot be considered with this version of the code. The first line of the mesh follows the ground profile, as represented in Fig. 2. The following lines in the direction parallel to the ground then slowly adapt towards a horizontal line at the top of the domain. The mesh lines in the *z*-direction remain at constant axial positions, which leads to the simplification of Eq. (1) with $\xi x=1$ and $\xi z=0$.

High-order, optimized numerical schemes are used to accurately capture acoustic fluctuations. A fourth order Runge-Kutta scheme including six sub-steps is used in time (Berland *et al.*, 2006). The equations are discretized in space using fourth order finite difference schemes over 11 points (Bogey and Bailly, 2004; Berland *et al.*, 2006). Selective filters of sixth order for the interior points (Bogey *et al.*, 2009) and of second order for the boundary points (Berland *et al.*, 2006) are used. They aim at removing grid-to-grid oscillations that are not solved by the finite-difference schemes and can lead to numerical instabilities. A shock-capturing method (Bogey *et al.*, 2009; Sabatini *et al.*, 2016) is applied to handle shocks, that can be generated during sonic boom propagation. It consists in adding artificial viscosity only near the shocks. Note that as an explicit time-integration scheme is used, the time step and the grid size are related via the Courant-Friedrichs-Lewy (CFL) number.

The grid size in both directions d*ξ* and d*η* remains constant in Cartesian coordinates throughout the computational domain and is the same, i.e., $d\xi =d\eta $. While the mesh coordinates and Jacobian terms are computed over the whole domain, flow variables are only simulated in a moving frame, allowing one to reduce the computational cost significantly. As shown in Fig. 1, this moving frame corresponds to a smaller domain which advances in time along the ground profile and at the flight Mach number *M*, in order to follow the sonic boom wave emitted by the supersonic aircraft.

A rigid wall boundary condition is set at the bottom of the domain, with normal velocity equal to zero. While ground absorption is expected to attenuate the high frequency content of the reflected boom significantly, increasing the rise time and reducing the peak pressure, it is not taken into consideration in order to focus on topographic effects. Note that specific numerical methods suitable for high-order time-domain solvers have been developed to account for locally- or extended-reacting ground surfaces (Dragna *et al.*, 2015; Troian *et al.*, 2017). Modeling ground absorption on sonic boom reflection is thus readily feasible and is left for future work. As previously explained, the use of curvilinear coordinates allows this boundary to be curved to take into account variations of terrain elevation. The ground profile is read from a file, so that any curved complex terrain may be considered. A non-reflective boundary condition of convolutional PML (perfectly matched layer) type is applied at the top of the domain (Komatitsch and Martin, 2007). As illustrated in Fig. 2, it uses a sponge zone over which flow variables are attenuated increasingly with altitude, up to the top boundary where perturbations are zero. This damping is implemented by applying the following coefficient:

where $dt$ is the time step of the simulation, $zPML$ is the altitude in the beginning of the PML zone, and $dPML$ its size. The attenuation is therefore zero at $z=zPML$ and it reaches a maximum of $1/dt$ at the top boundary of the computational domain. Note that since the moving frame is attached to the aircraft, it advances at supersonic speed. Therefore, the acoustic fluctuations naturally leave the computational domain at its left boundary, without the need of a non-reflective condition. The sonic boom wave is injected through the right-hand side of the domain by interpolation. Linear interpolation is used, which proves sufficient, especially when the original sonic boom wave is much more discretized than the mesh. Usually, only the acoustic pressure waveform $p\u2032$ is available as an input, for instance, from ray-tracing simulations, which are typically used to propagate the wave from the aircraft near-field through the atmosphere and towards the ground. Thus, the relation between pressure and the other variables must be specified. Assuming isentropic perturbations without a mean flow (Thisse *et al.*, 2015), the boundary conditions on the right of the domain are written as

where $c0=\gamma p0/\rho 0$ is the sound speed and the wave is injected at an angle *θ* from the horizontal, equal to $\theta =sin\u22121(1/M)$ for a homogeneous atmosphere at rest. For sufficiently small values of $p\u2032/p0$, these relations become the classical plane wave relations between the acoustic pressure and the acoustic density and velocity. In particular, for the input waveforms considered thereafter, the differences between the nonlinear relations in Eqs. (8)–(10) and the plane wave relations were found to be negligible. It is, however, recalled that the full Euler equations are solved with this method, rather than the linearized Euler equations to account for possible nonlinear effects during the propagation of the sonic boom, for example, in the case of focusing.

The flow variables are initialized with their ambient values. The sonic boom wave progressively develops from the right boundary and, after a transient period, extends over the entire moving frame. To reduce computation time, the simulations are performed in two steps. First, a simulation is run for a flat ground and is stopped once the transient period is over and a steady sonic boom wave is obtained. Second, the simulation with the irregular terrain is launched using the flow at the last iteration of the transient computation as the initial condition. This allows one to perform simulations for different ground profiles without having to compute the transient period each time. Note this is only the case if the grid size remains unchanged.

The code is parallelized using the OpenMP parallel-computing language (Mattson *et al.*, 2019) to reduce restitution times. The flow variables are obtained in the temporal domain, and the resulting signals are post-processed, notably to compute the perceived noise. To do so, two different metrics are used: PL (perceived level) (Stevens, 1972), which is well suited to various sonic boom signatures (Leatherwood *et al.*, 2002), and CSEL weighting, which gives emphasis to lower frequencies. Other metrics, namely, ASEL, BSEL, DSEL, ESEL, and ISBAP, were shown to be highly correlated to human perception (Loubeau *et al.*, 2015), like PL. For conciseness, they are not discussed hereafter. However, Appendix B shows these metrics vary with ground elevation in a similar way as PL and CSEL.

### B. Case description

#### 1. Input parameters

The configurations investigated are described in this section. Two sonic boom waves are considered. The first is an N-wave, which is the classical sonic boom wave expected from a conventional supersonic aircraft. The chosen wave has a rise time of 0.0011 s. The second is the low-boom C25D wave, that originates from the near-field signature of a notional configuration aiming at representing a sonic boom demonstrator class vehicle. It was in particular used in the 2nd AIAA Sonic Boom Workshop (AIAA, 2017; Rallabhandi and Loubeau, 2019). This near-field signature was then propagated down to the ground by ONERA using BANGV nonlinear ray tracing code in a stratified and absorbing atmosphere (Loubeau and Coulouvrat, 2009). The rise time of this low-boom wave is equal to 0.014 s which is much larger than that of the N-wave. Both waves are presented in Fig. 3(a). They are injected through the right boundary of the domain, as illustrated in Fig. 3(b) which shows the reflection of the N-wave in the baseline case of a flat terrain through a pressure fluctuation map.

Throughout this study, the atmosphere is kept homogeneous and at rest following the conditions in Table I, which allows one to focus on topographic effects. Note the code allows the use of varying atmospheric conditions, which will be taken into account in future studies. Furthermore, the flight Mach number is set to 1.6.

ρ_{0}
. | c_{0}
. | γ
. | M
. |
---|---|---|---|

1.22 kg/m^{3} | 340 m/s | 1.4 | 1.6 |

ρ_{0}
. | c_{0}
. | γ
. | M
. |
---|---|---|---|

1.22 kg/m^{3} | 340 m/s | 1.4 | 1.6 |

Simple ground profiles are investigated for easier identification of mechanisms associated with topographic effects on sonic boom reflection. A single terrain depression, a hill and a sinusoidal ground profile are chosen, and discussed in Secs. III and IV. A maximal elevation variation *H = *50 m is chosen, associated with a terrain irregularity characteristic length *L = *300 m, leading to a large slope and allowing to highlight topographic effects in these simplified configurations.

#### 2. Computational domain parameters

For all simulations, the moving frame is shifted along the *ξ*-direction by one grid size every two time steps. The CFL number defined by $CFL=c0dt/d\xi $ is thus equal to $CFL=1/(2M)$. This yields a CFL number of 0.3125, well below 1 throughout the domain in curvilinear coordinates.

A mesh convergence study is performed to choose a grid size leading to accurate sonic boom propagation. Table II gives the evolution of peak pressure and perceived noise with grid size for both N and C25D sonic boom waves reflected over a flat surface. Compared to the wave initially injected into the computational domain, the peak pressure variation appears negligible for all grid sizes, with only 1.1% variation at worst with the N-wave and a grid size $d\xi =0.25$ m. Perceived noise levels prove more sensitive to grid size, and variations heavily depend on the type of metric considered. PL and CSEL metrics are chosen for this study, as explained in Sec. II A. Indeed, with the N-wave and the coarsest grid ($d\xi =0.25$ m), the error compared to the initial wave reaches −4.0 PLdB and only −0.7 dBC. This is due to the nature of these metrics, as CSEL weighting takes into account the entire range of audible frequencies, including low frequencies, while PL gives much more weight to high frequencies than to low ones. These variations are greatly reduced with a grid size of $d\xi =0.1$ m, although a difference of –0.5 PLdB persists, becoming negligible at $d\xi =0.06$ m. The C25D wave has a lower frequency content, in particular, with a longer rise time than the sharp N-wave. This is reflected on the evolution of perceived noise with grid size, with a –1.3 PLdB difference compared to the initial C25D wave with the coarse grid and a negligible error with $d\xi =0.1$ m, while the difference is negligible even with the coarsest grid using CSEL.

Wave . | dξ (m)
. | $ppeak$ . | PL (dB) . | CSEL (dB) . |
---|---|---|---|---|

N | initial | 47.8 Pa | 98.2 | 103.2 |

N | 0.25 | −1.1 % | −4.0 | −0.7 |

N | 0.10 | −0.05 % | −0.5 | −0.1 |

N | 0.06 | +0.01 % | −0.01 | +0.02 |

C25D | initial | 38.8 Pa | 83.2 | 97.2 |

C25D | 0.25 | −0.008 % | −1.3 | −0.05 |

C25D | 0.10 | −0.003 % | −0.3 | −0.002 |

Wave . | dξ (m)
. | $ppeak$ . | PL (dB) . | CSEL (dB) . |
---|---|---|---|---|

N | initial | 47.8 Pa | 98.2 | 103.2 |

N | 0.25 | −1.1 % | −4.0 | −0.7 |

N | 0.10 | −0.05 % | −0.5 | −0.1 |

N | 0.06 | +0.01 % | −0.01 | +0.02 |

C25D | initial | 38.8 Pa | 83.2 | 97.2 |

C25D | 0.25 | −0.008 % | −1.3 | −0.05 |

C25D | 0.10 | −0.003 % | −0.3 | −0.002 |

The spectra of these waves are represented in Figs. 4(a) and 4(d). They show the $d\xi =0.1$ m mesh captures high frequencies well compared to the initial wave, as well as to the $d\xi =0.06$ m mesh in the case of the N-wave. Also, note the vast improvement compared to the coarsest mesh. This confirms the good convergence with a grid size of $d\xi =0.1$ m. The other plots presented in Fig. 4 show similar spectra, this time for reflection over irregular terrain: at *x* = −40 and 60 m, corresponding to positions in the downward and upward slopes of a terrain depression as illustrated in Fig. 2. Note the initial wave corresponds to a flat surface case and it is shown as reference only. Concentrating on the N-wave first, Figs. 4(b) and 4(c) display slightly degraded spectrum with $d\xi =0.1$ m compared to $d\xi =0.06$ m. This is due to the curvilinear transformation of the mesh, which, in the case of a terrain depression, leads to slightly larger grid sizes in the vicinity of the ground irregularity than with the Cartesian mesh. In the case of the C25D wave, the $d\xi =0.1$ m spectra frequency content appears satisfactory at both positions, as shown in Figs. 4(e) and 4(f). Note the initial wave does not account for topographic effects which explains spectrum shape differences. In view of this result and taking into account computational cost, a grid size of $d\xi =0.1$ m is used for the C25D and sinusoidal cases. However, considering the sensitivity of PL noise level to high frequency content, a grid size of $d\xi =0.06$ m is used for the N-wave terrain depression and hill cases in order to ensure good accuracy. See Appendix A for a discussion on the sensitivity of PL with the N-wave.

In the case of local terrain irregularities (terrain depression, hill), the moving frame domain size is set to 800 × 400 m, while with sinusoidal ground profiles which affect a larger region, it is 1000 × 400 m. The moving frame mesh sizes are thus for the local terrain irregularities of 32 × 10^{6} and 90 × 10^{6} points for grid sizes of $d\xi =0.1$ m and $d\xi =0.6$ m, respectively, and for the sinusoidal ground profile of 40 × 10^{6} points.

The simulations were led with 32 processors. For the cases with locally irregular terrain, computational times $d\xi =0.1$ m are about 1000 and 730 CPU hours for the transient period and the useful signal simulation over 1.6 km respectively, and 4400 and 3400 CPU hours $d\xi =0.6$ m, while they are of 1200 and 2200 CPU hours in the case of a sinusoidal ground profile to obtain a signal over 4 km.

## III. IMPACT OF ISOLATED TERRAIN IRREGULARITIES

The objective of this section is to characterize the sonic boom reflection mechanisms induced by topographic variations and their effect on the noise perceived at ground level in the case of simple local elevation variations, in order to isolate and explain topographic effects.

### A. Terrain depression

A terrain depression is considered first. It is defined by the relation

where *H = *50 m and *L = *300 m. The ground profile is visible on the pressure fluctuation maps shown in Fig. 5, which are represented at different iterations and for both the N and C25D waves. Their evolution along the full profile is available for the N-wave in Mm. 1. Focusing on the N-wave, first note the pressure map is heavily affected locally where the wave is reflected on the ground. As the ground profile slope reduces, expansion occurs, affecting the reflected sonic boom, in particular by the reduction of its peak pressure. This is visible in the latter part of the terrain irregularity in Fig. 5(b). The opposite is true in the case of an increasing slope such as in Fig. 5(a), where compression occurs and the peak pressure increases. Second, a strong variation of the reflected wave can be observed, both in angle and amplitude, which is directly related to the change in elevation at ground level. Third, a focal zone is present in Fig. 5(c), leading to the most significant change in the pressure fluctuation field, with an additional high amplitude zone at the focal point, but also a component of the wave propagating back down towards the ground. This phenomenon can be explained by analyzing the wavefronts represented at different instants in time in Fig. 6(a). They are computed using a ray-tracing method (Candel, 1977; Gainville, 2008; Scott *et al.*, 2017). They are folded in the upward part of the terrain depression, at *x > *0, which leads to the focal zone observed on the pressure map. The caustic curve (red) corresponding to this focal zone is also plotted. In this terrain depression case, it is reflected on the upward part of the profile at *x = *53 m. In addition, a ray is plotted (green). After reflection (*x* = –100 m), it grazes the top of the terrain irregularity near *x* = 150 m. Subsequent rays will be reflected on the upward part of the depression, so that this ray delimits a shadow zone. Indeed, reflected rays cannot reach the zone below this ray, to the right of the terrain irregularity at ground level, for *x > *150 m. Figure 6(b) shows the pressure fluctuation map presented in Fig. 5(c), on which the wavefronts and caustic curve are superposed. Note the pressure map and the wavefront corresponding to the same instant (full black line) are in good agreement. At this macroscopic level, the same mechanisms are found in the case of the low-boom C25D wave in Figs. 5(d)–5(f). Indeed, on these pressure maps, the local variations at ground level, the evolution of the reflected wave, and the focalization phenomenon are globally comparable with the N or C25D sonic boom waves.

Now, let us focus on how these phenomena affect the signals obtained at ground level. Waveforms simulated at different axial positions and resulting from the reflection of the N-wave are presented in Fig. 7. They are represented following a relative time $\tau =t\u2212tflat$, where $tflat$ is the time at the beginning of the waveform in the corresponding flat case. The signal that is the most affected compared to the smooth surface case was computed at position *x = *70 m, just after the folding of the wavefront. This results in a U-wave, which is characteristic of focal zones and which originates from the additional contribution propagating back towards the ground. The resulting additional peaks are of larger amplitude and sharper than the first one. At position *x = *150 m, which corresponds to the end of the terrain depression, the U-wave is still present, but it has already been heavily attenuated. The rapidity of this attenuation can be explained by the presence of a shadow zone resulting from the elevation variation, as explained in the previous paragraph. Rays cannot reach this region and the energy of the signal is reduced. Finally, the signal at *x* = −40 m, before wavefront folding occurs, is much less affected, although a low amplitude secondary bump can be observed for $\tau =0.25$ s. In addition, the N-wave peak appears more rounded at this position. In the same way as for the pressure fluctuation fields, the waveforms obtained with the C25D wave follow a similar trend as with the N-wave, although there are some differences. Indeed, the secondary peaks caused by wavefront folding have lower amplitudes than the first one. They also remain much more rounded with C25D than with the N-wave, corresponding to a lower frequency content of the original wave.

The signals discussed in the previous paragraphs are used to estimate the noise perceived at ground level with PL and CSEL metrics. It is normalized by the noise levels obtained with a flat terrain in each case, in order to highlight the influence of ground elevation. Its evolution over the ground profile is given in Fig. 8(a) in the N-wave case, along with the terrain slope $dh/dx$. In the first part of the terrain depression ($x=[\u2212150;0]$ m), both PL and CSEL variations due to topography remain negligible. After the minimum elevation point of the terrain depression, where the slope is zero and continues to rise ($x=[0;75]$ m), both PL and CSEL are significantly amplified by ground elevation variations, by 5.5 PLdB and 4.9 dBC, respectively. This corresponds to the focal zone and the U-wave resulting from wavefront folding, as described previously. As the elevation gradient reduces again between the upward part of the terrain depression and the subsequent flat ground ($x=[75;150]$ m), the perceived noise collapses back towards levels obtained over flat ground. This corresponds to the attenuation of the U-wave illustrated in Fig. 7(a).

A similar plot is shown for the C25D wave in Fig. 8(b). PL noise variations due to topography remain negligible in the downward part of the terrain depression, while CSEL variations are present, which is consistent with the lower frequency content of the C25D wave. Noise amplification reaches 1.1 dBC where the slope rises, towards the bottom of the terrain depression ($x=[\u221275;0]$ m). This can be explained by the rise in signal energy caused by the superposition of the forming U-wave with the main waveform. This leads to an increase in wave amplitude in particular, as shown in Fig. 9 at position *x* = −22 m. The peak-to-peak amplitude reaches 72.3 Pa, that is 20% more than the initial C25D wave. The CSEL amplification due to topography then sharply reduces at the position *x = *24 m, where the forming U-wave is superposed with the minimum of the original wave, cancelling out the negative peak. Once the U-wave has fully separated from the primary peak, noise spikes up before reducing back to the flat case levels at the end of the terrain depression, in a similar way as with the N-wave. The maximum perceived noise amplification reaches 4.6 PLdB and 3.5 dBC. This is lower than with the N-wave, by 0.9 PLdB and 1.4 dBC, respectively.

### B. Hill

This section is dedicated to sonic boom reflection over a hill. Similar to the ground profile discussed in the Sec. III A, it is defined by the relation

where *H = *50 m, *L = *300 m. The ground profile is visible on the fluctuating pressure field maps presented in Fig. 10. They show both N and C25D waves are affected locally at the point of reflection on the ground and that their reflected component is strongly impacted, in a similar way as with the terrain depression. The evolution of the N-wave pressure fluctuation maps along the full profile is available in Mm. 2. Figures 10(c) and 10(f) show two focal zones are generated this time: the first is due to the transition between flat ground and the upward slope ($x=[\u2212150;\u221275]$ m) and it results in a low amplitude additional contribution at ground level; the second, generated by the downhill part ($x=[75;150]$ m) leads to stronger levels on the ground. Figure 11 depicts the same pressure map as Fig. 10(c), in the case of the N-wave, on which wavefronts and caustic curves are superposed. It shows wavefront folding over the upward part of the hill leads to the first focal zone, which moves upwards with time, while the second is generated as the wave moves from the downward part of the hill to the flat surface. This time, the focal zone moves towards the right and the caustic curve remains close to the ground. This explains why the additional contribution at ground level is of larger amplitude for the second focal zone, as the decay of a signal away from a caustic is related to the distance from the caustic [see, e.g., Salomons (1998)].

As for the terrain depression case, a U-wave is observed in Fig. 12 where focal zones are present, for both the C25D and N-waves. It is of low amplitude at *x* = –30 m, a position in the first half of the ground profile, where the additional contribution at ground level due to wavefront folding is weak. It is more significant at *x *= 270 m, where this contribution is stronger, but the secondary peak amplitude remains much smaller than in the terrain depression case (see Fig. 7), where the reflection of the caustic leads to increased amplitude. However, with the hill, this secondary peak is barely attenuated at positions further along the ground profile, and the signals at *x = *340 m display similar secondary peak amplitudes as at *x = *270 m. This is explained first by the absence of a shadow zone, which accelerates the decay in the case of the depression, and second by the proximity of the caustic curve with the ground after the hill, because of the relation between signal attenuation and distance from the caustic. As previously, a small secondary bump is observed at *x = *100 m, just before wavefront folding occurs, and the main peak is blunter than the original N-wave.

The perceived noise levels computed from these waveforms, obtained at ground level in the N and C25D cases and normalized by the corresponding flat case, are plotted in Fig. 13, along with the elevation gradient. Variations with topography remain small up to the second and stronger wavefront folding, towards the end of the hill (*x > *150 m). Then, as detailed in Sec. III A for the C25D wave, as the U-wave is formed and superposed to the original wave, CSEL increases before decreasing sharply as this superposition cancels out part of the signal. Finally, CSEL rises again as the U-wave detaches from the original wave. These CSEL variations are larger with C25D than with the N-wave, due to the C25D wave's thicker peak. In both cases, PL is much less affected by the superposition of the U-wave with the original wave, but it reaches similar amplification once the U-wave detaches, leading to additional peaks in the waveform. Unlike with the terrain depression, this focal zone leads to a plateau in the perceived noise plot in all cases. This is due to the slow attenuation of the U-wave along the caustic line, as explained in the previous paragraph. The maximum amplification due to topography along these plateaus are very close with the N and C25D waves, reaching 1.2 PLdB and 1.3 PLdB, respectively, and 0.9 dBC for both waves. These are significantly lower than in the case of a terrain depression (Fig. 8), which can be explained by the smaller slope variation between the hill's downward part and the flat terrain compared to the terrain depression.

## IV. IMPACT OF SINUSOIDAL TERRAIN

While Sec. III was devoted to isolated terrain features, this section deals with sonic boom propagation above a sinusoidal terrain that can be considered as a succession of hills or, equivalently, of terrain depressions of the same size. The ground profile is given by

with the same values of *H* and *L* as those used in Sec. III.

An example of a fluctuating pressure field map after propagation over more than ten terrain periods is presented in Fig. 14 for the N-wave. First, it can be observed that the wavefront structure of the incident and reflected boom waves near *x = *3.6 km is similar to that for the terrain depression shown in Fig. 5(b). Second, after these two waves, a large number of wavefronts are observed. These are due to the generation of U-waves at each terrain irregularity and to their subsequent diffraction on the following irregularities.

Figure 15 shows a pressure waveform at the ground level in the case of both the N and C25D waves. In accordance with Fig. 14, it does not present a single contribution, as for a flat ground, or two main contributions as observed for isolated terrain irregularities, but many more. The two first contributions, considered separately, are very close to the waveform for the terrain depression at *x = *70 m in Fig. 7 and they are the dominant ones. The additional contributions are however far from negligible: thus, for the N-wave, the peak-to-peak amplitude is of 23 Pa for the third contribution around $\tau =0.6$ s, 15 Pa for the fourth one at $\tau =0.9$ s and 11 Pa for the fifth one at $\tau =1.2$ s. A similar order of magnitude is obtained for the C25D wave. Compared to the main U-wave, each successive additional contribution appears more and more rounded, and the sharp peaks of the U-wave are absent. Its duration seems to increase as well, which shows diffraction tends to cumulatively reduce the high-frequency content of the U-wave.

The evolution of the PL and CSEL metrics along the ground, relative to the flat ground case, is depicted in Fig. 16. It presents a periodic behavior, associated with the periodic variation of the ground profile. For comparison, the evolution of PL and CSEL along the isolated terrain depression (see Fig. 8) is also represented in Fig. 16. Interestingly, a close match is observed. This shows that the effects due the terrain depression prevail over those caused by the hill. This also indicates that, while having a large impact on the pressure waveforms, the cumulative effects related to the topography have little influence on the noise perceived, which is thus mostly governed by the local evolution of the ground profile, at least using the metrics employed in this work.

## V. CONCLUSION

In this paper, a numerical method for the prediction of sonic boom reflection over irregular terrain based on the Euler equations is presented. The physical mechanisms associated with varying elevation are characterized using academic ground profiles (terrain depression, hill, periodic). The effect of slope variation is observed locally in the reflection zone, on the reflected wave, and on the wavefronts. Furthermore, wavefront folding is observed where slope variations are largest, leading to focal zones, caustics, and U-waves, which strongly affect the signals and lead to the most significant variations in perceived noise at ground level. Similar mechanisms were observed with the classical N-wave and the low-boom C25D wave, but the different natures of these waves lead to differences in how these mechanisms affect the waveforms and the perceived noise at ground level, depending on the chosen metrics.

Between the terrain depression and the hill ground profiles, the first was found to have the strongest effect, in particular, with stronger focalization at ground level, resulting in higher amplitude U-waves and larger levels of perceived noise. This terrain depression focalization was also found to be dominant in the case of a sinusoidal ground profile, which yielded similar perceived noise levels. Furthermore, the sinusoidal ground profile allowed us to highlight the presence of cumulative effects that affect the signal tail and could prove significant with real terrain.

In order to estimate the significance of the mechanisms highlighted in the present study, the investigation of topographic effects on sonic boom reflection using real ground profiles is under way. In addition, while this study focuses on topography, the tools used allow us to consider various atmospheric conditions, so that they can be used to investigate the combined impact of topographic and meteorological effects on sonic boom propagation in the future. This study has highlighted the influence of the topography based on a two-dimensional model. While similar effects are expected, the three-dimensionality of both the ground surface and of the sonic boom could affect the quantitative results that remain to be evaluated.

## ACKNOWLEDGMENTS

This project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 769896 (RUMBLE, 2020). This publication reflects only the authors' view and the Innovation and Networks Executive Agency (INEA) is not responsible for any use that may be made of the information it contains. It was performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program “Investissements d'Avenir” (ANR-16-IDEX-0005) operated by the French National Research Agency (ANR). This work was granted access to the HPC resources of PMCS2I (Pôle de Modélisation et de Calcul en Sciences de l'Ingénieur et de l'Information) of Ecole Centrale de Lyon, PSMN (Pôle Scientifique de Modélisation Numérique) of ENS de Lyon and P2CHPD (Pôle de Calcul Hautes Performances Dédiés) of Université Lyon I, members of FLMSN (Fédération Lyonnaise de Modélisation et Sciences Numériques), partner of EQUIPEX EQUIP@MESO.

### APPENDIX A

This appendix aims at completing the discussion on mesh convergence in Sec. II B 2 including the effect of elevation variations in the case of the N-wave reflecting over the terrain depression ground profile. Furthermore, this appendix primarily concentrates on PL variations with grid size, as this metric is particularly sensitive in the presence of the N-wave. This is due to PL's sensitivity to high frequencies, present in the N-wave due to its sharp peak and short rise time.

Figure 17 shows waveforms obtained using three different mesh sizes at positions *x* = −40 and 70 m along the terrain depression, i.e., on either side of the point of minimum elevation at *x = *0 m. At *x* = −40 m, the three waveforms are very close, although one may notice peaks are slightly more rounded with the coarse grid size $d\xi =0.25$ m. The waveform obtained at *x = *70 m is the least favourable case, with particularly sharp secondary peaks. These secondary peaks are somewhat affected by the grid size. Using grid sizes of 0.25, 0.1, and 0.06 m, the second peak reaches 49.6, 64.3, and 72.2 Pa, respectively, and the third peak reaches 39.0, 50.5, and 57.9 Pa.

Overall, these waveforms are thus very close, especially with the two finer grid sizes. This explains why CSEL, which accounts for a wide range of frequencies, yields accurate results with both these meshes, as detailed in Sec. II B 2. Furthermore, the perceived noise variations along the terrain depression are plotted for different grid sizes using CSEL in Fig. 18(b), which shows levels obtained with these two meshes remain close even in the presence of ground elevation variations. On the contrary, PL is more sensitive to high frequencies than to low frequencies and consequently the convergence of the mesh with an N-wave is not so simple with this metric, as illustrated in Fig. 18(a). First notice the substantial improvement in accuracy between a grid size of 0.25 and 0.1 m, compared to $d\xi =0.06$ m. The difference between the two finer meshes remains small over flat ground ($x<\u2212150$ m and *x > *150 m), as detailed in Table II, but some significant deviations appear with elevation variation, especially between *x* = −122 and 35 m. This corresponds to the downhill part of the ground profile, where as the curvilinear mesh adapts to the surface, some cells are larger than the nominal size. On the other hand, the uphill section, where wavefront folding occurs, is favourable in terms of mesh convergence, and some cells are smaller than the nominal size in this part of the mesh.

In view of this, a grid size of 0.06 m is used to simulate the propagation of the N-wave over both the terrain depression and hill cases. Our confidence in the accuracy of the results with this mesh is supported by the good accordance with the reference initial wave in the flat case (see Table II), as well as by the spectra in Fig. 4 which show this grid size allows to capture frequencies in the higher frequency range, including with non-flat ground. The downside of using the finer mesh is increased computational cost. While it remains reasonable over short ground profiles, it becomes a limiting factor when longer profiles or a large number of cases are considered. A grid size of 0.1 m is therefore chosen for all cases with the C25D wave, which has a lower frequency content than the N-wave (see Sec. II B 2). Furthermore, the discussion in Appendix B shows the major variations in perceived noise induced by topographic effects are captured by each of the metric types which have been tested, including CSEL. A grid size of 0.1 m will therefore also be used for the sinusoidal terrain case, as well as other long ground profiles which will be studied in the future. Then, CSEL is used to calculate the perceived noise in the case of the N-wave, while both PL and CSEL are used with the C25D wave.

### APPENDIX B

The PL and CSEL metrics have been used in Secs. III and IV to evaluate the impact of terrain on perceived noise of sonic booms. This appendix presents corresponding results for five other metrics, that were shown to correlate highly with human perception (Loubeau *et al.*, 2015). These metrics are ASEL, BSEL, DSEL, ESEL, and ISBAP.

Figure 19 shows the evolution of the seven metrics along the terrain depression ground profile, normalized by their values for a flat terrain. Concerning the N-wave in Fig. 19(a), ASEL, BSEL, DSEL, and ESEL have similar values and are closely related to PL. As explained in Sec. III A, CSEL follows a comparable trend, but its maximum is smaller by about 1 dB. For the C25D wave in Fig. 19(b), notice the curves for CSEL and DSEL are almost superimposed and that ASEL, BSEL and ESEL have a similar evolution to PL. For both boom waves, the ISBAP metric, which is composed of PL and CSEL, is located in between the PL and CSEL curves. Analogous results (not shown for conciseness) are obtained for the hill ground profile.

In conclusion, the seven metrics follow the same evolution along the two ground profiles. It is thus expected that the variations related to terrain elevation determined for one of these seven metrics would also be representative of those of the other ones.^{1}

^{1}

See supplementary material at https://www.scitation.org/doi/suppl/10.1121/10.0003816 for videos of the reflection of the low-boom C25D wave along both the terrain depression and hill ground profiles using pressure fluctuation maps ($[\u221250;50]$ Pa).

## References

*InterNoise 2015*, 9–12 August,