Applying the diffusion equation to urban scenarios: Computational analysisof the diffusion coefﬁcientRaúl PAGÁN MUÑOZ (1) , Juan Miguel NAVARRO RUIZ (2) , Maarten HORNIKX(3)(1) Eindhoven University of Technology, The Netherlands, r.pagan.munoz@tue.nl(2) Universidad Católica de Murcia, Spain, jmnavarro@ucam.edu(3) Eindhoven University of Technology, The Netherlands, M.C.J.Hornikx@tue.nl

ABSTRACT
Prediction of sound propagation is playing a key role in the planning and design process of urban areas and during
the last decades different techniques have been developed for the computation of the sound ﬁelds in cities. Amongthem, the diffusion equation, based on the propagation of sound energy, is a simple and attractive tool for certain
scenarios where the diffuse ﬁeld is predominant, such as inner city environments. The diffusion equation is a wellknown and efﬁcient method to compute sound ﬁelds in rooms with a low amount of absorption. However, thismethod, which includes the diffusion coefﬁcient describing the diffusivity of the energy propagation, requires further
understanding for urban environments. In this work, the diffusion coefﬁcient has been evaluated using wave-based
numerical techniques for simpliﬁed two-dimensional urban scenarios. The parameter is computed spatially from theacoustic intensity and sound energy density ﬁelds in every direction for different geometries. This paper investigates
the directional components of the diffusion coefﬁcient.
Keywords: Diffusion equation, urban sound propagation

1

INTRODUCTION

The diffusion equation describes sound propagation by a diffusion process of the sound energy. The method relies
mainly on a single parameter - the diffusion coefﬁcient D - that requires further understanding for urban environments.
The hypothesis assumed in this paper is that the sound ﬁeld can be approximated by a diffusion process in canyon-type
street conﬁgurations. For these cases, the diffusivity of sound energy has different properties in different directions
and as a consequence, the diffusion coefﬁcient is not uniform (1). Although some authors have derived analyticalexpressions or methodologies to estimate the coefﬁcients for urban cases, for instance (1) and (2), there is still alack of knowledge regarding how the diffuse sound energy propagates in these cases. The main contribution ofthis paper is to investigate the directional components of the diffusion coefﬁcient by using a novel methodology tospatially compute D using solutions of a wave-based method for simpliﬁed urban environments. Section 2 presents
the principles of the diffusion equation for sound propagation with special emphasis on the presentation of the nonhomogeneous and directional diffusion coefﬁcient for urban cases. In this work, the computation of the coefﬁcientsis done by using the sound pressure and acoustic velocity as calculated by the discontinuous Galerkin (DG) wavebased time domain method; the approach is also valid for any numerical method that allows the steady state acoustic
intensity and sound energy density to be computed at discrete positions of the domain under investigation. The
full methodology of obtaining the diffusion coefﬁcients is detailed in Section 3. The approach has been used in
two simpliﬁed two-dimensional scenarios, described in Section 4: a straight, long street and an intersection of two
perpendicular long streets. The solutions of the diffusion coefﬁcients for these scenarios are presented in Section 5
and the ﬁnal conclusions can be found in Section 6.

3756

2

THE ACOUSTIC DIFFUSION EQUATION



−2 is
Fick’s ﬁrst law, adapted to a diffuse sound energy ﬁeld, is presented in Equation 1, where w(r,t) kg m−1s
the sound energy density. This law states that the ﬂow of acoustic energy (sound intensity vector J kg s−3 ) occurs
from regions of high energy to areas of low energy concentration, with a proportionality relation between the acoustic
intensity and the sound energy density gradient (∇w). The proportionality factor in the equation is the diffusion
coefﬁcient D [m2 s−1 ]. The derivation of Fick’s second law - or more popularly the diffusion equation (DE) - fordiffuse sound ﬁelds was presented by (4) (see Equation 2). The absorption by the domain surfaces is included in σ[s−1 ], which represents the probability rate of a sound particle to be absorbed during one second. The DE presented
in Equation 2 represents sound energy density propagation in the time domain. In order to obtain frequency solutions,the equation needs to be resolved for each frequency band of interest separately.
J(r,t) = −D∇w(r,t).

