Propagation of bird vocalizations in the Alpine environment
Didier DRAGNA(1) ∗ , Loïc BERGER(1) , Sébastien OLLIVIER(2) , Frédéric SEBE(3)(1) Univ Lyon, Ecole Centrale de Lyon, INSA Lyon, Université Claude Bernard Lyon I, CNRS, LMFA, UMR 5509,

36 Avenue Guy de Collongue, F-69134, ECULLY, France(2) Univ Lyon, Université Claude Bernard Lyon I, Ecole Centrale de Lyon, INSA Lyon, CNRS, LMFA, UMR 5509,

43 boulevard du 11 novembre 1918, F-69100, VILLEURBANNE, France(3) Univ Lyon, Université Jean Monnet, CNRS, ENES/Neuro-PSI, UMR 9197,

23 rue Michelon, F-42023, SAINT-ETIENNE, France

Abstract
Rock ptarmigan is a bird living in the mountains tops, whose population tends to decrease in Europe over the last
decade, notably due to climate change, habitat loss and human disturbance. Censusing of number of birds in the
Alps is realized via point-count monitoring using birds vocalization. However, detection of bird vocalizations
by the human hear is somewhat tricky in these non-homogeneous landforms, due to sound propagation effects.
In order to improve the counting method, it is necessary to understand how bird vocalizations propagate in the
Alpine environment. With this aim, a propagation model based on a ray-tracing approach is developed andincludes topography, ground reﬂection and absorption and atmospheric refraction and absorption. The model isﬁrst validated in both frequency and time domains against reference results obtained with a numerical solution of
the linearized Euler equations for simple cases. Propagation of the bird vocalization in the Alpine environmentis then illustrated. It is shown that the active space of the bird dramatically depends on the environment. The
modeling of active and detection space will provide precise information for the birds monitoring, enabling
managers to deﬁne more precisely their counting plans.
Keywords: Ethology, Outdoor propagation, Ray tracing

1

INTRODUCTION

The rock ptarmigan is a bird (see Fig. 1 (a)) living in the northern parts of Eurasia and North America. In
France, it can be found in subalpine and alpine habitats of the Pyrenees and the Alps at altitude greater than1800 m [1]. Its characteristic vocalization consists of a train of pulses, with frequency content between 900 Hz
and 3700 Hz see (Fig. 1 (b)-(c)). Due to climate change and habitat transformation, its population is decreasing.
Rock ptarmigan is often considered as a bioindicator, meaning that its population is a global indicator of thebiodiversity in alpine environments.
While acoustic monitoring starts to be used [2], the current counting method is based on a human-based approach. The counting area is covered by several observers, spaced apart each other by 500 m, which corresponds
to an empirical distance at which rock ptarmigan can not be heard anymore. This protocol has several biases;
a single individual can be counted multiple times or can not be detected at all. Indeed, propagation of thebird vocalization in Alpine environment is not simple to due to a complex landform (topography, ground effect)
and meteorological conditions. In order to improve current counting methods, it is therefore required to better
understand and predict propagation in this environment.
The objective of this paper is to develop a propagation model, suited for the Alpine environment. The paper isorganized as follows. In Sec. 2, possible approaches to predict sound propagation are examined and motivations
for the use of geometrical methods are discussed. The ray-tracing approach is described. Test-cases to validate∗ corresponding author, didier.dragna@ec-lyon.fr

3810

the numerical model are performed. Application of the model for propagation of rock ptarmigan vocalization ina mountain environment is illustrated in Sec. 3.

1.5

60

1

50

0.5

40

0

30

-0.5

20

-1

10

-1.5
0

0.05

0.1

0.15

0.2

0.25

0.3

0
10 2

10 3

Figure 1. (top) A pair of rock ptarmigans in spring. [3] (bottom) Vocalization of the rock ptarmigan: (left)
waveforms and (right) spectrum.

2

RAY-TRACING APPROACH

