A granular energy-controlled boundary condition for discrete element simulations of granular flows on erodible surfaces

Entrainment of the substrate, which can cause the volume of geophysical flows to increase by many times and greatly affects its mobility, occurs frequently. Even though many large-scale studies have been conducted to investigate entrainment during flows, particle-scale studies are urgently needed to understand this complex process. We propose a novel configuration with an energy-controlled boundary to model inclined granular flows on erodible surfaces with the aim of investigating the rheological features of granular flows with entrainment. The erosion degree is controlled by adjusting the energy of the bottom granular bed. The bulk velocity, stress, microstructure, and constitutive behavior are systematically analyzed. Our simulations capture the exponentially decaying velocity in the creeping region, confirming that this configuration can be applied to study geophysical flows with entrainment. The velocity at the free surface increases with erosion depth, which means entrainment can greatly enhance the mobility of granular flow. However, the flow in the shear band still belongs to the inertial regime. The novel energy-controlled configuration provides effective and efficient way to investigate granular flows on erodible surfaces, serving as an initial step to understand geophysical flows involving entrainment, which is essential to develop new geo-disaster mitigation strategy.


Introduction
Geophysical flows, such as avalanches, debris flows, and pyroclastic flows, which are extremely hazardous because of their high mobility and can lead to huge casualties and property losses, frequently occur in mountainous regions worldwide (Hungr, 1995;Lai et al., 2021;Roche et al., 2013). Field investigations have found that the final volume of geomaterials in such flows is significantly larger than the initial volume (Kang et al., 2016;Liu et al., 2019;Lo and Chau, 2003). For example, the initial volume of failed rock at the source region of the Xinmao avalanche in 2017 in China was only approximately 3.0 × 10 6 m 3 ; however, the final volume reached 6.3 × 10 6 m 3 , twice the initial volume (Liu et al., 2019). Another representative case is the Tsingshan debris flow, whose initial and final volumes were 400 m 3 and 20,000 m 3 , respectively (Lo and Chau, 2003), with the final volume being approximately 50 times larger than the initial volume. While the fragmentation and expansion of geomaterials are possible reasons for this phenomenon, the entrainment of the erodible substrate along the path of a flow may also be a critical reason for the increased volume (Hungr and Evans, 2004;Iverson and Ouyang, 2015;Mangeney et al., 2007). Entrainment-induced enlargement of the volume can strongly affect the mobility of a geophysical flow, which may ultimately result in more serious damages. It is therefore important to study the flow behavior of geophysical flows involving entrainment.
There are two typical entrainment mechanisms including basal erosion and frontal entrainment or ploughing in geophysical flows. Basal erosion, which scrapes materials from the basal surface along its flow path, happens when geophysical flows run across slopes composed of well consolidated geomaterials with high shear strength (Sovilla et al., 2006). However, when the geomaterials along the flow path are featured by low-density (e.g. less consolidated, recently deposited or partially fluidized), frontal entrainment or ploughing, where the flows front may impact the loose materials to increase volume, will happen. Front entertainment is the dominant entrainment mechanism in snow avalanches (Sovilla et al., 2006;Köhler et al., 2016), but basal erosion is a more common process in geophysical flows. Hence, basal erosion is the main focus of this paper, while front entertainment is beyond the scope of this work.
Because entrainment can cause the volume of a flow to grow by many times and greatly affects the mobility of geophysical flows, studies investigating entrainment during flow are cutting edge and have gained a considerable amount of attention worldwide in the field of engineering geology. Hungr and Evans (2004) studied two typical rockslides (i.e., the Eagle Pass Slide and the Nomash River Slide) using both field surveys and numerical simulations. A large volume of the substrate was entrained in both cases based on field mapping, and the dynamic evolution processes of the slides and the effect of entrainment on their mobility were analyzed and evaluated using the dynamic analysis of landslide (DAN) code. Roche et al. (2013Roche et al. ( , 2016 systematically investigated entrainment in pyroclastic density currents (PDCs) using both model tests and theoretical analyses. They found that the long-lived upward pore pressure gradient is the main driving force for entrainment and proposed an equation to calculate the velocity of a PDC. They also found that a large initial velocity is not a necessary condition for a long runout of a PDC, instead the entrainment of the substrate during the flow matters. A depth-averaged continuum method with an entrainment model was employed to reproduce the flow process of the Xinmo landslide and to assess the entrainment effect (Liu et al., 2019). The entrainment and geomorphic characteristics of the Chada avalanche were studied using geological and aerial surveys and drilling (Lai et al., 2021). A three-dimensional numerical method was also adopted to simulate its evolution process. The above large-scale studies are significant with respect to understanding the basic processes and effects of entrainment in geophysical flows. However, there is still a lack of particle-scale studies, which can reveal both the macro-and microcharacteristics of entrainment, and whose results can significantly advance large-scale studies considering the inherent granular attributes of geophysical flows.
From the perspective of granular rheology, geophysical flows involving entrainment are granular surface flows on heaps featuring an exponentially decaying velocity profile; this is a challenging problem because of the coexistence of both liquid-and solid-like regimes. Daerr and Douady (1999) reported two distinct types of granular avalanches when a static layer of glass beads on an inclined plane is perturbed. For a thin layer, an avalanche propagating downhill and laterally is observed, while uphill propagation is also presented for a thick layer. Pouliquen and Forterre (2002) experimentally studied the spreading of particles along rough surfaces and also developed a depth-averaged model with a friction law for this complex process. Gray et al. conducted systematically experimental research of granular flows on erodible beds using a stable inflow to study erosion-deposition waves (Edwards and Gray, 2015) and the release of different number of particles to study formation of levees (Edwards et al., 2017) and travelling waves (Viroulet et al., 2019). Based on these experimental results, a depth-averaged framework was also developed to describe the complex behaviour. The depthaveraged framework was also applied to investigate self-channelization and levee formation of monodisperse granular flows on non-erodible slopes (Rocha et al., 2019). A depth-integrated equations by balancing the kinetic energy is proposed to model entrainment and detrainment without erosion law and verified using experimental data (Capart et al., 2015). de Ryck et al. (2010) employed both numerical and theoretical methods to study granular flow on both 'frozen' and 'non-frozen' heaps. The frozen heap has static particles, while the non-frozen heap refers to a heap whose granules even in the deep of the heap present creeping motion. They found that the erodible bed in non-frozen case can greatly increase the mobility of a granular flow. However, their study only involved two extreme cases and did not consider different erosion degrees, making the study insufficient with respect to systematically understanding the effect of an erodible bed on the granular mobility. Moreover, the setup of a whole heap model requires a considerable number of particles, which will entail a huge computational cost. An inclined surface flow model with periodic boundaries in both the streamwise and lateral directions can simulate granular flow in an infinite domain with a finite number of grains (Silbert et al., 2001;Zhu et al., 2020). However, this setup can only access either liquid-or solidlike flows by adjusting the angle of inclination. To model the coexistence of both liquid-and solid-like regimes, two frictional sidewalls have been applied on the lateral boundaries in previous studies (Berzi et al., 2020;Bi et al., 2005;Jop et al., 2005). Even though both liquid-and solid-like regimes can be obtained in such a case, the complex boundary effect from the frictional sidewalls makes analyses difficult and the narrow dimensions in the lateral direction further limit the application of the results to real geophysical flows. In this paper, we propose a novel method to control the erosion depth of the bottom granular bed by tuning its energy. We then systematically study the bulk velocity and stress of the entire bed and evaluate the effect of the erosion degree on the rheological properties of the flow.
The rest of the paper is arranged as follows. The discrete element method (DEM), numerical setup of the inclined granular flow, controlling method of the erodible bottom granular bed, and post-processing method are introduced in Section 2. The proposed numerical scheme is verified and validated in Section 3. The numerical results are then presented and discussed in Section 4. Finally, our conclusions are summarized based on the analysis and future remarks are presented in Section 5.

Numerical methodology and conditions
We conducted numerical simulations using the DEM implemented in the open-source Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code (Thompson et al., 2022). The motion of each particle is tracked based on the Newtonian motion equations by considering the contact forces between the particles. In this study, only repulsive, finite-range contact forces described by a linear springdashpot model are considered. The normal F n ij and tangential F t ij forces between the particles i and j are determined such that given elastic constants k n and k t , viscoelastic damping constants γ n and γ t , an overlap distance δ ij , relative velocities v n ij and v t ij , a unit vector n ij , a tangential displacement vector u t ij , and an effective mass m eff . The Coulomb friction criterion (i.e., To generalize the implementation of our numerical findings, the results are presented in nondimensional quantities. The particle diameter (d), density (ρ), and gravity (g) are selected as control variables. The scalings for time, velocity, stress, and the elastic constant are , gρd, gρd 2 , respectively (Silbert et al., 2001). In this study, the parameters are set to k n = 200, 000 gρd 2 , k t = 2/7k n , γ n = 50 ̅̅̅̅̅ gd √ , and μ p = 0.5.
Bulk granular flow on a rough inclined plane (i.e. inclined chute flow) is frequently encountered in geophysical flows (GDR MiDi, 2004;Silbert et al., 2001). The configuration of inclined chute flow is simple: a granular bed flows on frictional inclined bottom boundary. The simple geometry and well-controlled condition make the inclined chute flow a widely-used configuration for studying fundamental rheological properties of granular materials (GDR MiDi, 2004;Pouliquen, 1999;Pouliquen et al., 1997;Silbert et al., 2001;Zhu et al., 2020). The inclined chute flow, a simulation snapshot of which is shown in Fig. 1, is thus adopted for this work. The three-dimensional model is established with periodic boundaries in the x-and y-directions. The implementation of periodic boundary condition can be regarded as surrounding the domain with replicates of the configuration. The particles close to periodic boundary will interact with particles of the replicas and particles will reappear on the opposite side when they cross the model box. The periodic boundary condition enables the simulation of an infinite domain using a finite number of particles without solid boundary effects. The bottom boundary in the z-direction is a rough layer with a thickness of 2d formed by particles whose motion was 'frozen' (the black particles in Fig. 1). The upper boundary in the z-direction is a shrink-wrapped boundary, where the limit of the box is adjusted to encompass the whole particles. However, there is no interaction between the particles and the shrink-wrapped boundary. The shrink-wrapped boundary provides a way to determine the location of the free surface of the granular flow at steady state. The dimensions are 20d and 10d in the x-and ydirections, respectively. The inclination angle is θ = 26 • , and the gravity is vertically downward. In geophysical flows, the particles underneath the free-flowing surface layer are shown to be perturbed by the surface flow, but dissipate the kinetic energy from the flow so that they transit from being quasistatic to being static at the bottom. To capture the behaviour of those particles and their interaction with the flowing particles usually requires simulating the full system or introduce additional wall boundaries to constrain the particles in a smaller domain (Berzi et al., 2020). Inspired by the dissipative nature of the bed particles and to reduce the computational cost, we employ a domain with periodic boundaries to simulate one part of the flow on an erodible bed, as shown in Fig. 1 and control the kinetic energy,E, of the bottom region (the blue region in Fig. 1). The kinetic energy E is defined as where N is the number of particles at blue region, m i and v i is the mass and velocity of particle i. In each timestep, the particles in the blue region will be detected and their kinetic energy are then calculated. When there is a difference between measured E m and target or controlling energy E c (i.e. |E m − E c |〉E w for E w = 1 %E c ), the velocity of particles in the blue region will be rescaled accordingly (Eq. (4)) to ensure the desired energy is achieved.
where v scaled i is scaled velocity of particle i. The specific flowchart is presented in Fig. 2. The simulations are run until a steady state is reached. The steady state is determined by checking the time history of the kinetic energy across the entire granular bed. To the best of the authors' knowledge, this is a new strategy to model heap flow using energy-controlled protocol.
The thickness of the energy-controlled region is 50d, which is demonstrated to be sufficiently deep to capture the salient feature of an exponentially decaying velocity profile in a heap under the surface flow (GDR MiDi, 2004;Komatsu et al., 2001), as shown in Fig. 3. The creeping bulk can be recognized in the region where μ < μ c (Liu and  Henann, 2017). μ is stress ratio and μ c is critical stress ratio and is only affected by the microscopic particle friction. A non-zero particle velocity in the creeping region is observed for all cases, and the magnitude of the velocity decreases with the controlling energy E c . An exponential function, V x = 0.2785 V 0 e 0.0425H , where V 0 is the boundary velocity between the upper inertial flow and the bottom creeping flow, was used to fit the velocity profiles with high fitting goodness when the controlling energy E c /E max is between 2.56E-2 and 2.56E-4. However, deviations can be found in the bottom region close to the solid boundary with thicknesses of approximately 10d. Such deviations are induced by the solid boundary via a cooperative effect and have been observed in inclined chute flows (GDR MiDi, 2004). Therefore, heap-like flows are well recovered using this novel configuration with an energy-controlled region. The exponential equation fails to fit cases when E c /E max < 2.56E-5, under which conditions a roughly plug flow occurs in the middle region. This is reasonable because the cooperative effect becomes stronger under such conditions, enabling more particles to feel the boundary. The confirmed heap-like properties in our configuration when E c /E max > 2.56E-4 indicates that this work provides a reliable tool to understand geophysical flows on an erodible bed.
We used two specific cases to set the upper and lower limits of our simulations. In case (1), the kinetic energy is free to evolve without external control. The flow is a standard inclined chute flow with an initial height of 90d. This is a case that has been well studied (Silbert et al., 2001;Zhu et al., 2020). All the particles are involved in a granular flow, and a Bagnold flow is reached at the steady state. In case (2), the kinetic energy is set to 0 (i.e. E c = 0). The particles in the blue region are frozen and form a rough boundary. The upper particles flow as a Bagnold flow with a smaller number of particles than in case (1). The particles in the blue region evolve from a frozen state to a free flow state by tuning its kinetic energy; these particles present heap-like behavior, enabling us to study the effect of different bottom heap-like flows on the upper inclined granular surface flow.
The rheological data at the steady state are extracted using a timeaverage and space-average method. The space average is conducted using bins with a thickness of 1d in the z-direction across the x-y-plane. Then, the time average is taken to output the final profiles. The stress tensor σ is calculated such that given a particle centre-to-centre vector r ij , a particle contact force F ij , a particle mass m i , a particle velocity v i , an average velocity v in the sampling region, and a sampling volume V.

Verification and validation
To verify and validate the numerical scheme used in this study, we compare our numerical results with theoretical data for the case without energy control at steady state. The stresses of the granular assembly balance its gravity at steady state. Hence, the theoretical stresses are where σ zz , σ xz , μ, g, ρ, φ, θ, H s , and z are the normal stress, shear stress, stress ratio, gravity, density, volume fraction, inclination angle, height of the granular bed, and coordinate vertical to the flow direction, respectively. The normal and shear stress profiles are linearly distributed across the entire bed with 0 at the free surface based on the above equations. Moreover, the stress ratio is constant at steady state. The velocity profile V x (Silbert et al., 2001) is where A Bag is the strain rate amplitude, which is set to 0.345 in this case. The velocity and stress ratio profiles are compared in Fig. 4, where the blue circles and the red line indicate the DEM simulation results and the theoretical results, respectively. Good agreement across the entire granular bed is observed, which confirms the accuracy of the numerical scheme used in this study.

Results and discussion
The case without a controlling energy E c , whose energy at blue region E max = 585 ρgd 4 , free surface velocity v xmax = 120.6 To verify our numerical protocol, we also design a case by imposing the energy of 585 ρgd 4 (E c /E max = 1.0). The profiles (pink lines) of velocity, volume fraction (Fig. 5), stress and stress ratio (Fig. 6) for the E c / E max = 1.0 case agree well with those (black solid lines) in the Bagnold flow (i.e., the Free case), which proves the reliability and accuracy of our energy-controlled scheme. We can clearly observe that the free surface velocity increases with the bottom controlling energy. The free surface velocity can increase up to three times compared with the solid bottom (i.e., E c /E max = 0). The increasing trend of the velocity agrees well with a previous study (de Ryck et al., 2010), in which only the two special cases of E c /E max = 0 and Free were simulated. The above results indicate that the erodible bottom boundary can greatly enhance the mobility of a granular flow, which agree well with previous studies (de Ryck et al., 2010;Iverson and Ouyang, 2015;Roche et al. 2013Roche et al. , 2016. Moreover, the velocity in the bottom energy-controlled region approaches zero (i.e., a solid-like state) when E c /E max < 2.56E-3. However, the upper region always presents fluid-like behavior with large velocities. This phenomenon indicates that the regime of the granular bed evolves from a fluid-like regime to a solid-fluid-like coexisting regime. Finally, note that the steady height of the granular bed H s also monotonically increases with the energy (Fig. 5  (c)), indicating an increasing number of particles in the bottom region gaining velocity and beginning to flow.
Indeed, the volume fraction profiles (Fig. 5(b)) confirm that more particles are engaged in the flow because the bottom volume fraction clearly decreases with energy. The volume fraction, which is an important parameter used in granular rheology, is defined as the ratio of the total solid volume to the sample volume. The volume fraction is constant across the entire granular bed when the bottom energy is not controlled (the black line in Fig. 5(b)); this has also been reported in previous studies (Silbert et al., 2001;Zhu et al., 2020). Interestingly, the volume fraction profiles above the energy-controlled region collapse to a single line. This phenomenon indicates that the rheological features of the upper region are not affected by the bottom region even though the upper velocity is greatly increased. However, the bottom volume fraction roughly decreases with height and gradually evolves to the constant volume fraction of upper region at the interface when energy is controlled, except at extremely low energy (e.g., E c /E max = 0 and 2.56E-6) where a jump in volume fraction is observed. With the increase of E c / E max , the bottom region becomes more erodible, and the depth of liquidlike part also increases. The volume fraction evolving from the bottom dense solid-like region to the top liquid-like part thus presents a decreasing trend at the interface. The bottom granular bed is frozen like a solid without any movement when E c /E max = 0, in which case the nearly constant volume fraction (the blue line in Fig. 5(b)) confirms the initial homogenous state of the granular bed. Meanwhile, a low energy of E c /E max = 2.56E-6 leads to limited movement of the particles, which means that a much longer time is required for particles to rearrange and reach the final steady state. This may be the main reason for the irregular volume profile (the red line in Fig. 5(b)) when E c /E max = 2.56E-6.
To better understand the bulk behavior of granular beds under different conditions, the scaled stress profiles calculated using Eq. (5) are presented in Fig. 6. For the E c /E max = 0 case, we ignore the stress in the bottom region, which is determined by the initial packing state. The stress profiles in the region with and without a controlling energy are remarkably different. In the upper region, where the energy is not controlled (i.e., H > 50d), the steady stress monotonically increases with a constant gradient and the stress ratio remains the same under the different energy conditions. This phenomenon can be understood by referring to the theoretical analysis. For the upper region, gravity is the only external factor affected the bulk behavior. Therefore, the profiles of the shear and normal stresses and the stress ratio at the steady state follow Eqs. (6)-(8), respectively. The different magnitudes of the stress at the same height in the different cases is caused by the differing location of the free surface. The similar bulk stress properties further confirm the similar bulk rheological behavior in the upper region.
However, the profiles of the bulk stress in the bottom region (H < 50d) present different characteristics. A monotonically increasing trend of the shear stress is clearly observed for the different cases. However, there is no strong linear distribution, as in the upper region. Moreover, the shear stress at the same height significantly increases with the controlling energy. However, the bulk normal forces in the bottom region roughly collapse to a single line. The different stress ratio profiles are therefore primarily induced by shear stress. This phenomenon signals that the effect of the energy on the shear and normal stresses appears to be different. Specifically, the bottom shear stress is greatly affected by the controlling energy, while the normal stress is not obviously influenced. To understand this, we need to keep in mind that both controlling energy and gravity are external factors in the bottom region and can jointly affect the stress. At steady state, both the vertical (V z ) and lateral (V y ) velocities are very small compared with the flow (V x ) velocity. Therefore, the controlling energy has a greater effect on the steady flow velocity than on the other two velocities, significantly altering the stress in this direction, i.e., there is a significant effect on the shear (τ) stress but a limited effect on the vertical (σ zz ) and lateral (σ yy ) stresses. The different bulk stress profiles indicate that the orientation of the mechanical contact is different for the different cases.
To better understand the bulk stress profiles at the particle level, we further explored the microstructure of the granular bed. The distribution of the pairwise contact force (defined in Eqs. (1) and (2)) (Ness and Sun, 2015), which illustrates the radial direction of the force in the x-z-plane, is presented in Fig. 7. The angle between the pairwise contact force and the x-axis represents the orientation of the particle contact force. Fig. 7 (a) shows the radial contact force distribution of the entire granular bed at steady state. There is an obvious alignment of the contact force with obvious differences between the different cases. Specifically, the percentage of pairwise contact force at main orientation is decreasing with the decrease of controlling energy, but it increases at the direction normal to main orientation (see the arrows in Fig. 7(a)). Therefore, the isotropy of the microstructure increases with decreasing bottom energy, except in the E c /E max = 0 case, in which only the microstructure of the upper region is considered. The steady inclined granular flow (i.e., E c / 6 E max = free) is the Bagnold flow (Zhu et al., 2020). Hence, the above alignment of the contact force is expected (Da Cruz et al., 2005). With decreasing bottom energy, a greater proportion of the bottom region presents solid-like or quasi-static features, whose microstructure is more isotropic (Ness and Sun, 2015). Therefore, the isotropy of the microstructure of the entire bed increases with decreasing energy. We separately plot the microstructures of the upper and bottom regions in Fig. 7 (b) and (c), respectively. In the upper region, the different lines collapse to a single microstructure profile (Fig. 7(b)) because gravity is the sole external factor in this region. This consistent microstructure again confirms the similar stress profiles in upper region (Fig. 6). Moreover, the different microstructures in the bottom region (Fig. 7(c)) are expected because the coexisting fluid-and solid-like regimes and the proportion of these different regimes vary with the bottom controlling energy.
To quantitatively identify the fluid-and solid-like regions in the granular bed, we use the shear zones to find the boundaries between the regions. In this study, the shear zone where fluid-like behavior is observed is defined as the region where the stress ratio μ is larger than a critical value μ c ; this type of criterion has also been used to recognize shear zones in vertical chute flows (Barker et al., 2022). μ c is only affected by the microscopic particle friction and is 0.382 when μ p = 0.5 (Chialvo et al., 2012). The particles in the upper red region do not experience a controlling energy; therefore, they definitely present fluidlike properties. We therefore focus on investigating the shear zone height H shear in the bottom energy-controlled region; the normalized shear zone height H shear /H bottom is shown in Fig. 8. H bottom = 50d is the height of the bottom energy-controlled region. We can clearly see that the bottom shear zone height increases with the controlling energy. When E c /E max = 0, the entire bottom granular bed is in a solid-like state; this state switches to a fluid-like state when E c /E max > 0.256. Between these two points, the state is that of a solid-fluid coexisting in the bottom region with the shear zone height monotonically increasing with the energy. The aforementioned phenomenon is reasonable because the bottom region flows more easily with increasing controlling energy, ultimately enhancing the velocity of the granular bed as shown in Fig. 5. To study the rheological behavior of the shear zone, we analyzed the relationship between the dimensionless pressure pd |Φ− Φc| 2/3 k (i.e., p * /k) and the dimensionless shear rate given the pressure p, the volume fraction Ф, the critical volume fraction Φ c , and the shear rate γ in Fig. 9. Chialvo et al. (2012) defined three regimes (i.e., the quasistatic, intermediate, and inertial regimes) of steady granular flow at all volume fractions and shear rates and proposed a constitutive model (the master curve in Fig. 9) to bridge these different regimes. The rheological data of the shear zone for all cases collapse to the inertial regime line, except at several points close to the free surface. This phenomenon confirms that the constitutive behavior of the shear zone is not affected by the controlling energy and can be well described using the classical model proposed by Chialvo et al. (2012). The rheological data of the entire granular bed under different conditions are presented in the inset of Fig. 9. We can clearly observe that these data collapse well and present solid-like characteristics (i.e., the quasistatic regime) where the pressure is independent of the shear rate. However, the magnitude of the dimensionless pressure is below the master curve. Moreover, the rheological data of vertical chute flow (squares in the inset) (Barker et al., 2022) can well collapse to the energy-controlled inclined chute flow. The lower dimensionless pressure is still under investigation and perhaps due to the coexist of both liquid-and solid-like regimes.
To further demonstrate the usefulness of the energy-controlled method in studying the rheological behavior, a power-law scaling in granular rheology involving the stress ratio μ, the inertial number I and a dimensionless temperature Θ = ρT/σ zz (Kim and Kamrin, 2020) is introduced, where T is the granular temperature T = δv 2 D with δv being the velocity fluctuation and D being the space dimension. For the Bagnold flow (E c /E max = 1.0), the dimensionless temperature Θ profile (pink line in Fig. 10 (a)) is constant across granular bed. With the decrease of E c /E max , Θ decreases at energy-controlled region, but keeps constant at upper region (H > 50d) and collapses to Bagnold flow. Moreover, roughly exponential decreasing trend of Θ with height is observed at energy-controlled region when E c /E max < 2.56E-3. The relation between μ, Θ and inertial number I is shown in Fig. 10(b), where the data collapse well into a single curve, consistant with the previous findings for nonuniform flows (Kim and Kamrin, 2020). This indicates the granular temperature is useful for construction of constitutive models for flows with coexisting liquid-and solid-like regimes.

Conclusions
In this paper, we conducted DEM simulations to study granular flows on erodible surfaces. The energy of the bottom heap-like bed is innovatively controlled. The well-recovered exponentially decaying velocity profile confirms the heap-like features of the bottom bed; such velocity profiles are commonly observed in avalanches and landslides with erodible substrates. Analyzing the bulk behavior of the granular bed, the bottom erodible bed can greatly enhance the mobility of granular flow. With increasing bottom bed energy, the height of the shear band increases and the velocity at the free surface is enhanced. However, the    9. Dimensionless pressure versus dimensionless shear rate for various controlling energy at the shear zone. The inset shows the rheological data for the entire granular bed, the squares are rheological data of vertical chute flow with channel width of 110d and initial volume fraction of 0.56 and 0.59 for blue and red squares, respectively (Baker et al. 2022). The master curve (black line) at steady state is that proposed by Chialvo et al. (2012). rheological behavior at the shear band still belongs to the inertial regime, indicating that the previous constitutive model can also be applied to describe inclined surface flows on heaps. Moreover, granular temperature can be a critical parameter for future development of constitutive model for multi-regime granular flow.
Although only monodisperse spherical particles are used in this work, the dispersity and shape of particles will not affect the effectiveness of proposed configuration modelling granular flow on erodible surfaces. Therefore, our results have several implications with respect to studies of geophysical flows involving entrainment. First, attention should be paid to the effect of entrainment on geophysical flows in hazard zoning and forecasting which are main ways of risk reduction. Entrainment can greatly increase the volume and mobility of geophysical flows, but the current models for practical hazard zoning and forecasting haven't fully considered entrainment yet, which will underestimate runout distances and destructive power. Second, the rheology in the solid-like regime deviates from steady-state rheology obtained using simple shear, which calls for further research to investigate its physical mechanism and to develop suitable constitutive model. The energy-controlled configuration proposed in this work provides an effective and efficient tool to study the complex rheological behaviour of geophysical flows with entrainment and can be extended to scenarios with different dispersity, shape and interstitial fluid involved in real geophysical flows.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
All data used within this publication can be accessed at: https://doi. org/10.7488/ds/3778.