(1)∂w(r,t) = D∇2 w(r,t) − σ w(r,t).
∂t(2)

The diffusion coefﬁcient is deﬁned for room acoustics problems by Valeau et al. in (3) as a term of the diffusion
equation that takes into account the room morphology through the mean free path λ [m]. This classical deﬁnition
corresponds with the theoretical derivations presented in (4) for rooms of proportionate dimensions, where D takes
a constant value depending on the geometrical properties of the space as D = λ c0 /3, where c0 is the sound velocity
and λ = 4V /S, with V and S the volume and surface area of the proportionate room. However, when the space ismore complex, as for instance in a long room or in a canyon-type street, the diffusion coefﬁcient is not as simple.
Visentin et al. (5) carried out an investigation to study the validity of Fick’s law of diffusion in room acoustics froma numerical evaluation of the gradient Equation 1. One of the main conclusions of the study was that D = λ c0 /3,
is only valid in areas close to the source in long rooms, while linearly increasing with distance from the source.
Therefore, they proved that the sound ﬁeld in long rooms is described by non-homogeneous diffusion, i.e., with aspatially varying diffusion coefﬁcient. In this type of space, where one dimension is much bigger than the others,
the mean free path is not homogeneous and spatially depends on the location of the source and receiver. In receiverareas close to the source in an inﬁnite long room, the distance between successive reﬂections is shorter (limited by
the distance between parallel surfaces when the sound waves travel as a standing wave in the corridor), than in areas
far from the excitation, where the sound waves have on average a larger angle with respect to the source (assuming
zero degrees angle is perpendicular to the side walls). The hypothesis assumed in the present paper is that the soundﬁeld can be approximated by a diffusion process in canyon-type streets conﬁgurations by using the DE model aspresented in Equation 2 without fundamental modiﬁcations. In this type of street, the diffusivity of sound energy
has different properties in different directions and, as a consequence, the diffusion coefﬁcient is not uniform anddepends on the direction of propagation (1). The spatially dependent directional diffusion coefﬁcients are presented
in Equation 3 for three-dimensional rectangular streets in Cartesian coordinates. The diffusion gradient Equation2, assuming the diffusion coefﬁcients D(r) given in 3, is employed to calculate the spatially dependent directional
diffusion coefﬁcients of urban cases.
⎛

Dx (r)
D(r) = ⎝ 0
0

⎞
0
0
0 ⎠.
Dy (r)0
Dz (r)

(3)

3757

3

METHOD TO OBTAIN THE DIFFUSION COEFFICIENT FOR URBAN APPLICATIONSIn this Section, the methodology for obtaining spatially dependent directional diffusion coefﬁcients D(r) is presented.
The approach begins with the frequency domain J( f ) and ∇w( f )1 , calculated from solutions of a full wave-based
method at a large number of equally-spaced discrete locations in the physical domain under investigation. The method
is detailed for two-dimensional rectangular cases. The approach is based on solving Equation 1 by evaluating the ratio
between the solutions of the steady state components of J and the corresponding energy density gradients ∇w.
One assumption in the derivation of the diffusion Equation 2 is that the variations per mean free path of w, J, andthus D, must be small. However, the solutions of a full wave-based method, as proposed here, implicitly include
the interference effects due to the wave nature of sound propagation and therefore, the results are not as smooth as
the diffusion equation would demand. In order to reduce this effect, calculating the energy density and the intensity
in broader frequency bands, as indicated in Section 3.1, contributes to smoothing the results, as well as the local
spatial averaging, performed according to Section 3.2. After the calculation and local averaging, the sound energy
density gradient is computed as presented in Section 3.3. In Section 3.4, the gradient and the acoustic intensity
are then averaged over the narrow dimension of the street. Before obtaining the diffusion coefﬁcients, a smoothing
process over the acoustics intensity and sound energy density is applied according to Section 3.5. Finally, the diffusion
coefﬁcients are obtained according to Section 3.6. Additionally, the variations per mean free path of the energy density
and the acoustic intensity are not expected to be smooth along the narrow dimension of a street canyon. Therefore, themethodology presented in this section is only employed for computing the diffusion coefﬁcients along the directionof the biggest dimension of the street.
3.1 Broad band frequency calculation
The energy density and the active part of the acoustic intensity, indicated generically in Equation 4 by Ψ(r, f ), arecalculated for a broad frequency band fbb by a summation of the results from the lowest frequency of interest fl to
the highest fu , as indicated in Equation 4.
fu