2.1 MOTIVATION
They are several ways to predict sound propagation outdoors. Analytical approaches are the fastest ones toobtain reliable results. Standards to compute noise maps are thus based on these approaches. A lot of commercial softwares are available for this purpose. They often take as input data GIS (geographic information
system) ﬁles, which allow one to account easily for complex landforms. Typical calculation times for an areaof several square kilometers is only of a few minutes. There are however two main drawbacks for the use of
these softwares in a bioacoustic context. First, the standards implemented in these softwares are developed foroctave or one-third octave frequency bands, which does not make it possible to obtain the bird vocalization after
propagation. Second, the formulas employed in these standards are often derived from empirical considerations,especially for ground reﬂection and refraction by the atmosphere. In addition, they were mostly developed forurban environment and their validity for mountain environment is not clear. Note that here has been recently
a sound propagation model [4] for anthropogenic noise in wildland environment that has been implemented in
the ArcGIS platform. Numerical approaches are also numerous. Wave-based methods are appealing as they
can account accurately for most of the phenomena encountered in atmospheric acoustics. Their computational
cost however remains prohibitive for large 3D domains. Geometrical methods are an interesting compromise
between accuracy and efﬁciency. They are based on a high-frequency assumption, or more precisely on theassumption that the characteristic length scale of variations of the ambient variables is large compared to thewavelength. They can handle specular reﬂection on an irregular terrain and refraction by inhomogeneities of theatmosphere. Accounting for diffraction or diffuse reﬂection is possible, but is usually not included because it
signiﬁcantly complicates the approach. There are however considerable limits to ray-tracing. There are thus twoparticular regions in which predictions given by the ray-tracing approach are expected to be erroneous. First
ones are shadow zones, which are the regions in which no ray penetrates and are due to atmospheric refraction or topography. Geometrical approaches would thus predict a null amplitude for the acoustic ﬁeld. This is3811

of importance in Alpine environment, for which terrain irregularities can induce large shadow zones near the
ground. The other critical region are caustics, which are points or surfaces on which the cross-sectional area ofthe ray tube tends to zero. On caustics, the ray-tracing approach predicts an inﬁnite amplitude.
2.2 DESCRIPTION
20

6

0

4

-20

2

-40

0

-60

-2

-80

-4

-100
0

500

1000

1500

10 -3

-6
0.02

0.025

0.03

0.035

0.04

Figure 2. (top) Energy spectral densities and (bottom) time-domain signatures for a source at (0, 0, 1 m) and a
receiver at (300 m, 0, 3 m) and for a homogeneous atmosphere.
In the ray-tracing approach [5], the acoustic pressure is sought as p̂( f , x) = A( f , x)eiωτ(x) , with A the amplitudeof the pressure which is assumed to be slowly varying and τ the travel time. In an inhomogeneous atmosphereat rest with sound speed proﬁle c̄(x) and density proﬁle ρ̄(x), the ray equations are given by:dx= c̄ n(1)dtdn= (n · ∇c̄)n − ∇c̄(2)dtwith n the unit vector tangent to the ray. From an acoustic point of view, a bird can be considered as apoint source, generating a spherical wavefront. The rays are therefore parametrized with two angles ϕ and θ ,corresponding to the elevation and azimutal angles, respectively. Introducing the geodesic elements rϕ = dx/dϕand rθ = dx/dθ as well as the variables qϕ = dn/dϕ and qθ = dn/dθ , it can be shown from Eqs. (1) and (2)that they satisfy the equations:drϕ= (rϕ · ∇c̄)n + c̄qϕ
(3)
dt
drθ
= (rθ · ∇c̄)n + c̄qθ
(4)
dt
dqϕ
= (n · ∇c̄)qϕ + (qϕ · ∇c̄)n − vϕ + (n · vϕ )n
(5)
dt
dqθ
= (n · ∇c̄)qθ + (qθ · ∇c̄)n − vθ + (n · vθ )n(6)
dt
with vϕ = (rϕ · ∇)∇c̄ and vθ = (rθ · ∇)∇c̄. The cross-sectional area of the elementary ray tube along the ray
can then be obtained with S(t) = |rϕ × rθ |. The amplitude of the pressure along a ray is determined with the
Blokhintsev invariant, which states that A2 (t)S(t)/(ρ̄(t)c̄(t)) is constant along a ray tube. The constant is usuallyobtained by comparing the solution of the ray equations with that of the Helmholtz equation for a homogeneous
unbounded medium [6]. For a point source, the amplitude of the pressure is A = ŝ( f )/r, with r the distanceand ŝ( f ) is the source strength. Finally, the amplitude of the pressure along a ray is given by:
cos ϕs ρ̄(t)c̄(t) 1/2(7)
A(t) = ŝ( f )S(t) ρ̄s c̄s3812

where ϕs is the elevation angle at launch and ρ̄s and c̄s are the density and sound speed at the source. The
parameter ŝ( f ) can also be interpreted as the acoustic pressure at 1 m from the source and can account for a
possible source directivity. The pressure at an observer is therefore the sum of the contribution of all the rays
that reach the observer, which are referred to as eigenrays. Accounting for phase shifts at caustics and possible
ground reﬂections, it can be written as:

