Numerical Modeling of Infrasound Energy Radiation by Debris Flow Events

Debris flows constitute a severe natural hazard in Alpine regions. Studies are performed to understand the event predictability and to identify early warning systems and procedures. These are based both on sensors deployed along the channels or on the amplitude of seismic and infrasound waves radiated by the flow and recorded far away. Despite being very promising, infrasound cannot be used to infer the source characteristics due to the lack of a physical model of the infrasound energy radiated by debris flows. Here the simulation of water flow along a simple channel is presented, experiencing the fall from a dam, performed within the open source simulation code OpenFOAM. The pressure perturbation within the atmosphere produced by the flow is extracted and the infrasound signature of the events as a function of the flow characteristics is defined. Numerical results suggest that infrasound is radiated immediately downstream of the dam with amplitude and period that scale with dam height and water level. Modeled infrasound waveform is interpreted as being produced mostly by waves at the water free surface developing downstream of the dam. Despite the effect of sediments is not considered in this first study and will be implemented in future investigations, numerical results obtained with this simple model are in general agreement with experimental results obtained from array analysis of infrasound data recorded at Illgraben, Switzerland. Results highlight how numerical modeling can provide critical information to define a source mechanism of infrasound energy radiation by debris-flow, that is required also to improve early warning systems.


Introduction
Debris flows are water/sediment mixtures typically mobilized in steep mountain torrents affecting highly erodible areas. Given the large volumes (up to [ 200; 000 m 3 ), high velocity (up to [ 10 m/s), high mass loading, with density exceeding 2000 kg/ m 3 and unpredictability, debris flows are among the most hazardous processes affecting mountain territories and constitute a severe natural hazard for downstream communities (Jakob 2005). This leads to the instrumentation of many debris flows prone valleys in Alpine environment, adopting different solutions for automatic alert and notification of flow occurrence, as well as the realization of expensive confinement structures (Huebl and Fiebiger 2005).
Models and studies focusing on seismic energy radiation from rivers as well as debris flows are quite well established [see, ''e.g''., Burtin et al. (2008), Gimbert et al. (2014), Kean et al. (2015), Lai et al. (2018), Tsai et al. (2012)]. Lai et al. (2018) developed a mechanistic physical model for seismic energy radiation associated to debris flows. Moreover, similarly to what predicted for rivers (Burtin et al. 2008;Gimbert et al. 2014;Tsai et al. 2012), they showed that recorded seismic frequency depends only on the source-to-receiver distance, while the amplitude of seismic ground shaking is proportional to boulder size and to flow velocity as well as to boulder snout extension. This model was confirmed by the direct seismic observations of the Montecito debris flow occured on 9th January 2018 (Lai et al. 2018), despite additional seismic sources that might be related to water flow interaction with channel bed and banks and standing waves (Schmandt et al. 2013).
Differently, studies focusing on the infrasound energy radiation are still very few and any self-consistent model has been proposed yet. Up to now most of the studies in the literature (Huebl et al. 2013;Schimmel et al. 2018;Kogelnig et al. 2014;Li et al. 2012) focused primarily on event detectability for monitoring purposes. Very recently, Marchetti et al. (2019) performed array analysis of infrasound observations at Illgraben (CH). Based on signal characteristics, they suggested that infrasound is likely radiated by multiple sources (or by an elongated source) acting simultaneously along the channel, resulting into an incoherent signal with high amplitude. This is in general agreement with Schmandt et al. (2013), who studied seismic and infrasound energy radiation at Colorado river and suggested that infrasound energy was most likely radiated by waves at the water free surface. However, Marchetti et al. (2019) also showed that infrasound is most efficiently radiated by the flow interaction with dams and topography changes, in agreement with experimental observation at rivers (Kudo 1993;Feng. et al. 2014). The latter in particular, following the analytical solution developed by Ostrovsky and Bedard (2002) for a falling object hitting the water and inducing in the water surface a circular wavefront, modelled infrasound as being generated by non-stationary water surface variations within dam absorption pools that are acting as a dipole source.
Considering debris flows at Illgraben with variable volumes (25,000-70,000 m 3 ) and flow velocities (4-7 m/s), Marchetti et al. (2019) and Belli et al. (2021) showed how spectral amplitude and peak frequency of recorded infrasound seem to scale with flow volume.
Due to the lack of a physical model of the infrasound energy radiated by moving fluid, infrasound observations cannot be used to infer informations about the source characteristics. Therefore, while increasing number of experimental observations have been performed in the last decade and a large dataset of events is now available, numerical physical models of the infrasound energy radiation by debris flows and rivers are required to link acoustic waves to the initial physical parameters of the phenomenon, that is, density, velocity, terrain's characteristics.
In this paper the simulation of fluid flow along a channel is presented, experiencing the fall from a dam. We consider a two dimensional domain. The fluid phase is constituted by water and the geometry is very simplified. The solver (compressibleIn-terFoam) within the open source simulation code OpenFOAM (Open Field Operation And Manipulation) is used. A two phase mixture of water and air, with a water flow entering with constant velocity the numerical domain filled with air and flowing across a dam, is considered. Despite debris flows exhibit a rich variety of flow regimes, ranging from a solid-like behavior characterized by small deformations and frictional dissipation mechanisms to a gas-like behavior with strong collisions among the grains, the presence of sediment is neglected in this first study, to investigate the basic relation between flow parameters and infrasound wave characteristics. Their effects will need to be considered in future investigations in order to properly address debris flows. The Navier-Stokes equations are numerically solved and by extracting the pressure perturbation within the atmosphere the infrasound signature of the events is defined as a function of the flow characteristics. Numerical results are eventually compared with experimental measurements of the infrasound wavefield produced by debris flow events recorded in the spring of 2017 at Illgraben in Switzerland (Marchetti et al. 2019).

The Mathematical and Simulation Model
The acoustic waves emitted by debris flow falling down a vertical dam are studied. The dynamics of the phenomenon and the generation and propagation of acoustic waves in air is reproduced by a numerical model. The debris flow is described in a two dimensional simplified geometry that contains a single vertical step. The numerical solution is performed by using the compressibleInterFoam solver available in the open source CFD software OpenFOAM (https://www.openfoam.com/). It is a cell-centered, co-located finite volume method for multiphase flow. The Volume of Fluid (VOF) method is used as basis of the solver. The coupling between pressure and velocity in transient flows is performed by the Pressure Implicit with Splitting of Operators (PISO) algorithm (Issa 1986). Water and air are considered as a mixture of two compressible, non-isothermal, immiscible and not reacting phases. A phase-fraction is associated to each phase of the mixture. The function a denotes the fraction of water in the fluid and 1 À a the fraction of air. The regions of the domain where a ¼ 1 (a ¼ 0) correspond to the presence of water (air). The physical conditions are such that water and air can be considered as immiscible. In this case a should change abruptly from zero to one, the discontinuity indicating the surface separating the two phases. From a numerical point of view, the solver implements dedicated correctors able to capture the sharp interface between the two phases.
The density of the mixture is denoted by q ¼ aq w þ ð1 À aÞq a , where q w , q a are the density of water and air, respectively. To a high degree of approximation, the density of water and air can be modeled by the equations of state for a compressible perfect fluid (Martínez Ferrer et al. 2016) where p is the pressure field, T the temperature and R w , R a are the fluid constant and the gas constant, for water and air, respectively, and q w;0 is the water density at room temperature in the incompressible case.
The dynamics of the water-air mixture is governed by the following system of balance laws for mass, linear momentum and energy. The continuity equation for the mixture density q is where U is the mixture velocity vector. As indicated before, the phase fraction a is used in order to trace the surface separating the liquid and the gas phase of the mixture and obeys to the following mass balance equation where D Dt o ot þ U Á r is the material derivative. According to the Multidimensional Universal Limiter for Explicit Solution (MULES) procedure the standard balance equation is modified by the presence of the term U c , called compression velocity and equal to the relative velocity between water and air U c ¼ U w À U a . The compression velocity term ensures that the interface between the phases remains sharp and bounded between zero and one; for its numerical treatment we referee to Berberović et al. (2009). The evolution of the velocity of the twophase fluid is described by the following momentum equation Here, l eff ¼ 4ðal w þ ð1 À aÞl a Þ=3 is the mixture viscosity, l w and l a are the water and air viscosity, respectively, and e y is the unit vector along the ydirection. The last term of the equation describes the contribution due to the surface tension r between the two phases. It holds f r ¼ rkra; where k ¼ r Á ra jraj is the mean curvature of the interface (Brackbill et al. 1992).
The energy balance equation, assuming constant specific heats, is written with respect to the temperature, which explicitly enters the Eqs. (1)-(2), and reads as where E ¼ jUj 2 =2 is the specific kinetic energy and c V;w and c V;a the specific heat at constant volume for water and air, respectively. For any further details the reader is refereed to Moukalled et al. (2016) and Moukalled and Darwish (2000).

Simulation Results
An example that illustrates the simulations that have been performed is presented in the snapshots of Fig. 1. The simulation starts with the water entering the modeling grid filled with unperturbed air at t ¼ 0 s (Fig. 1a), thus simulating a flood wave. In order to simulate a constant slope of the terrain, a constant input velocity of the water is imposed, assumed to be 5 m/s. The depth of the fluid entering the simulation domain is equal to 0.5 m. This is resulting into a constant flux of 2.5 m 2 /s. The water flow falls down from the step, Fig. 1b), and impacts the terrain, Fig. 1c). After the dam, waves at the water surface develop downstream with a length scale between \ 1 to few meters.'' Eventually, the water reaches the end of the modeling grid, where water is allowed to exit, around at t ¼ 6 seconds, Fig. 1d). The flow evolution is controlled both by the constant inlet velocity and by the gravity, that cause a thinning of the water depth.
The dimensions of the domain are ð42 Â 100Þ m 2 , with ð12 Â 100Þ m 2 on the lest side of the step and ð30 Â 100Þ m 2 on its right, and they have been chosen in such a way that the most part of the reflection waves in air due to the boundaries disappears. The simulation domain is discretized with a uniform square grid for the lower part ð42 Â 5Þ m 2 , with Dx ¼ Dy ¼ 0:1 m , then Dy is graded linearly with a grading factor equal to 10; the grading factor being equal to the ratio between width of the final cell and the initial one. There are no appreciable variations on the results by further refining the grid.
The following initial conditions are taken: aj t¼0 ¼ 0, Uj t¼0 ¼ 0, pj t¼0 ¼ qgh þ p atm , where p atm is the atmospheric pressure at the top of the domain, g the constant gravity acceleration, h the height in the y direction, and Tj t¼0 ¼ 293 K . The inlet initial velocity U in is set equal to 5 m/s for all simulations. The total simulation time is equal to 100 s . For the temporal part, the mixed Euler-Cranck Nicolson scheme is used, with a blending factor of 0.5, and default upwind schemes for the space operators. The Courant number is equal to 0.3, as well as the interface Courant number. The physical parameters of the mixture used in the simulations are indicated in Tab. 1, where MW, H f , Pr are the molecular weight, the fusion point and the Prandtl number, respectively.
In order to investigate the atmospheric perturbation and the infrasound waves associated with the fluid flow, the value of the calculated atmospheric pressure is extracted within a point of the modeling grid located 20 m above the point where the water impacts downstream the dam; this impact points depend on the variables of the specific simulations, as the water depth, the dam height, the water inlet velocity. Spurious numerical effects were excluded by comparing the pressure extracted from different nodes of the computation grid, and by taking spatial average of the pressure. See Sect. 5.2 for more details.
The formal derivation of the correlation existing between the infrasound waves emitted by a moving fluid and the physical characteristics of the fluid (such as for example the composition, density, initial velocity and pressure) is still an open problem. In order to gain insight in this field, the modification of the pressure signals obtained by varying the height of the dam is analyzed. The snapshots of the evolution of the water flow at different times are showed in Fig. 2 where four cases have been simulated, each one corresponding to different dam height h (from up to bottom h ¼ 0, 1.0 m, 1.5 m, 2.0 m). Different snapshots are considered: at 2.0 seconds, when water is falling from the dam, at 6.0 seconds, when the flow is exiting from the modeling domain, and at 12.0 seconds, when stable features of the flow for most of the modeling tests are observed. These timings are kept constant for all the figures throughout the paper for consistency.
The corresponding pressure time series is depicted in Fig. 3a, b, where the red dashed vertical lines identify the timing of the snapshots shown in Fig. 2. Panel (a) shows the value of the absolute pressure extracted from the simulations in the absence of the dam (blue trace in Fig. 3) and with increasing dam height. In order to remove the long term fluctuation, related to the water inlet in the modeling grid, and appreciate the perturbation effects, in panel (b) the signal is plotted after application of a filter that removes the frequency components inferior to 2 Hz. At the very beginning, the simulations show low frequency oscillations (i in Fig. 3a) of the signal which are produced by the input of material in the modeling grid. This is followed, around 2 s in the pressure time series by a high frequency transient produced by the impact of water downstream from the dam (ii in Fig. 3). This is not synchronous for all cases and increases with the dam height, as more time is required for the water to first impact the ground.
The high frequency transient is consistent with the theoretical modeling and experimental analysis for infrasound radiated at waterfalls for rivers (Kudo 1993;Feng. et al. 2014) and debris-flows (Marchetti et al. 2019), that are identified as preferential sources of high energy signals. In particular, this was explained by Feng. et al. (2014) in terms of dipole like acoustic sources produced by waves induced in the absorption pool of the dam by the falling water. Immediately after the water impact (Fig. 3b), when the flow spreads downstream the dam, a lower frequency, almost monochromatic pressure oscillation starts to form, with emergent onset and cigar shape, reaches its maximum amplitude between 7 and 9 seconds (Fig. 3b) and eventually decreases. No roughness of the terrain is considered, except from the effect of the dam that causes the water flow to fall down and impact with the ground. While the modeled flow keeps a flat surface in the absence of a dam (Fig. 2a), waves at the water free surface develop downstream for the presence of the dam (Fig. 1b-d), with amplitude and wavelength that appear to scale with the dam height.
The pressure signals in frequency domain are plotted in Fig. 3c, d. In particular, the Power Spectral Density (PSD) of the pressure, once normalized for maximum amplitude, is depicted in linear scale in the panel c and for the sake of comparison, in the right panel, the same signals are displayed in a single plot in log scale. The first 2 seconds of each simulation (i in Fig. 3a) is neglected in order to eliminate the unphysical effects related to the initial condition on the pressure field which is chosen to be uniform on the entire domain. Spectral analysis of modeled waveforms shows that larger amplitude and lower frequencies are recorded with increasing dam height (Fig. 3c, d).
In Fig. 4 the effect of varying the volumetric fluxes of water entering in the domain is studied. In particular, simulations with variable water depth w (w ¼ 0.5, 1, 1.5 m) have been performed. The case with the height of the dam of 1 m and the input velocity of 5 m/s is discussed. The modeled pressure time series (Fig. 5), suggest that the water impacts the ground downstream the dam at around 2 seconds (ii in Fig. 5b), triggers low frequency oscillations afterwards, and exits the modeled grid between 4 and 5.5 seconds (iii in Fig. 5b). Such a time difference is a consequence of the different water depths, being flow evolution partly controlled by gravity leads indeed to the thinning and horizontal spreading of the water depth along the flow path.
Similarly to what observed in case of variable dam heights, the low frequency signal induced by the water drop has a variable frequency and spectral amplitude scaling inversely with water depth (Fig. 5c, d).  The behaviour of the frequency spectrum varying the water depth is confirmed also for other values of the dam height as shown in Fig. 6, where the dependence on the inverse of the water depth multiplied by the constant inlet velocity is highlighted as well.