Ψ(r, fbb ) = ∑ Ψ(r, f ).

(4)

f = fl

3.2 Receivers and local spatial averaging process
A large number of discrete locations r in the spatial domain is needed in order to obtain a meaningful representationof the sound ﬁeld to get the diffusion coefﬁcients. The total number of receivers in the whole domain is Nr and the
distance between receivers is denoted as Δr [m] for every Cartesian dimension. Δr is chosen to be constant everywherein the domain for the ease of the averaging process and the later calculation of the spatial energy gradient. The spatial
average of the broadband w and J are computed at a number of equally spaced recording positions in each Cartesian
direction nr (same number in each direction). They are arithmetically averaged and the result is assigned to the central
position. The total number of averaged positions is nr = n2r . The total number of evaluation positions Nr is reduced to
Nr = Nr /nr , and the distance between positions is now Δr. The averaged frequency solutions are referred as Jx (r),Jy (r) and w(r) for a certain frequency range.
3.3 Energy density gradient calculation
The next step is to compute the energy density gradient over the whole domain by using central differences for interior
data points, while using single-sided differences along the edges of the domain. The following notation is used forthe spatial discretisation of the energy density in the two-dimensional domain: wi, j = w(iΔr, jΔr). The gradient
1 Unless otherwise indicated, the reference to the frequency domain in the solutions has been omitted for subsequent references to these quantities, i.e. J( f ) = J and ∇w( f ) = ∇w.

3758

calculation in the x-direction is described in the following Equations 5 for interior data points and 6 for the positions
located at the edges (equivalent expressions are used for the y-direction), where Nx is the total number of averaged
points in the x-direction:∇x wi, j ≈∇x w1, j ≈w2, j − w1, j,
Δr

wi+1, j − wi−1, j,
2Δr∇x wNx , j ≈wNx , j − wNx −1, j
.
Δr

(5)

(6)

3.4 Averaging from 2D to 1D solutions
Since the streets considered in this approach are rather narrow, a second spatial averaging process is performed over
the results Jx (r), Jy (r), ∇x w(r) and ∇y w(r). Essentially, the two-dimensional solutions are transformed into a set ofone-dimensional results for each Cartesian dimension by spatially averaging the results over the other direction. In
this approach, it is assumed that the solutions are constant along the narrower dimension of the street. When one of
the results is expressed, for instance, by Jx (rx ), it means that the values have been averaged and only spatial values
along the x-dimension remain, i.e. the values in the y-direction have been averaged.
3.5 Smoothing process over the sound energy density and the acoustic intensity
Despite the spatial averaging process and the calculation of solutions in a broad frequency band, the results stillpresent strong variations within the spatial domain. In order to reduce these local variations, a moving averaged
process is applied to the solutions obtained in the previous Section 3.4, that is equivalent to low-pass ﬁltering the data.
The response of the smoothing process is given by the Equation 7 for a generic solution ψx (rx ). The ﬁlter smooths
the data by replacing each data point with the average of the neighbouring data points deﬁned within a certain spandeﬁned as nspan = 2Nsm + 1.

ψx,sm (i) =1[ψx (i + Nsm ) + ψx (i + Nsm − 1) + · · · + ψx (i − Nsm )],(2Nsm + 1)

(7)