1/2
iωti
iωti cos ϕs, i ρ̄(ti )c̄(ti )(−i)n(i) ∏ Ri j
(8)p̂( f ) = ∑ Ai (ti )e = ∑ ŝ( f )eS(ti )
ρ̄s c̄sjiiwhere n(i) is the number of times that the eigenray i has crossed a caustic and Ri j is the reﬂection coefﬁcientassociated to the j-th reﬂection of the eigenray i.
(a) ac = 0.1 s-1

8(b) ac = 0.5 s-1
20

6

15

4

10

2

5

0

0

50

100

150

200

250

300

0
0

20

20

0

0

-20

-20

-40

-40

-60

-60

-80

-80
-100

-100
0

1

500

1000

1500

10 -3

50

0

1.5

100

150

500

200

250

1000

300

1500

10 -3

1

0.5

0.5

0

0
-0.5

-0.5
-1

-1
0.875

0.88

0.885

0.89

-1.5

0.87

0.88

0.89

0.9

Figure 3. (top) Eigenrays, (middle) energy spectral densities and (bottom) time-domain signatures for a sourceat (0, 0, 1 m) and a receiver at (300 m, 0, 0.3 m). The sound speed gradient is set to (left) ac = 0.1 s-1 and(right) ac = 0.5 s-1 .
The computation of the sound pressure with the ray-tracing approach is performed in two steps. First, theeigenrays at each receiver are determined. Second, the amplitude along the eigenrays is computed and the3813

pressure at the receiver is determined with Eq. (8). The ﬁrst step is by far the more computationally intensive.
In particular, the search of eigenrays is somewhat tricky, especially in a large 3D domain. In this study, we
use a brute force approach. The launch angles are discretized on a grid ϕ ∈ [ϕ min ; ϕ max ] and θ ∈ [0; 2π],with corresponding grid size Δϕ and Δθ . For each set of launch angles (ϕ, θ ), the ray trajectory is obtainedby integrating Eqs. (1) and (2), with a Runge-Kutta algorithm. The minimal distance between the ray andthe receiver dmin (ϕ, θ ) is then determined. After having computed all the ray trajectories, the eigenrays aredetected by computing the minima of the function of two variables dmin and by checking that for each local
minimum, the minimum distance to the ray is smaller than a threshold, set to 1 m. Once the eigenrays havebeen obtained, Eqs. (3)-(6) for the geodesic elements are integrated and the amplitude of the pressure along
the ray is determined with Eq. (7). The pressure in the frequency domain at the receiver is ﬁnally obtained bysumming coherently the contributions of the eigenrays with Eq. (8). The pressure in the time domain can thenbe computed by calculating the inverse Fourier transform.
The angular steps Δϕ and Δθ must be sufﬁciently small. Indeed, at long range, the launch angles of eigenrays can be very close, in particular in a downward-refracting atmosphere (see, e.g., [7]). Therefore, a large
number of rays has to be launched to determine all the eigenrays. For the simulations performed in Sec. 3,
approximately 107 rays were launched. The code is written in MATLAB and is run on a standard workstation.
2.3 VALIDATION
In order to validate the ray-tracing approach, comparisons with a numerical solution of the linearized Euler
equations (LEE) by a ﬁnite-difference time-domain approach [8] are performed. For this, we consider soundpropagation above a ﬂat rigid ground in an atmosphere whose sound speed proﬁle varies as: c̄ = c0 + ac x3 , withc0 = 340 m s-1 . The density is constant with ρ̄ = ρ0 = 1.22 kg m-3 . The source is located at (0, 0, 1 m) andgenerates a short pulse.
In the ﬁrst test-case, the atmosphere is homogeneous with ac = 0. The rays are straight lines and there are
only two rays, corresponding to the direct and reﬂected contributions. Fig. 2 shows the energy spectral densities(ESD) obtained with the ray-tracing approach from Eq. (8) and the corresponding time-domain signals calculatedby an inverse Fourier transform. Analytical solution and the numerical solution of the LEE are also shown forcomparison. It can be seen that the three solutions perfectly agree in both frequency and time-domain.
3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0
0

20

40

60

80

100

0
0

50

100

150

200

