Key Points: 8 • Bed heights of bedload-dominated rivers modeled by Distinct Element Method (DEM) 9 simulations follow a Gaussian distribution. 10 • The standard deviation of bed height, s η , increases as the shear stress increases. 11 • Peak entrainment of bed particles occurs at a distance 2s η above the average bed 12 height. Abstract 14 We investigate the statistics of bed height variability and particle entrainment height un-15 der steady state bedload transport conditions using distinct element method (DEM) sim-16 ulations. We do so in the context of a theoretical probabilistic formulation derived to 17 better capture spatial variation in sediment exchange between bed material load and al-18 luvial deposits (Parker et al., 2000). Using DEM simulations, we set the foundation for 19 a physics-based closure of this probabilistic framework toward its practical implemen-20 tation. Towards this, we perform DEM simulations for bedload transport under simi-21 lar boundary conditions to those of Wong et al. (2007) laboratory experiments: a bed 22 of gravel particles of median grain size 7.1mm with lognormal grain size distribution trans-23 ported under bed shear stresses ranging from τ 0 = 8.70 to 13.7 Pa. We first validate 24 these simulations by demonstrating that they capture measurable transport and height 25 variations from experimental measurements. We then compute the statistics of both the 26 bed height and entrainment height as they vary with bed shear stress. We find that vari-27 abilites in both bed height and entrainment height variabilities follow Gaussian distri-28 butions, for which: (1) the standard deviation of bed height variability s η increases with 29 shear stress, and (2) the peak entrainment height occurs a distance of twice the stan-30 dard deviation of bed height variability (2s η) above the mean bed height. We discuss 31 implications of these results and next steps for understanding these transport statistics 32 under a broader range of conditions. 33