where ψx,sm (i) is the smoothed value for the ith data point and Nsm is the number of data points at each side of point i.
The span must be an odd integer number. Moreover, the end points of the function are not smoothed. Finally, the span
is adjusted for data points that cannot accommodate the speciﬁed number of neighbours on either side by reducing
the number of points.
3.6 Calculation of the diffusion coefﬁcient
The diffusion coefﬁcients are ﬁnally computed at the averaged positions, for the frequency range of interest, using
Equation 8 for the x-dimension and an equivalent expression for the y-dimension.
Dx (rx ) = −4Jx (rx ),∇xW (rx )(8)WAVE-BASED URBAN SCENARIOSTwo simpliﬁed urban problems have been investigated in this paper: a straight, long street (LS) and the perpendicularintersection of two long streets (LX). The two-dimensional domains are shown in Figure 1a. The façades of the
scenarios are built using quadratic residue diffusers (QRDs) in order to obtain high values of the surface diffusion(dΨ ) and scattering (s) coefﬁcients. The arrangement of the QRDs followed the indications in (7) and the design3759

has been done for frequencies higher than 500 Hz. The geometrical details of the QRD are given in Figure 1b.
The DG scenarios have been designed for a frequency range of interest from 500 to 2500 Hz. The QRDs represent
façade irregularities with a constant low sound absorption coefﬁcient, while the beginning and the end of the streets
are simulated as completely absorbent openings. The calculations in DG are initiated with a broadband pressure
distribution and zero acoustic velocity.

0.66m
0.06m
0.03m

0.09m

0.16m

0.12m

(b)

(a)

Figure 1. a) Dimensions of the two-dimensional urban domains including the source position (black dot) and theopenings (red lines): upper domain straight, long street (LS), lower domain perpendicular intersection of two longstreets (LX). b) Geometrical detail of the two-dimensional QRD used at the boundaries of the scenarios to simulate
the façade’s irregularities.
The results of the absorption coefﬁcient for random sound incidence αr , and the scattering and the diffusion coefﬁcients (averaged value for all source positions dΨ ) according to (7) for the used QRD diffuser are presented in Figure
2. These values were obtained from a numerical evaluation using the DG model.
1

0.8

0.8

0.8

0.6

0.6

0.6
r

s [-]

[-]

1

d [-]

1

0.4

0.4

0.4

0.2

0.2

0.2

0

0

500

1000

1500

2000

2500

0

0

500

1000

1500

2000

2500

0

0

500

1000

1500

Frequency [Hz]Frequency Hz

Frequecy [Hz]

(a)

(b)

(c)

2000

2500

Figure 2. Results of the QRD used in this study for 1/1 (black big circles) and 1/3 (red small circles) octave bands: a)diffusion coefﬁcient dΨ , b) scattering coefﬁcient s, and c) random absorption coefﬁcient αr .

3760

4.1 Removing the direct sound and specular reﬂections
The nodal discontinuous Galerkin (DG) numerical time-domain wave-based method (6) was employed in this work
for solving the linearised Euler equations representing sound propagation in urban scenarios. The impulse responses(IRs) of the acoustic solutions of DG are recorded at equally spaced multiple locations r in the domain. These
solutions represent the steady state condition when transformed to the frequency domain. The source and receivers
locations in the two-dimensional coordinate system are indicated as rs = [xs , ys ] and rnr = [xnr , ynr ], respectively. The
time step is computed depending on the maximum frequency of interest fmax according to Δt < 1/(2.5 fmax ) [s]. The
end of the time signals is tapered by a single-sided Gaussian window with a length ≈ 6 ms to avoid Gibbs effects. As
by deﬁnition, the diffusion equation does not include the coherent part of the sound ﬁeld, the beginning of the IRs istapered by windowing the ﬁrst part of the IRs. The window length is constructed by computing the number of time
steps (n1a ) the direct sound takes to arrive to every receiver location in the case of receivers located in the line-of-sightfrom the source. The cases where the receiver location is not in the line-of sight from the source, the number of time
steps until the ﬁrst arrival (n1a ) has been computed from the distance between the source and the closest diffracted
corner, plus the distance from that corner to the receiver position. The windowing function Γw is constructed using
Equation 9, where the total number of recorded time steps is denoted by Nt with indexes nt = (1, 2, . . . , Nt ) and αw = 7
and βw = 4 are constants controlling the shape of the exponential function. The length of the exponential window Nw
is ﬁxed in this work (Nw Δt ≈ 2.4 ms) to remove the direct sound and assuming that all reﬂections are non-coherent
for frequencies above 500 Hz, according to the design of the QRD.
⎧⎪
⎨ 0
nt −Nw 2βw
Γw =e−αw log(10) Nw
⎪
⎩
1