Discussion and Conclusions
The comparison between the modeled flow evolution, as a function of dam height (Fig. 2) and water depth (Fig. 4), and the related pressure time series (Figs. 3, 5) allows us to define some basic features of infrasound energy radiated by a flow. Modeling shows that the generation of waves at the water free surface appears as a consequence of the water drop and is enhanced for larger water depths. Wavelength and amplitude of the waves are controlled by the dam height and water depth. Such waves are not produced in the absence of the dam (Fig. 2a), resulting into a smooth water surface. Waves at the water surface could be gravity waves produced by the impact of water within the dam absorption pool, Ostrovsky and Bedard (2002), and propagating downstream faster than the flow, or they could result from internal turbulence of the flow, investigated experimentally by Horoshenkov et al. (2013), which emanate from the bed where the flow interacts with the rough boundary, represented by the dam in our case, and are transported and shaped by the flow.
Despite the geometry of the modeled fluid is strongly simplified and the presence of sediment and boulders is not considered, the key results appear to be consistent with experimental observations of infrasound radiated by rivers and debris flow. In accordance with the analytical and experimental investigation by Feng. et al. (2014), who recorded infrasound from a dam that appeared to be consistent with surges induced in the absorption pool by falling water, modeled infrasound signal shows a clear transient induced immediately after the water impact downstream the dam. Additionally, modeling shows a longer lasting, and lower frequency, infrasound signal, that is scaling with flow parameters (dam height and water depth) and that develops only afterwards, when waves develop at the water free surface. Based on the specific timing of such signal and on the correspondence with modeled flow features, we suggest that it is likely related to the waves at the water surface produced by turbulent flow evolution downstream the dam, that represent an irregularity along the path (Horoshenkov et al. 2013). This interpretation is in general agreement with the experimental investigation by Schmandt et al. (2013), who explain indeed the infrasound produced by Colorado rivers in term of waves at the water free surface enhanced in irregular section of the flow path of the Colorado river (Hance Rapid), or by explanation provided by Marchetti et al. (2019) for debris flow, who explained the almost monochromatic, cigar shaped infrasound signal recorded during debris flow events at Illgraben as being produced by an extended source or multiple sources acting simultaneously, with dams representing preferential, but not unique, sources of infrasound energy.
Despite the proposed modeling is not considering sediment and boulder, and therefore still far from being able to predict debris flows, results appear to be able to explain to some extent the relation, observed experimentally, between recorded infrasound and discharge rate for debris flows. Infrasound experimental observation of debris flows shows indeed a linear relation between discharge rate and amplitude (Schimmel et al. 2018;Marchetti et al. 2019;Belli et al. 2021), consistently with modeling results. More recently, building on the measurements of debris flows at Illgraben recorded between 2017 and 2019, Belli et al. (2021) showed a clear inverse relation between peak frequency of recorded infrasound and discharge rate, with lower frequencies recorded for flows with larger velocity or flow depth. The same might not hold when modeling of seismic energy radiation is desired. Indeed seismic analysis performed with geophones positioned nearby the channel at the Gadria basin, Italy, (Coviello et al. 2019) clearly showed that peak energy is recorded during the passage of a debris-flow surge front, highlighting the role of the boulder rich front in seismic energy radiation. Thanks to the consistency between modeling results and existing experimental data we suggest that the main features of infrasound energy radiation by the flow are identified. This is critical both to investigate the source mechanism of infrasound by debris flows and to use infrasound to remotely characterize the source. This last aspect, in particular, will be required when the magnitude of an event will be inferred from infrasound data. Future studies will require to include in the modeling the transported solid phase (boulders) and the flow characteristic.

