## Appendix A  
In the present study, the modelling and detection of particle-particle collisions was computed using a deterministic collision model proposed in Reference [47]. This model assumes that the collisions are binary and quasi-instantaneous. Furthermore, the elastic contact between two hard spheres occurs in a single point. In order to determine whether a collision occurs, the time-sequenced collision detection algorithm described in detail in Reference [38] was applied. Such algorithm is briefly summarised here.  
Firstly, all particles are tracked to the end of the present Lagrangian particle time step, \(\Delta t\). Inter-particle collisions have occurred if particles overlap after the actual time interval. In the case of a collision, the distance between their centers is equal to or less than the half sum of the two particle diameters. The determination of all existing collision pairs is realised under application of the cell-index method [47]. The potential collision partners of a given particle are found among the other particles in the cell plus the particles in all 26 neighbouring cells. The list of particles in each cell is set up and updated using linked lists. Assuming that the particles move at constant velocities during the particle time step, the collision time for each overlapping pair, \(\Delta t_{ijc}\), is calculated by solving the following quadratic equation:  
\[\left|x_{pij}(t+\Delta t)+\mathbf{u}_{pij}(t+\Delta t)\Delta t_{ijc}\right|=0.5 \Big{(}d_{pi}+d_{pij}\Big{)} \tag{1}\]  
Here, \(x_{pij}\) and \(\mathbf{u}_{pij}\) are the relative position and velocity of particles \(i\) and \(j\) which have diameters \(d_{pi}\) and \(d_{pij}\), respectively. Equation (1) suffers from the limitation that it will not detect a particle-particle collision if complete interpenetration occurs during the given particle time step. Subsequently, the determined collision times (i.e., the time between the beginning of the actual time step and the time of particle collisions) are sorted in descending order. In the time-sequenced implementation, the first overlapping pair is moved backward in time to the position of collision. The determination of the post-collision velocities of the collision couple is based on the laws of impact for hard spheres provided in Reference [48]. A detailed derivation of the used collision operators can be found in Reference [38].  
After the collision is computed, both particles are again moved forward in time to complete the actual particle tracking time step. New overlaps may be added to the list of colliding particles, i.e., in the case of multiple particle collisions within one tracking time step, while processed or non-existent ones are deleted. These collision computations are carried out until no overlaps are left.  
The computational cost of the searching collision partners can be estimated as follows. In the employed Cell Index Method, the simulation box is divided in a regular array of \(M^{3}\) cells; on average, there are \(N_{p}/M^{3}\) particles (\(N_{p}\) is the total number of particles) within each cell, where 27 neighbouring cells have to be searched per particle. Therefore, avoiding double counting, the search algorithm requires \(\left(\frac{27}{2}\right)N_{p}/M^{3}\) operations [47].