5

for 1 ≤ nt ≤ n1a ,for n1a < nt ≤ n1a + Nw ,for n1a + Nw < nt ≤ Nt .

(9)

RESULTS OF THE WB SCENARIOS AND CALCULATION OF D

In this section, the results of the acoustic intensity, the energy density and the gradient of the energy density arepresented as calculated from the solutions of the WB method. Additionally, the calculated diffusion coefﬁcients areshown for the computed cases, LS and LX. The set of graphs shown in Figure 3 presents the results of w, Jx , ∇x w and
Dx of both scenarios. Excluded from this initial analysis the areas in the proximity of the source and street openings.
The one-dimensional solutions of the LS case in Figure 3 show an increase of the diffusion coefﬁcient along the street,
while the energy density resembles an exponential decay as the distance from the source increases. This is in linewith the diffusion coefﬁcients reported by Visentin (5) for long rooms. As a consequence of the decay of the energy
density, the gradient smoothly approaches zero as the slope of the energy decay decreases (see Figures 3a and 3c).
The intensity follows a similar exponential pattern as the energy density, smoothly decaying with distance from thesource as shown in 3b. Close to the right end of the main streets, there is a predominant direction of the intensity, as
shown in Figure 4. However, in this area there is a drop of the value of the diffusion coefﬁcients that can be explained
due to the fact that the back-scattering in nearby areas of the openings is highly reduced. Moreover, the intensity ﬁeld
around the source location is not presenting a predominant direction, as shown in Figure 4, that can be additionallyobserved by the change of sign of Jx in Figure 3b. In the same area, the energy density presents its maximum value,
decaying towards both openings (see the area around the source location in Figure 3a). Overall, it can be concluded
that the openings of the streets affect the results in their proximity, questioning the validity of the methodology forcomputing the diffusion coefﬁcients in these areas.
In the LX scenario, the energy density is distributed along the main street, but part of the energy is distributed to theside streets as shown in Figure 4b. The ﬂow of energy occurs from the source towards the ends of the main and sidestreets. In this scenario, the results of the acoustic variables (w, J, ∇w and D) are presented in Figures 3 and 5, forthe main and side streets, respectively. Each graph includes red dotted vertical lines indicating the position of thecrossing street. Of particular interest in the main horizontal street is the analysis of the intensity ﬁeld direction aroundthe crossing area that can be observed in Figure 4b. The energy density presented in Figure 3a for this scenario is

3761

1.5

2
1

7000
6000

wx [Kg m-2 s-2]

3

10-12

1

1

Jx [Kg s-3]

w [Kg m-1 s-2]

4

10-10

0.5
0
-0.5

0.5

Dx [m2 s-1]

10-12

0
-0.5

-1

0
0

20

4000
3000
2000
1000

-1.5
0

40

5000

20

x [m]

-1
0

40

20

x [m]

(a)

0
0

40

x [m]

(b)

20

40

x [m]

(c)

(d)

y [m]

Figure 3. Results of the LS scenario (grey solid line) and the main street of the LX scenario (black broken line). a)w, b) Jx , c) ∇x w and d) Dx . The red dotted vertical lines represent the position of the crossing street (only for the LX
scenario).

3
0
0

5

10

15

20

25

30

35

40

25

30

35

40

x [m]

(a)
6

y [m]

4
2
0
-2
-4
0

5

10

15

20

x [m]

(b)

Figure 4. Direction of the acoustic intensity ﬁeld of the street scenarios: a) LS and b) LX.
10-13

1.5

10-13

2000

1

Jy [Kg s-3]

3
2
1

wy [Kg m-2 s-2]

1
0.5
0
-0.5

1500
0.5

Dy [m2 s-1]

4

w [Kg m-1 s-2]

10-11

0
-0.5

1000
500
0
-500

-1
-1

0
-20

0

y [m]

(a)

20

-1.5
-20

