1D grain-size-specific morphodynamic model for sand bed rivers

Morphodynamic model research code for sand bed rivers used in Leonard and Wilcock (2024). Morphodynamic models evolve the bed surface grain size and topography from initial conditions to a steady state using conservation of mass and momentum for open channel flow and sediment mass conservation. We use a simplified version of channel geometry with constant width and no floodplain such that the only adjustment is in bed texture and bed aggradation or degradation. The model starts with a specified slope and bed grain size at each model node and calculates bed shear stress from a specified discharge to predict the transport rate and grain size. Bed topography and grain size adjust based on the difference between the mass of each sediment size fraction delivered to and transported from each node at each timestep. The 1D shallow water equation of mass and momentum conservation describes flow within the rectangular channel. Flow resistance and the skin friction portion of the total boundary stress are specified using the Wright and Parker (2004b) formulation, which accounts for the effects of density stratification and flow resistance over dunes. The grain-size-specific formulation of the Exner equation (Parker et al., 2007), which conserves sediment mass for individual size fractions, determines bed grain size and topography. The channel bed is divided into an upper active layer that exchanges with the bed material load, a lower substrate layer that maintains a constant grain size, and an interface layer that exchanges sediment between the active layer and the substrate as the bed aggrades and degrades (Hirano, 1971). The active layer thickness is specified as the height of the bedforms, predicted as a function of flow depth using the relation of Julien and Klaassen (1995). The grain size of the interface layer evolves as the bed aggrades and erodes using the relation formulated by Hoey and Ferguson (1994) and Toro-Escobar et al. (1996). Grain-size specific volumetric bed material transport rates are calculated using a separate transport relation for bed and suspended loads. Total volumetric bed material transport is the sum of transport of each grain size fraction in the bed and suspended loads. We use the Wright and Parker (2004b) entrainment model (W-P) coupled with a Rouse profile and van Rijn (1984) initiation of suspension criterion to estimate suspended load transport. W-P is a modified version of the Garcia and Parker (1991) (G-P) entrainment model that accounts for reduced mixing due to density stratification in the presence of large suspended loads. G-P and W-P are the only entrainment models with a mixed-size hiding function tested against field data, making this relation ideal for predicting size-selective transport in the suspended load that drives the sorting of bed grain. Bed load is calculated from the Ashida and Michiue (1972) relation (A-M), which includes the Egiazaroff (1965) hiding function. A-M was developed from flume measurements of sand bedload, making this relation ideal for our modeling purpose.

本代码为用于Leonard与Wilcock（2024）研究的沙质河流地貌动力学模型（Morphodynamic model）研究代码。地貌动力学模型基于明渠水流质量与动量守恒以及泥沙质量守恒原理，将床面粒度与地形从初始状态演化至稳态。本研究采用简化渠道几何构型，即宽度恒定且无漫滩，仅床面质地与床面淤积或冲刷可发生调整。模型以每个计算节点处指定的坡度与床面粒度为初始条件，通过指定流量计算床面剪切应力，进而预测输沙率与粒度。床面地形与粒度随每个时间步内各泥沙粒级在节点处的输入与输出质量差进行调整。一维浅水质量与动量守恒方程（1D shallow water equation）描述矩形渠道内的水流运动。水流阻力与总边界应力中的肤面摩擦（skin friction）项采用Wright与Parker（2004b）公式确定，该公式考虑了密度分层与沙丘上水流阻力的影响。针对单个泥沙粒级守恒泥沙质量的埃克纳方程（Exner equation，Parker等，2007）的粒级特异性公式，用于确定床面粒度与地形。渠道床体分为三层：与床沙载荷（bed material load）交换的上部活动层（active layer）、维持恒定粒度的下部基质层（substrate layer），以及随床面淤积与冲刷在活动层与基质间进行泥沙交换的界面层（interface layer）（Hirano，1971）。活动层厚度设定为床面形态高度，该高度通过Julien与Klaassen（1995）提出的关系式，以水流深度为变量进行预测。界面层的粒度随床面淤积与冲刷，采用Hoey与Ferguson（1994）以及Toro-Escobar等（1996）提出的关系式进行演化。针对各粒级的床沙体积输沙率，分别采用床沙载荷与悬移载荷（suspended load）的输沙关系式进行计算。总床沙体积输沙量为各粒级床沙与悬沙输沙量之和。本研究采用Wright与Parker（2004b）的起动模型（entrainment model，简称W-P），结合劳斯剖面（Rouse profile）与van Rijn（1984）的悬移起动准则，估算悬沙输运量。W-P模型是Garcia与Parker（1991，即G-P模型）起动模型的改进版本，考虑了高悬沙浓度下密度分层导致的混合作用减弱。G-P与W-P模型是仅有的两种采用混合粒级掩蔽函数（hiding function）且经野外观测数据验证的起动模型，因此该关系式非常适合用于预测驱动床面粒度分选的悬沙粒级选择性输运。床沙载荷采用Ashida与Michiue（1972）的关系式（简称A-M）计算，该公式包含了Egiazaroff（1965）的掩蔽函数。A-M关系式基于水槽沙质床沙载荷实测数据推导得出，因此非常适配本研究的建模需求。