Numerical Modelling of Flow-Debris Interaction during Extreme Hydrodynamic Events with DualSPHysics-CHRONO

: Floods can transport debris of a very wide range of dimensions, from cohesive sediments to large ﬂoating debris, such as trees and cars. The latter increases the risk associated with ﬂoods by, for example, obstructing the ﬂow or damaging structures due to impact. The transport of this type of debris and their interaction with structures are often studied experimentally in the context of tsunamis and ﬂash ﬂoods. Numerical studies on this problem are rare, therefore the present study focuses on the numerical modelling of the ﬂow-debris interaction. This is achieved by simulating multiple laboratory experiments, available in the literature, of a single buoyant container transported by a dam-break ﬂow in order to validate the chosen numerical approach. The numerical simulations are carried using the open source DualSPHysics model based on the smoothed particle hydrodynamics method coupled with the multi-physics engine CHRONO, which handles the container–bottom interactions. The trajectory, as well as the velocity of the centroid of the container, were tracked throughout the simulation and compared with the same quantities measured in the laboratory. The agreement between the model and the experiment results is quantitatively assessed using the normalised root mean squared error and it is shown that the model is accurate in reproducing the ﬂoating container trajectory and velocity.


Introduction
Tsunami events such as the one occurred in Japan in 2011 in Tohoku [1] and the ones in Indonesia in 2004 and 2018 [2] resulted in catastrophic consequences on coastal communities and structures causing many casualties. Studies in the literature often focused on the impact of tsunamis on the coastal region [3][4][5][6][7]. However, recent research investigated the effects that debris transported by such tsunamis can have on existing structures, especially after the damage documented during the 2011 Tohoku tsunami. Using the evidence found after this event a debris hazard classification was developed by Naito et al. [8] by analysing the damages caused by different debris types. Among the types with most energetic estimated impact forces on structures were shipping containers, especially loaded ones. Forces between 1500 and 6400 kN were estimated by Naito et al. [8] on structures by applying the FEMA P646 [9] guidelines for the design of structures. These recommend a formula for the calculation of the maximum force exerted by floating debris based on the Haehnel and Daly [10] relationship. In the wake of the 2011 Tohoku tsunami, further relationships were incorporated in the American Society for Civil Engineers (ASCE) [11] design codes and provide a simplified model for large floating debris impact forces for structural design, accounting also for impact orientation based on the mean impact orientation of Haehnel and Daly [10]. However, important factors such as debris material are not yet taken into account in any design guideline to date.
In order to have more detailed models for large floating debris dynamics, a large number of studies focused on investigating their behaviour when entrained in a flow, e.g., trajectories, spreading, orientation etc. [1,[12][13][14][15] or on analysing their impact forces during a collision with a structure [16][17][18][19]. More recently Derschum et al. [20] and Stolle et al. [21,22] studied both processes in conjunction by using a single "smart debris" [23] shaped as a 1:40 geometrically scaled 20 ft standard container to analyse both the trajectories and the impact forces on a structure with the aim of improving the reliability of design guidelines for structures that may be exposed to such hazards and to serve as benchmark tests for future numerical simulations.
The numerical modelling of these processes is rare in the literature [24] with significant limitations in the numerical methods used. Initial studies on solid-fluid interactions for debris were conducted by Wu et al. [25]. They used Navier-Stokes equations coupled with Volume Of Fluid (VOF) to track the free surface, an approach very often used in modelling fluid-structure impacts (see, e.g., De Finis et al. [26]), and a large eddy simulation turbulence model. The motion of the floating body was modelled with a Discrete Element Method (DEM) model validating it with laboratory experiments. However, implementations of floating body mechanisms on mesh-based models are usually very challenging, requiring ad hoc formulations for specific cases and long processing times [27]. Due to these difficulties, Lagrangian particle-based methods, such as the Moving Particle Semi-implicit method (MPS) [28] or the Smoothed Particle Hydrodynamics method (SPH) [29] are often used as good alternatives to mesh-based methods [30,31]. SPH, due to its formulation, is well suited to investigate high energy phenomena since it is capable of modelling violent flows as well as flow particle detachments in the context of tsunamis [32,33]. However, even if these methods produced accurate results for solid-fluid interactions [31,34] they do not inherently account for solid body contact laws, i.e., this contact requires a specific sub-or coupled model. Canelas et al. [35] and successively [36] expanded an SPH solver via coupling with a DEM model which was validated with PVC cubes subject to a flow, but showed limitations in correctly simulating the friction between moving and boundary elements of the simulations. This was further shown in Goseberg et al. [37] which modelled experiments from Nistor et al. [15]. The combined use of SPH and DEM clearly resulted in negligible variations of debris dynamics when changing the kinematic friction coefficient f c of the materials. Additionally, DEM implementation was also prone to numerical instabilities, also resulting in high computation times.
In order to address the aforementioned issues, the newly implemented DualSPHysics-CHRONO coupling in DualSPHysics v5.0 [38,39] is used to model the container trajectories measured in Stolle et al. [21,22]. This is, additionally, a first necessary step to study the feasibility of such approach and successively extending their use for flow-debris-structure interactions. DualSPHysics is chosen due to its advantage over other SPH numerical models of using the computational power of the Graphics Processing Unit (GPU) of a PC to speed up simulation times up to 146 times [38]. CHRONO [40,41] is a multiphysics engine that allows to solve rigid or deformable contacts and impacts between solids. In addition, it includes the possibility of integrating different types of constraints as recently done by Tagliaferro et al. [42] to simulate a planing hull. In contrast to the DEM implementation, CHRONO allows for a realistic frictional behaviour due to the inclusion of a full Coulomb sliding/sticking/rolling model. Additionally, the application of different types of constraints such as joints, hinges and springs allows to include complex structures allowing to realistically model the ones found in coastal regions or in harbours. Another advantage is the increased numerical stability gained when compared with DEM coupling.
The main objective of this work is to validate the coupling of DualSPHysics and CHRONO with the experiments presented in Stolle et al. [22]. First, the modelling of the hydrodynamics of the flow is validated, then a sensitivity analysis of the modelling of the details of the debris geometry and description of the material characteristics on the debris dynamics is carried out. Finally, the results of the calibrated models are presented and discussed.
The remainder of the manuscript is organized as follows. The explanation of the methodology, illustrating the numerical setup along with a brief introduction on the numerical methodology used is presented in Section 2. In Section 3 the main results of the initial hydrodynamics validation and subsequent calibration and validation of the debris dynamics are presented. Finally, in Section 4, the results are analysed, and the feasibility of this numerical approach is discussed in relation to potential applications for debris-structure interactions and conclusions of the work are given.

Methodology
For the present study, DualSPHysics v5.0 [29,38,39] coupled with CHRONO Engine are used to solve the hydrodynamics and simultaneously consider the collisions and interactions between container and floor.

DualSPHysics
DualSPHysics [39] is based on Weakly Compressible SPH (WCSPH) with the fluid phase governed by the Navier-Stokes equations, reduced to ordinary differential equations solved in a Lagrangian framework. The conservation of mass and momentum is expressed as [29,43]: where ρ is the fluid density, t is time, p is the pressure, v = v x , v y , v z is the velocity vector, g = g x , g y , g z is the gravitational acceleration vector (here g x = 0, g y = 0, g z = −9.81 m/s 2 is used) and Γ groups all the dissipative terms that can be defined by using two different formulations, i.e., artificial viscosity [44] or laminar with sub-particle scale turbulence formulations [45]. SPH discretises every part of the computational domain into sets of particles carrying different properties such as density, pressure, velocity, etc. In general, two steps are defined to apply Equations (1) and (2) to the SPH method, i.e., a kernel approximation and a particle approximation [46]. First, any variable f of a particle a, located at r a = (x a , y a , z a ), is represented by an integral at location r = (x, y, z) with Ω as the computation domain, W as a weighting function called smoothing kernel, which monotonically decreases with distance, and h p as the smoothing length which determines the size of the kernel support.
In the second step, the integral in Equation (3) is approximated by interpolating the characteristics of the surrounding particles as: where the summation is extended to all the particles inside the kernel. In Equation (4), W ab = W r a − r b , h p , and m b and ρ b are the mass and density of a neighbour particle b (located at r b = (x b , y b , z b )). In addition, any derivative of f can be expressed as where ∇ a indicates derivation with respect to the coordinates of particle a. Equation (1) is rewritten in the SPH framework for a particle a as A more extensive analysis of the SPH governing equations and model can be found in the reference literature [38,39,[46][47][48].

CHRONO Engine and Coupling with DualSPHyics
The multi-physics engine CHRONO [40,41], coupled with DualSPHysics, is here used to solve the interaction between the container and the bottom of the flume. In the problem at hand, i.e., the transport of a container by dam-break flow, the container first becomes floating during the initial mobilisation by the flow then it impacts with the bottom and, subsequently, slides along the flume floor.
CHRONO is able to consider multiple types of structural constraints. CHRONO solves collisions through two alternative formulations, (i) Non-SMooth Contacts (NSMC), which considers fully rigid impacts, and (ii) SMooth Contacts (SMC) which solves deformable contacts. For NSMC only the restitution coefficient K and the kinetic friction coefficient f c are used as variables when solving interactions between solids. SMC additionally considers Young's modulus E and the Poisson's ratio ν. In the present study (ii) is used to approximate as closely as possible the real behaviour for all materials and elements involved.
The coupling between the two models is defined with a message passing interface such that DualSPHysics first calculates the interaction between fluid particles and, subsequently, computes the fluid-floating interactions for a defined t. CHRONO is activated when mechanical constraints between floating bodies and/or material properties are defined as input of the simulation. In this case, both the linear and angular accelerations, together with the length of the time step ∆t in SPH are passed to CHRONO. This model computes and applies forces and momentum of the flow to the floating bodies by integrating the fluid contributions with the defined dynamics or kinematic restriction, including possible collisions. If the time step used for CHRONO is the same as for DualSPHysics, then the computation reverts to the latter with CHRONO communicating linear and angular velocities, as well as the position of the centre of mass of every floating body. In DualSPHysics the floating bodies particles positions, velocity and density are then updated considering CHRONO output. The entire system is then updated starting the cycle over again for t + ∆t.

Modelled Experiments
The experimental setup described in Stolle et al. [22] is modelled herein due to its focus on container dynamics. The experiments were conducted in a 30 × 1.5 m flume. A reservoir 21.55 m long was filled with an impoundment depth of h 0 = 0.4 m and the water was released via a swing gate to generate a dam-break flow on an 8.45 m long test area with the horizontal concrete floor elevated by 0.2 m from the flume bottom. Note that, only the first 6.50 m of the experimental test area were modelled because the impact with the structure at the end of the area described in Stolle et al. [21] is not the focus of the present study. The area modelled herein is shown in Figure 1. The swing gate structure consisted of two 0.05 × 0.05 m metal columns with an additional 0.03 m ledge towards the inside of the flume covered in rubber to ensure a watertight seal. This resulted in a 0.08 m protrusion on each side of the gate slightly obstructing the dam-break flow and generating three dimensional effects [49]. The container x-coordinate of the geometrical centre was positioned on the false floor, respectively, at x = 3.2 m from the end of the reservoir (see

Numerical Setup
The numerical setup closely reproduced the experimental conditions, with two exceptions: i.e., the swing gate was approximated by using an instantaneous release of the water in the reservoir at t = 0 s and the gate structure was modelled as two 0.08 × 0.05 m rectangular columns as shown in Figure 1. Note that the coordinate axis in the transversal direction of the flume is mirrored with respect to the one used by Stolle et al. [22]. First, the hydrodynamics was calibrated to define the best combination of parameters to accurately simulate the experimental hydrodynamic conditions. All the most relevant parameters defined after initial calibration are summarised in Table 1, for the remaining fundamental parameters default values were used. A final particle spacing (dp) of 0.01 m was chosen as 0.02 m resulted in a higher tip of the flow (due to the larger size of the particles). A convergence test showed that results for dp < 0.01 m converged to those for dp = 0.01 m at the expense of a significant increase in computational time by increasing from 28 × 10 6 particles with dp = 0.01 m to 54 × 10 6 particles with dp = 0.008 m. dp =0.01 m was also used as the dimensions of every element of the numerical setup is a multiple of dp, allowing correct element sizes. Here the relationship, h l = Coefh 3dp 2 , was used to calculate the smoothing length h l with Coefh being a calibration coefficient. Here Coefh = 1 was used since was already validated for dam-break flows. Note that, this value resulted in a ratio between the smoothing length and the particle spacing h l /dp = 1.73 (Table 1) which was also chosen by Tan et al. [32] to simulate landslide tsunamis and considered suitable for the problem at hand. In this work no sensitivity analysis on the effect of this ratio on the modelling was carried out.
In this study, to account for dissipation, i.e., to compute the Γ term in Equation (2), the artificial viscosity formulation [44] is applied since it was also successfully used to model other extreme hydrodynamic phenomena such as landslide-tsunami [52,53]. The artificial viscosity parameter between fluid particles α f f was defined as 0.04 after initial calibration to ensure a coherent velocity of the dam break in relation to the experimental values. An additional calibration was performed for the viscosity between fluid and boundary particles α f b = visc b f × α f f , with visc b f being a viscosity multiplication factor. In this study visc b f = 9 was used to approximate the roughness of the flume floor, which was covered in a paint/sand mixture, resulting in an estimated Darcy-Weisbach friction factor n = 0.014 [49]. visc b f was calibrated making sure that velocity with the distance and arrival time of the flow to each wave gauge was reproduced accurately. Additionally, a density diffusion term was added to the continuity equation, i.e., Equation (6) by using the formulation of Fourtakas et al. [51] on the fluid particles. This additional term is commonly used to deal with density oscillations that are distinctive of SPH often resulting in incorrect and unstable pressure on boundaries or floating bodies.
The container was simulated with dimensions of 0.15 × 0.06 × 0.06 m and a mass of 0.226 kg [22]. For the initial position in y-direction, the mean values of all experiments for each orientation were used. Since the mass distribution was not symmetrical in the experiments, the container was 3D modelled as shown in Figure 2. The container trajectories and velocities shown in the paper are referred to the centroid of the container. The inertia of the container was first calculated in DualSPHysics by modelling the container as in Figure 2 with a dp = 0.0005 m. This had the only purpose of letting DualSPHysics compute very accurately the inertia matrix (i.e., the 3 × 3 matrix of the moments of inertia). Since it was known that the mass distribution in the container was not symmetrical and the position of the centre of mass was not known from the experiments, this position was prescribed and calibrated with the experimental data. An eccentricity between 0.0025 m and 0.004 m in the longer dimension was considered in a series of calibration tests, with the latter value giving the best results and therefore used in the present study. The addition of the eccentricity resulted in a larger effect at the onset of the container motion with differences noticeable mainly in the y-direction; the differences in terms of centroid velocity along the y-direction (which have maxima of order of 10 −1 m/s) are up to 0.134 m/s at the motion onset and they become negligible after the container stabilised in the dam break flow. The centre of mass and inertia, taken from the high-resolution inertia matrix computation were subsequently prescribed at the start of each full-size simulation where a solid container was used due to particle size restrictions.

Properties of the Materials
Attention was given to the material properties (see Section 2.2) used in the simulations and summarised in Table 2, to ensure a realistic response to the debris-bottom interactions. High Molecular Weight Polyethylene (HMWPE) was used in the experiments [21] for which the design values of E and ν [54,55] are used in the simulations ( Table 2). The floor of the test area was made of concrete. However, since no information on the type of concrete was available, both E and ν were chosen to reproduce the characteristics of a low strength concrete, which is considered realistic due to its non-structural role in the experimental setup.
The remaining two parameters K and f c were varied in the ranges specified in Table 2. The reason behind this choice was twofold. First, these two parameters might vary with specific conditions of the experimental setup (e.g., lubricated friction) and second to understand their effects on the container dynamics.

Model Performance Assessment
The accuracy of the numerical simulation results of the flow depth (h) of the x and y trajectories of the container was assessed by using the root mean square error normalised with h 0 (nRMSE l ) as where ξ n,i and ξ e,i represent any i-th sample of one of the numerically modelled and experimental variables, respectively, and N is the number of samples. Furthermore, the accuracy of the simulated x and y velocity components of the container (v x and v y , respectively) was assessed with a root mean square error normalised with the shallow water velocity (nRMSE v ) as where v n,i and v e,i represent the numerical and experimental velocity components, respectively, in either xor y-directions. Two further metrics for assessing modelling performance are used in Section 4, namely the modulus of the mean error (|ME|) defined as |ME| = N ∑ i (ξ n,i − ξ e,i )/N (9) and the mean distance error (MDE) as where the subscript n and e, representing a numerical and experimental quantity, respectively, are used to further quantify the simulations performance by comparing them with the accuracy of the positioning system used in Stolle et al. [22]. The two parameters |ME| and MDE (following Equation (16) in Stolle et al. [56]) are chosen because they are the same ones used by Stolle et al. [56] to test the accuracy of the optical tracking algorithm used for the experimental measurements presented hereinafter.

Dam Break Hydrodynamics
The hydrodynamics modelling is validated by comparing the numerical time series of the flow depth with the mean experimental one (Section 2.3) measured by Stolle et al. [22]. Figure 3 shows the evolution of the normalised flow depth h/h 0 over t at every experimental WG position. The modelled emptying of the reservoir at x = −0.10 m (Figure 3a) follows closely the experimental mean, also reproducing the backward travelling perturbation due to the presence of the step and the gate structures, resulting in an nRMSE l = 0.010 for h. The evolution of h/h 0 in the remaining two WGs is less accurate with nRMSE l = 0.068 at WG2 (Figure 3b) and nRMSE l = 0.031 at WG3 (Figure 3c). At WG2 (Figure 3b) h/h 0 is better modelled for t ≥ 2.25 s whereas h/h 0 immediately after the arrival of the tip of the dam break flow shows a larger difference, this is due to the 3D flow effects, as mentioned in Section 2.3, not being fully captured by the numerical model. However, the "step" present in the experimental data of h/h 0 , caused by the gate structure protruding inside the flume cross section, is qualitatively reproduced. At the WG3 position, the close match between the experimental and numerical arrival time is essential for a good initialization of the container dynamics. This is because its centre in x-direction is positioned at the same location of WG3. The close match for 1.75 s ≥ t ≥ 2.5 s is also important as that part of the flow is responsible for the container entrainment in the dam break flow.
Additionally, the velocity of the tip of the flow is compared with the values calculated by Stolle et al. [22,57] by using the dam break arrival time. The first study found this velocity being 2.46 m/s when using a single container setup with the second finding it equal to 2.2 m/s but with multiple containers. The numerical flow velocity is calculated by directly using the output particle velocity of the fluid at WG3, resulting in a value of 2.53 m/s. Therefore, when compared to the experimental velocities only a percentage error of 2.8% is found. Finally, if compared to the analytical solution of Chanson [58] for which the bore velocity would results in 2.65 m/s only a 4.5% difference is found.

Sensitivity of the Container Dynamics to the Material Characteristics Parameters
A sensitivity analysis, as already introduced in Section 2.5, was performed for f c for the container and K to better investigate the effects of these material characteristics on the container dynamics, allowing to choose the best performing set of parameters for the case studied. When varying one of the two parameters the other was kept constant to exclude combined effects, i.e., f c = 0.3 and K = 0.7 were chosen for this purpose.
In every test involving the effect of f c only the value for the container was varied, maintaining f c = 0.3 for the flume floor. This choice was due to preliminary tests showing the flume kinematic friction coefficient not affecting the container dynamics, as expected. Figure 4 shows the difference in coordinate in x-direction (∆x) and y-direction (∆y) with respect to the container initial position for both investigated variables. This representation was chosen to better highlight the small differences between different values of the parameters that would not be otherwise as clear when using global coordinates. The effect of f c in xand y-directions is shown in Figure 4a  K is seen to affect the container trajectory less than f c for this study. This is due to impacts between the container and the flume floor only occurring at the onset of the container motion with only sliding between the two surfaces occurring for the rest of the motion. Note that the effect on the x evolution can be considered negligible in this case.

Container Kinematics
In this section, the container motion is analysed and compared with the experimental findings. Hereafter, only the results with the combination of the material characteristics that gave the most accurate representation of the experiments, after initial calibration, are presented. K = 0.7 was kept constant for both container and floor whereas a f c = 0.15 and f c = 0.30 was used for container and floor, respectively. Both OR1 and OR2 container orientations are presented. Only the results of trajectories projected on the x-y plane are shown as the objective of the paper is to validate the simulations with experiments focused on the container motion on a coastal area.

Trajectories
First, the container trajectories are analysed for both OR1 and OR2. Figure 5 illustrates the results for OR1 with Figure 5a showing the trajectory over the test area with the solid grey lines and solid black line representing all the experiments and their mean, respectively, compared with the numerical results represented with a solid red line. Note that the mean experimental values were calculated only up to when data was available for all experimental runs (see Figure 5b). The larger differences between the simulation and the experiments are found during the initial movement, for x ≤ 3.5 m, where the simulated container proceeds almost in a straight line in contrast to the real ones which are very often characterised by a transversal motion as shown also by the mean trajectory. However, at a later stage the simulated container also starts a meandering motion reaching a maximum value of y = +0.0447 m at x = 5.77 m, slightly further than the experimental one for which the same value occurs at x = 5.80 m. This behaviour of the container can be mainly associated with the not perfectly symmetrical distribution of the masses (Section 2.3) as the flow was mainly symmetrical with respect to the y-coordinate in the experiment (see Derschum et al. [20]) as well as in the numerical simulation. Figure 5b,c show the evolution of the trajectories components in x and y, respectively, over time. This helps to additionally identify the point of motion initiation which, although starting slightly before than in the experimental mean due to the slightly higher velocity of the flow, closely matches between simulations and experiments. In Figure 5b it can be noticed that the slope of the x (t) curve is nearly constant in the simulations, while in the experiments it varies, indicating differences in the velocity as discussed in the next subsection. The performance was measured with the nRMSE l , i.e., Equation (7) for both directions and summarised in Table 3. The results for OR2 are shown in Figure 6 and compared with the corresponding experiments, as done for OR1. The numerical trajectory (Figure 6a) follows the mean experimental trajectory but with a divergence from the mean values on the y-coordinate for x > 4.5 m. Additionally, the numerical model is not able to reproduce the experimental container trajectory at the onset of the motion, first towards and then away from the centre of the flume. These differences, which are particularly noticeable in Figure 6c, also affect the nRMSE l (Table 4) in y-direction with a higher value than that for OR1.

Velocity
Hereafter, the velocity of the centroid of the container is analysed to give a deeper insight into the accuracy of the used numerical method in simulating debris motion. As velocities were not directly measured in the laboratory, in this study they were calculated as the time derivative of the trajectory in both x and y directions. The performance of the modelling of the velocity was assess using Equation (8) and the results are shown in Table 5. Figure 7a,b shows the container velocity in x-(v x ) and y-directions (v y ) for OR1 whereas Figure 7c,d show the results for OR2. The velocity results are more accurate for OR1 than for OR2, even if slightly, as shown in Table 5. Here, the initial container v x increases almost instantly, probably also due to the deeper tip of the flow, contrarily to the more gradual increase measured in the experiments affecting its motion initiation. However, the final simulated velocity follows closely the experiments with only a relative difference of −4.0% with respect to the experimental value. The opposite is shown for OR2 (Figure 7c) where the gradual v x increase is similar to the experimental one, but the final velocity results in a larger relative difference of −9.8%, again with respect to the experimental value. Finally, v y shown in Figure 7d confirms what shown in Section 3.3.1 with an initially negative v y becoming positive at t =1.5 s in the experiments, not being captured by the numerical model where a positive v y is always computed. Finally, v y for both cases tends to be more accurately predicted after the container stabilisation with the flow, which occurrs after some time from its initial mobilisation.

Qualitative Description of the Interaction between Container and Flow
A visualisation of the interaction between flow and container is shown for OR1 and OR2, simulated using the validated model, to qualitatively highlight the key differences of the two orientations. Figure 8 shows a top-down representation of different instants during the simulation. For OR1 the dam break flow impacts right after t = 1.125 s (Figure 8a). The interaction after the initial impact starts to develop and, as shown for t = 1.5 s (Figure 8b), the container starts to rotate counterclockwise due to its eccentricity while the adjacent flow moves faster than the container creating a "u" shape around it. Afterwards, the flow and the container movements start to stabilise (Figure 8c, t = 3.0 s). Here, the flow starts to close around the container, which rotates slightly clockwise while transported. Note that these rotations affect the transversal trajectory of the container and, in turn, the symmetry and flow direction of the tip of the dam break.
Some differences are found when analysing the interaction between flow and container for OR2. The main difference with OR1 is that, due to the smaller impact area, the flow forms a "v" shape separating from both sides of the container. The container still rotates clockwise, starting slowly at around t = 1.5 s (Figure 8e) but continues to increase with its motion up to t = 3.0 s instead of changing direction as seen in OR1. Contrarily to what happens for OR1, due to the slightly slower container v x and different orientation, the dam break flow manages to close around the front part of the container more quickly (Figure 8f).

Discussion and Conclusions
This section discusses the findings of this study with the focus on analysing the future applicability of this numerical methodology in a broader context, as introduced in Section 1. First, the accuracy of the simulation results was compared with the accuracy reached by the container experimental measurement system as additional benchmark. Stolle et al. [22], as measurement system, employed an optical tracking algorithm introduced in Stolle et al. [57] where it is investigated its accuracy using the parameters introduced in Section 2.5. The mean values of the parameters calculated by Stolle et al. [57] for the experimental validation tests with a single container are summarised in Table 6. |ME| and MDE values for the simulations of the present study are additionally calculated for both OR1 and OR2. This helps to highlight the very similar accuracy of the present numerical model and the experimental tracking system in the description of the container trajectory giving further confidence on the employed numerical method. In conclusion, the coupled models DualSPHysics and CHRONO can accurately simulate the behaviour of a debris transported by a dam break flow. Note that the modelling accuracy is referred to the mean trajectory among the ones resulting from repetitions of the experiment. This is because the numerical model was calibrated using the mean flow and container kinematic conditions, which was necessary given the highly variable characteristics of this phenomenon. Additionally, the calibration process of the model has shown that the details of the container internal structure and, above all, mass distribution is essential to achieve accuracy due to the importance of the container rotation in its kinematics. This rotation is, in fact, responsible for the container drift from the flume axis as in its absence the container would otherwise follow an almost one-dimensional trajectory. The difference in drift from the flume axis between simulations and experiments during the initial stage of the container motion, in turn affecting the velocity, suggests a lower accuracy in the forces if an impact would occur during this stage. However, as the trajectory stabilises the accuracy significantly improves especially for OR1. This has important consequences in modelling the impact of these debris on structures, which is the natural next step of this investigation. It is expected that impact forces are more accurately predicted if the impact occurs at later stages of the transport.