0

y [m]

(b)

20

-20

0

y [m]

(c)

20

-1000
-20

0

20

y [m]

(d)

Figure 5. Results of the two-dimensional cross street of the intersection model (LX). a) w, b) Jx , c) ∇x w and d) Dx .

3762

smoothly decaying until the crossing area where a more signiﬁcant drop occurs, corresponding with the distributionof the energy towards the side streets. This can be observed as well in the results of Jx (see Figure 3b). Clearly, after
the crossing area, the energy density and the acoustic intensity decay at a slower rate than before the crossing. The
diffusion coefﬁcient results present a clear drop in the main street right before the crossing and a sudden increase
after the crossing. In the cross street, the analysis is equivalent. The energy density is distributed from the crossing
to the side streets. Figure 4b shows how the acoustic intensity propagates with opposite directions in each side street
and Figure 5b shows how the intensity distribution is approximately mirror-symmetric. The coefﬁcients present anincrease with distance when entering the side streets until the proximity of the openings, as shown in Figure 5d.

6

CONCLUSION

A novel methodology for the calculation of diffusion coefﬁcients for the diffusion equation has been proposed forurban applications, based on the solutions of a wave-based method. The approach is based on locally computing the
diffusion coefﬁcients as the proportionality factor between the acoustic intensity and the gradient of the sound energy
density given by Fick’s ﬁrst law of diffusion. The variations per mean free path of w, J, and thus D, must be small andtherefore, efforts of the proposed methodology are oriented towards reducing the ﬂuctuations of the solutions whenusing a full WB model. The presented approach has been employed for two-dimensional2 scenarios composed ofrectangular domains where one dimension is much bigger than the other. The method has been applied to a straight,
long street and to a more complex urban case, the intersection of two perpendicular long streets. The proposed
methodology to compute the diffusion coefﬁcients has been proven to be a valid approach in areas where the soundﬁeld is smoothly decaying in a predominant direction. However, the method fails when the gradient of the energy
density approaches zero. These cases are associated with areas of the domain where there is not a predominant andunique direction of the energy ﬂow as near the source, openings or the narrow direction of the streets. The diffusion
coefﬁcients are found to be non-homogeneous, equivalently to the ﬁndings for long room cases (5).

REFERENCES
1. Picaut J, Simon L, Hardy J. Sound ﬁeld modeling in streets with a diffusion equation. The Journal of the AcousticalSociety of America, Vol 106 (5), 1999, pp 2638–2645.
2. Le Pollès T, Picaut J, Colle S, Bérengier M, Bardos C. Sound-ﬁeld modeling in architectural acoustics by a transport theory: Application to street canyons. Physical Review E, Vol 72 (4), 2005, pp 1-17.
3. Valeau, V., Picaut, J., Hodgson, M. (2006). On the use of a diffusion equation for room-acoustic prediction. The
Journal of the Acoustical Society of America, 119(3), 1504–1513.
4. Picaut, J., Simon, L., Polack, J.-D. (1997). A Mathematical Model of Diffuse Sound Field Based on a Diffusion
Equation. Acta Acustica united with Acustica, 83(4), 614–621.
5. Visentin, C., Prodi, N., Valeau, V., Picaut, J. (2012). A numerical investigation of the Fick’s law of diffusion inroom acoustics. The Journal of the Acoustical Society of America, 132(5), 3180–3189.
6. Wang, H., Sihar, I., Pagán Muñoz, R., & Hornikx, M. (2019). Room acoustics modelling in the time-domain with
the nodal discontinuous Galerkin method. The Journal of the Acoustical Society of America, 145(4), 2650–2663.
7. Cox, T.J., Dalenback, B., Antonio, P., Embrechts, J.J., Jeon, J.Y., Mommertz, E., et al. (2006). A Tutorial on
Scattering and Diffusion Coefﬁcients for Room Acoustic Surfaces. Acta Acustica united with Acustica, 92, 1–15.

2 Due to the maximum frequency of interest and the dimensions of the domain, it is not possible to employ the proposed technique using athree-dimensional WB model. This issue remains as a challenge for future work.

3763