Figure 4. Eigenrays for a source at (0, 0, 1 m) and a receiver at [left] (100 m, 0, 3 m) and [right] (200 m, 0,
3 m). The sound speed gradient is set to ac = −0.1 s-1 .
A downward-refracting atmosphere with ac > 0 is now considered. The eigenrays for a receiver at (300 m,
0, 3 m), as well as the ESD and the time-domain signal at the receiver are shown in Fig. 3 for sound speed
gradients equal to ac = 0.1 s and ac = 0.5 s. For ac = 0.1 s, an additional eigenray is observed, compared toa homogeneous atmosphere. The ESD is slightly overestimated with reference to the LEE solution by approximately 3 dB over the whole frequency range. A single contribution is seen on the time-domain signal, as thearrival times of the eigenrays are very close. For ac = 0.5 s, there are 9 eigenrays detected by the algorithm.

3814

(a) ac = −0.1 s-120(b) ac = −0.5 s-1

0

0

-20

-20

-40

-40
-60

-60

-80

-80

-100

-100
0
4

500

1000

1500

0

10 -4
10

500

1000

1500

10 -5

2
5

0
0

-2
-4
0.58

0.585

0.59

0.595

0.6

-5
0.58

0.585

0.59

0.595

0.6

Figure 5. (top) Energy spectral densities and (bottom) time-domain signatures for a source at (0, 0, 1 m) and a
receiver at (200 m, 0, 3 m). The sound speed gradient is set to (left) ac = −0.1 s-1 and (right) ac = −0.5 s-1 .
The ESD obtained with the ray-tracing is similar to that of the LEE solution below 600 Hz and above 900 Hz
despite an overestimation of a few dB. The destructive interference pattern at a frequency of 800 Hz is also
shifted towards a lower frequency. A good agreement is also observed in the time-domain signals. In particular,
the three ﬁrst contributions before t < 0.88 s perfectly superimpose. The amplitude of the last contribution withan arrive rime of t = 0.88 s is largely overpredicted. Analysis of the raytube cross-sectional area (not shownhere) reveals that the two eigenrays with the longest arrival time are close to a caustic. The evaluation of theiramplitude is therefore inaccurate.
In the last case, an upward-refracting atmosphere with ac < 0 is examined. Eigenrays for receivers at (100 m,
0, 3 m) and (200 m, 0, 3 m) are plotted in Fig. 4. Close to the source, there are still two eigenrays. At larger
distances, a single eigenray corresponding the the direct contribution remains, while at even larger distances
there is no eigenray at all. Fig. 5 shows the ESD and the time signals predicted at (200 m, 0, 3 m) forsound speed gradients of -0.1 s-1 and -0.5 s-1 . For ac = −0.1 s-1 , the ray-tracing approach gives only one
contribution (see Fig. 4). The ESD and the time signal are in good agreement with those of the LEE solution.
For ac = −0.5 s-1 , the receiver is in the shadow zone. Geometrical acoustics then leads to a null pressure.
Actually, compared to the case ac = −0.1 s-1 , the amplitude of the pressure signal is reduced by a factor 4,
corresponding to a 12 dB reduction.

3

PROPAGATION IN THE ALPINE ENVIRONMENT

Results obtained with the ray-tracing approach on the Flaine area are illustrated. Fig. 6 shows the topographyof the Flaine area that has been obtained from GIS tools. It is observed that on this area of a few square
kilometers, the difference in altitude is about 800 m. A valley can also be seen at x1 = 3.5 km. It can therefore
be expected that sound propagation would be totally different depending on the source location.

3815

(a)

(b)

2400
2000
1600
0
1000
2000
3000
4000
0

1000

2000

3000

4000

(c)

(d)

1400

12
10

1200

1400

70
60

1200

50

8

1000

6

40

1000

30
4

800

2

600

800

20
10

600

0

0

2000

2200

2400

2600

2800

2000

2200

2400

2600

2800