Boundary Conditions
In our model, the choice of the boundary conditions is critical for the simulation of the fluid flow. The imposed boundary conditions are able to reproduce some basic features of the fluid motion: water and air are free to enter and to leave the simulation domain; pressure and velocity are continuous at the boundary of the simulation domain. In particular, the package FixedFluxPressure is used in order to eliminate numerical effects generating spurious reflection waves (https://www.openfoam.com/ documentation/user-guide). The simulation domain X with boundary oX is depicted in Fig. 7  boundary oX is divided into six parts: airleftWall (oX 1 ), waterleftWall (oX 2 ), lowerWall (oX 3 ), air-rightWall (oX 4 ), rightWall (oX 5 ), atmosphere (oX 6 ).
The following boundary conditions are set for the phase fraction a, velocity vector U, the pressure field p and the temperature T.

Pressure Detection Points
The pressure data are analyzed in different points of the domain. In particular, we have compared the value of the pressure signal extracted from points close to the impact and to the end of the flow at different heights in vertical direction. No appreciable differences in the spectral data are detected. Moreover, the signal shows the same frequency content whether we consider the values of a single point or we take the mean over a volume. The pressure values are extracted at different points along the vertical and the direction at p=4 with the terrain. The pressure Power Spectral Density is showed in Figs. 8 and 9 for the cases related to the impact point.