Figure 6. (a) Relief map of the Flaine area and (b) zoom on the are of interest. (c) Number of eigenrays
received at each receiver on the map. (d) Map of the sound presure level (SPL).
In this illustration, the source is located at x1 = 1.2 km and x2 = 2.3 km and its height is 1 m above theground. Its location is represented by a dot in Fig. 6. The calculation domain is limited to the blue rectangle
in Fig. 6 (a). A zoom on the relief map restricted to the calculation domain is shown in Fig. 6 (b). The
terrain slope is largely varying. Near the source, it is almost constant, while for x1 ≤ 1 km the slope is more
pronounced. At top right corner of the domain, the terrain is almost ﬂat. Finally, near the source, three small
elongated hills oriented along the x1 -axis can be distinguished. Two of them are located 50 m on the left ofthe source in Fig. 6 (b). The third one is located at x2 = 2.6 km on the right of the source.
Receivers are placed on the computational domain at a height of 1.5 m and are spaced regularly, separated by
30 m. The atmosphere is inhomogeneous. The sound speed is varying linearly with the altitude. The soundspeed gradient is 0.2 s−1 and the reference temperature is 5 ◦ C at x3 = 2350 m. An impedance model of a
freshly fallen snow has been used for ground reﬂection. Fig 6 (c) indicates the number of eigenrays at each
receiver. It is signiﬁcant at the left and at the top of the source, with a maximum number equal to 12. Thereis no eigenray behind the two small hills at the left of the source. In addition, the third hill prevents ray
propagation to the right of the source. The decrease of the slope induces a shadow zone at the receivers locatedat the top of the domain. Similarly, an other shadow zone is generated due to the increase of the slope forthe receivers at the bottom of the domain. There are therefore only few receivers, for which the ray-tracing
algorithm gives information on the sound pressure. The corresponding SPL map is shown in Fig 6 (d). It gives
information on what is referred to as the active space which is the area around the bird where its vocalization
can be detected. It is seen that the topography has a large inﬂuence on the SPL. In particular, the SPL does
not depend only on the distance from the source. The higher SPL is obtained in the region at the left and at
the top of the source for which the number of eigenrays is the more important.

3816

4

CONCLUSIONS

Numerical prediction of bird vocalization propagation in the Alpine environment has been investigated. For
this, a ray-tracing approach was chosen. While there are inherent limitations to geometrical acoustics (shadow
zones and caustics), the model developed accounts for most of the physical phenomena encountered in soundpropagation in a mountain environment. Outputs are narrow-band spectra as well-as waveforms, which could
allow the analysis of bird vocalizations after propagation. Test-cases above a ﬂat rigid ground have been done
to validate the model. Finally, application of the model to predict the active space of rock ptarmigan has been
illustrated and the feasibility of the approach has been demonstrated.
In future works, further analysis of the active spaces of rock ptarmigan will be done. In particular, inﬂuence ofthe meteorological conditions on the active spaces will be studied. The inverse problem that is determining fora given receiver the area on which rock ptarmigan can be detected will be investigated. Moreover, the study
will be extended to other bird species.

ACKNOWLEDGEMENTS
This work was performed within the framework of the Labex CeLyA of the Université de Lyon, within the programme ‘Investissements d’Avenir’ (ANR-10-LABX-0060/ANR-16-IDEX-0005) operated by the French National
Research Agency (ANR).

REFERENCES[1] Marin-Cudraz, T.; Muffat-Joly, B.; Novoa, C.; Aubry, P.; Desmet, J.-F., Mahamoud-Issa, M., Nicolè, F.; Van
Niekerk, M.H.; Mathevon, N.; Sèbe, F. Acoustic monitoring of rock ptarmigan: A multi-year comparison
with point-count protocol, Ecological Indicators, Vol 101, 2019, pp 710-719.
[2] Marin-Cudraz, T.; Muffat-Joly, B.; Novoa, C.; Mathevon, N.; Sèbe, F. The use of bioacoustics to improverock ptarmigan (lagopus muta) spring counts in the French Alps and Pyrenees, Grouse News, Vol 53, 2019,
pp 23-28.
[3] Rock ptarmigan.
5/29/2019).

Wikipedia,

https://en.wikipedia.org/wiki/Rock_ptarmigan(Last

viewed[4] Reed, S.E.; Boggs, J.L.; Mann, J.P. A GIS tool for modeling anthropogenic noise propagation in natural
ecosystems, Environmental Modelling & Software, Vol 37, 2012, pp 1-5.
[5] Candel, S. Numerical solution of conservation equations arising in linear wave theory: application to aeroacoustics, Journal of Fluid Mechanics, Vol 83 (5), 1977, pp 465-493.
[6] Jensen, F.B.; Kuperman, W.A.; Porter, M.B.; Schmidt, H. Computational ocean acoustics, Springer, second
edition, 2011.
[7] Dragna, D.; Cotté, B.; Blanc-Benon, P.; Poisson, F. Time-domain simulations of outdoor sound propagation
with suitable impedance boundary conditions, AIAA Journal, Vol 49 (7), 2011, pp 1420-1428.
[8] Dragna, D.; Blanc-Benon, P.; Poisson, F. Time-domain solver in curvilinear coordinates for outdoor soundpropagation over complex terrain, Journal of the Acoutical Society of America, Vol 113 (6), 2013, pp 37513763.

3817