第四纪研究  2017, Vol.37 Issue (2): 368-379   PDF    
沙丘形态动力学数值模拟:以真实空间元胞自动机模型(ReSCAL)为例
张德国 , 杨小平②,①     
(① 中国科学院地质与地球物理研究所, 新生代地质与环境重点实验室, 北京 100029;
② 浙江大学地球科学学院, 杭州 310027)
摘要:沙丘研究不再局限于现有沉积状态的描述,而是延伸到其形成发育机理的研究。本文通过ReSCAL沙丘模型的数值模拟分析,获得如下结论:1)ReSCAL沙丘模型实现了风况因素影响下复合型沙丘的形成发育过程,揭示了复合沙丘表面次生沙丘形成与星状沙丘翼角延伸的动力机制,提出了沙丘形态与近地表风况之间的互馈机制对近地表风沙流的影响是次生沙丘出现和星状沙丘翼角延伸的重要原因。2)示踪了沙粒的移动轨迹,发现了沙粒每次在沙丘背风坡的埋藏深度具有随机性,使得沙颗粒的平均埋藏时间不同于以往计算的更新时间,继而发现新月形沙丘沙平均埋藏时间等于纵向最大剖面面积与其对应单宽输沙率的比值。3)基于沙丘脊线走向数学方法的发展,并结合ReSCAL沙丘模型的数值模拟结果,发现双向风况模式下,沙源供应能力的变化不仅能塑造不同的沙丘形态,而且能改变沙丘的生长发育模式,甚至是沙丘臂膀的延伸方向(走向)。
主题词沙丘形态动力学     ReSCAL沙丘模型     次生沙丘     星状沙丘     平均埋藏时间     风况     沙源供应能力     沙丘脊线走向    
中图分类号     P534.63+2;P512.2+1;P931.3                     文献标识码    A

1 引言

近年来,沙丘数值模拟方法的发展为沙漠研究提供了新机遇。沙丘不仅是现代气候环境变化的信号载体[1],也是研究早期地球气候环境变化的重要素材[2]。在风力、水力[3]、下附地形[4]等外部因素的影响下,沙丘的形态、形体、空间位置、内部沉积结构和内部沉积物的埋藏时间均会发生不同程度的后期改造。因此,沙丘为地貌学和地质学的研究提供了宝贵的资料。沙丘外部形态和内部沉积环境的变化及其与环境变化之间的相互联系,引起了深入的讨论和研究[5~8],不过对其形成发育的机理认识有限。沙丘的形成发育是一个复杂、多时空尺度的自组织过程[9]。野外观测[10, 11]、实验室模拟[12, 13]等重要和常见的研究手段在探究沙丘形成发育机理时,通常受时空尺度、初始条件设定等因素的限制,较难获得全面和规律性的认识。计算机数值模拟方法的快速发展不仅突破研究方法的瓶颈,而且通过模型模拟与实验数据结果的对比研究,将更好地帮助我们理解沙丘乃至沙漠的形成发育过程。

沙漠领域中,数值模拟内容包括风沙流[14]、沙波纹[15]和沙丘[16]等。在沙丘动力学模型方面,Parteli等[17]、Duran和Herrmann[18]以及Hersen[19]通过描述沙丘表面流场与输沙率的关系,建立了基于连续介质理论的线性模型(图 1a),此类模型将来流风场和近地表风沙流分别处理为稳定不可压的边界层流动和类流体,合理地模拟了沙丘发育和移动过程,阐明了沙丘特征波长与不同流体密度和沙粒密度之间的关系。但在此类模型中,沙丘的高度和输沙率等重要参数都按连续变量处理,而事实上,风沙系统是一个离散系统,沙粒的运动路径、沙丘尺寸的空间分布及其对应的表面流体运动状态等系统物理变量随时间变化的规律较难用连续函数描述,因此连续模型在解释沙丘形成过程的有效性方面值得商榷。

图 1 沙丘连续介质模型(a)与元胞自动机模型(b)示意图 Fig. 1 Schema of continuous dune model (a) and cellular automaton model (b)

元胞自动机模拟方法(图 1b)却能较好的处理此类问题[20~23],这类离散模型模拟的对象不是庞大的沙颗粒和流体分子个体,而是数目减少的沙颗粒和流体分子“微团”。这样一来,不仅模型计算的效率得到提高,而且沙颗粒和流体分子“微团”在一定的作用规则约束下,能够呈现出与自然界沙丘相似的形态,以及沙漠中沙丘间的融合、分离、再生等过程。但是此类模拟方法存在以下不足之处:首先,模型设计者所假设的元胞相互作用规则不能完全真实地反映微观尺度下的物理动力过程;其次,空间尺度确定的随意性(例如:元胞个体的单位长度,元胞个体的跳跃距离等)并不能反映实际沙漠形成和发育过程。因此,单位时空尺度的确定是迫切需要解决的难题。

本文研究中,将首先介绍一种新的沙丘模型,即真实空间元胞自动机实验室模型(Real-space cellular automaton laboratory,简称ReSCAL)[24]。此模型较好地再现了复合型沙丘表面的次生沙丘形态。通过次生沙丘在模型与自然界的对比研究,确定了模型的单位时空尺度,实现了模型的量化研究。本文继而通过ReSCAL模型的结果,探讨了风况条件下沙丘外部形态发育机理和内部沙粒的平均埋藏时间。此外,在ReSCAL模型模拟应用的基础上,进一步探讨了双风向模式下,沙丘形态及其脊线走向的发育机理。

2 模型

ReSCAL模型是一种混合模型,它由用以模拟沙粒运动过程的元胞自动机模块(图 2)和用以模拟流体运动的格子气自动机模块(图 2)组成。

图 2 ReSCAL沙丘模型组成及其内部关系示意图 Fig. 2 Schema of ReSCAL dune model's components and its relationship
2.1 元胞自动机模块(CA)

在元胞自动机模块中,Narteau等[24]依据风沙系统中的风蚀、搬运、沉积、扩散和因重力效应导致的“沙崩”等物理过程,设计了模型中沙粒“微团”之间的相互作用规则及其相对应的作用强度(Λe、Λt、Λd、Λg和Λc)(图 3a)。每一个“微团”都是一个单位长度为l0的立方体,而每一对相邻“微团”构成双“微团”单元。对于出现在模型中每一个独立的物理过程,都有一组相对应的双“微团”转换过程。图 3a左上角的立体图表示流体流动方向与每一个“微团”单元毗邻的6个相邻“微团”单元。

图 3 沙丘模型作用规则示意图 (a)元胞自动机模块相互作用关系图;(b,c)格子气自动机模块中流体粒子作用规则包含:移动(b)与碰撞(c)两个阶段,其中,箭头代表流体粒子的速度向量。每一个元胞单元中存在一个黑色节点,流体粒子只能在节点与节点之间的连线上运动 Fig. 3 Transition rules of sand dune model (a)Transition rules of cellular automaton; (b, c)Transition rules of fluid particles in LGCA model include migration (b) and collision (c). Particles and their velocity vectors are represented by arrows. Each dot represents the nodes of the LGCA and the center of the cells of the model of sediment transport, fluid particles move only along the connections between these nodes

关于沙堆的临界状态,Bak[25]有这样的描述:“只要有一粒沙滚落,都无法预料会出现什么样的结果,也许什么都不会发生,也许只有很少的沙粒会滑落,也许一个很小面积的沙粒滑落正好导致一场连锁反应”。在Bak的沙堆晶格模型理论之上,我们重新设定了沙粒“微团”的崩塌休止角,即沙丘的演化过程中,沙丘迎风坡与背风坡坡度分别不能大于20°和33°,当堆积角大于其对应的休止角时,沙粒“微团”因其自身重力发生崩塌,崩塌的沙粒“微团”会随机流向相邻的沙粒“微团”。以此类推,就发生了类似沙粒滑落时引起的自组织崩塌现象。“自组织”是指该状态的形成主要是由系统内部组份间的相互作用产生,而不是由任何外界因素控制或主导所致。

2.2 格子气自动机模块(LGCA)

在离散类模型中,Werner[9]、Nishimori和Ouchi[20]以及Baas和Nield[23]的模型多针对风沙二相流中固体颗粒的运动过程进行了模拟,却忽略了其中流体运动状态的模拟。在ReSCAL沙丘模型的设计之初,我们引入了格子气自动机模块(LGCA),用于流体运动的模拟。LGCA模块中,每个“微团”中心存在一个格点,流体粒子只能在格点之间的格线上以相同的离散时间步长运动,流体粒子与模型上壁面之间的碰撞遵循法向反弹法则,与沙粒“微团”、下壁面之间的碰撞遵循镜面反弹法则,与其他流体粒子之间遵循散射法则。当不同流体粒子在同一格点相遇时,根据散射规则,进行相互碰撞而交换速度。碰撞过程中粒子之间并不相容,即在同一时刻同一网格点上,每一个速度方向最多允许有一个流体粒子。以上模型的设计,不仅保证了整个流场的动量、能量守恒,而且较好的实现了流体与沙丘形态之间的相互作用关系。流体模型中获得的近地表风速,被用于计算近地表风蚀应力,同时,在临界起沙应力的约束下,我们建立了一个简单的,具有约束条件的线性风蚀起沙方程。最终较真实的实现了风沙二相流的数值模拟。

流体粒子在每一个时间步长的演化包括两部分:1) 移动(图 3b),粒子沿移动方向向距离最近的节点运动;2) 碰撞(图 3c),当不同的粒子同时到达某个节点时,按照一定的碰撞规则发生碰撞并改变运动方向。格子气模块具有两重意义:1) 建立一个尽可能简单的模型,能够用来模拟一个由大量粒子组成的系统;2) 反映粒子真实碰撞的本质,经较长时间模拟后可以获得流体的宏观特性。

LGCA模块从微观的尺度上真实的反映了流体的运动状况,例如:沙丘迎风坡上流体速度的加速。另外,临界摩阻风速uc概念的引入,将CA模块与LGCA模块很好地耦合在一起:当实际摩阻风速u*i大于临界摩阻风速uc时,沙粒开始运动,随着沙粒总量的增加,最终会达到一个平衡状态,即饱和输沙率(Qsat),而达到饱和输沙率所需的路径长度为lsat;当u*i < uc时,沙粒不运动[26]

2.3 沙丘模型单位时空尺度

Narteau等[24]成功模拟了沙漠中各类沙丘的形成发育过程,再现了沙丘表面次生沙丘的形成发育过程(图 4a),这为研究者认识沙丘发育过程的空间多尺度性以及时空尺度的量化奠定了基础。沙丘表面的次生沙丘的特征波长λ=20m[27](图 4b)。通过与模型中次生沙丘特征波长λmax=40l0的比较(图 4c),确定了模型中沙粒“微团”的长度l0=0.5m。在模型单位长度确定的基础之上,我们通过模型与实验中测量的平沙地水平饱和输沙率的比对,确定了模型的单位时间尺度t0的计算公式如下[24]

(1)
图 4 ReSCAL沙丘模型时空尺度的确定 (a)新月形沙丘形态及其表面气流状态;(b)沙漠中次生沙丘波长个数统计图[27];(c)不稳定波长随风速与临界起沙风速比值的变化关系[24] Fig. 4 Setting the length and time scales of ReSCAL dune model (a)Morphology and velocity streamlines of a barchan dune; (b)Histograms of the wavelength measured of superimposed dune[27]; (c)Evolution of the most unstable wavelength with respect to the ratio between the flow velocity and the threshold velocity for sediment transport[24]

公式(1) 中Qsat(τ1)是模型中计算得到的单宽饱和输沙率,〈Qsat〉是由实测风况数据估算预测的平均单宽饱和输沙率。数学公式如下:

(2)

公式(2) 中,其中u*i分别代表不同时间t1titNi∈[1, N]的风速和风向数据。另外,根据Iversen和Rasmussen[28]的公式,可以估算Qsat(u*i),公式如下:

(3)

公式(3) 中,ρfρs分别代表空气和沙粒的密度,ρf=1.205kg/m3ρs=2650kg/m3u*uc分别指示摩阻风速和临界摩阻风速,d是沙粒直径,d=180μm,g重力加速度,g=9.8m/s2。临界摩阻风速uc可依据Iversen等[29],由公式(4) 得到:

(4)
3 结果 3.1 单一风况条件下复合新月形沙丘

从沙丘形态上来看,沙丘形态的常见类型为新月形、横向、纵向和星形(金字塔形)[30]。从沙丘尺寸上看,沙丘的波长从几十米到几千米不等[31]。不同尺寸的沙丘在相同外部条件的作用下,会以不同的速率运动,继而会呈现出由于沙丘之间融合、分离等过程之后的复杂、复合沙丘形态。复合新月形沙丘的平面形态具有新月的外形,新月的两个尖端指向下风向,形体高大,在巨大沙丘的迎风坡上,层层叠置着次生沙丘,形成了“双层结构”;沙丘复合体的排列方向和优势风向大致相垂直[32]

科学家对复合沙丘形态的研究主要是基于数值模拟和遥感影像、野外实测。其中,数值模拟方面,Zhang等[33]利用ReSCAL沙丘模型模拟了复合新月形和横向沙丘的形成发育过程。单一风况下风速的增大,加强了沙丘表面的风蚀过程,这一过程的增强足以克服因重力和水平输沙作用引起的沙颗粒向下和四周的运动,从而约束了沙丘表面次生沙丘的形成发育,使沙丘形态沿主导风向和垂直向上的方向发展,即沙丘纵向剖面的坡形比值和移动速率增大,次生沙丘的数量减少。相反,即沙丘坡形比值和移动速率减小,次生沙丘数量增大。

进一步研究显示沙丘的形态也能反作用于流体的运动状态。具体来讲,沙丘迎风坡坡度的增加迫使风速加速,继而导致沙丘顶端与底部的输沙率比值呈线性增长关系,这与Lancaster[34]的野外观测以及Jackson和Hunt[35]的数学理论模型预测结果一致。然而,他们却忽略了沙丘沿主导风向移动时,发育于沙丘表面的次生沙丘的形成发育机制和移动速率的研究(图 5a5b)。Zhang等[33]的模拟研究结果显示,发育于沙丘表面的次生沙丘与发育于平坦沙质床面沙丘的形成发育机制是一致的,次生沙丘的特征波长以及分布的密度均受控于沙丘表面流体流速与携沙能力的大小。尽管次生沙丘与主体沙丘的形体大小、移动速率各异(图 5c),但是它们的移动速率均能被其顶端的输沙率归一化到同一条曲线之上(图 5d),说明它们遵循着同样的物理动力过程。其关系如公式(5) 所示[33]

(5)
图 5 单一风况下次生沙丘的形成发育机理 (a,b)复合型新月形沙丘的数字模拟效果图[35]与卡塔尔卫星影像图(24°59′43.20″N,51°22′18.63″E);(c)沙丘与次生沙丘移动速率与其高度关系图[33];(d)3种不同风力条件下沙丘及其次生沙丘的移动速率依据其对应的沙丘顶端输沙率归一化关系[33] Fig. 5 Dune morphodynamic under unidirectional wind regime (a, b)Observations of superimposed dune patterns in model[35] and nature of Qatar(24°59′43.20″N, 51°22′18.63″E) for megabarchan; (c)The propagation speed of dunes and superimposed dunes as function of theirs heights[33]; (d)The propagation speed of dunes and superimposed bedforms rescaled with respect to the sediment flux at the crest for three different flow strengths[33]

其中,c是沙丘移动速率,a是依赖于沙丘的长高比的常数,Qsatflat是沙丘上风向平沙地测量的单宽饱和输沙率,H是沙丘高度,H0对应的是最小新月形沙丘的临界高度。

3.2 单一风况条件下新月形沙丘沙的埋藏时间

沙丘内部的风沙沉积物蕴含重要的年代信息,是认识沙丘乃至沙漠形成、演化和沙漠地区环境气候变迁的“时间标尺”[36~38]。然而,由于沙丘运动过程中沙颗粒的移动路径、被埋藏过程和在沙丘内部停留时间等观测数据难以获取,以及沙丘外部形态受到外部营力的破坏和改造等原因,确定沙粒在沙丘内部的平均停留时间是极其困难的。一方面,沙颗粒的运动轨迹是一个随机的过程,每个沙颗粒的停留时间是不一致的;另一方面,现有的风成沙沉积物光释光年代数据只能代表沙颗粒距离上次曝光后的埋藏时间,不能代表其在进入沙丘内部后经历多次暴露-埋藏循环过程后,最终离开沙丘所需的停留时间。因此,通过数值模拟研究沙颗粒在沙丘内部的停留时间和移动轨迹,是认识沙丘沙平均埋藏时间和帮助理解地层中风成沙沉积物光释光年代数据的前提和基础。

已有的研究[39]认为沙丘移动速率与其形态之间的关系可被表示为:

(6)

公式(6) 中Qcrest为沙丘顶端单宽输沙率,H为沙丘高度,c为移动速率。假设移动过程中新月形沙丘的形态处于恒定状态,其更新时间Td近似等于沙丘自身长度L与其移动速率c的比值[40],数学公式如下:

(7)

此外,沙丘更新时间也是沙丘内部所有沙粒运动一次所需的时间。由公式(7) 得知,沙丘的移动速率决定了沙丘所有沙粒循环一次所需时间的长短。沙丘的移动速率与其表面的风沙输沙程度是呈线性增长的关系。换言之,沙丘更新时间受控于沙丘表面风沙活动的剧烈程度。然而,沙丘表面的输沙过程受到多个外部因素(水平输沙过程、风况和沙丘迎风坡坡形比等)的影响[41],因此弄清沙丘移动过程中,其表面输沙过程与沙丘形态之间的关系对于解析沙丘更新时间是至关重要的。Andreotti等[40]在考虑沙丘更新时间Td时,沙粒的运动轨迹被假设限定在一个二维纵向水平面上,而实际沙粒受到气流、颗粒碰撞、重力和沙崩等作用后,移动轨迹会在水平方向上发生偏移,因此会导致实际更新时间要不同于数学公式估算的更新时间。另外,在水平输沙过程的作用下,沙粒将经历多次的暴露-埋藏循环过程,最终经由新月形沙丘的翼角离开沙丘主体[33]。那么,沙丘更新时间Td与沙丘的平均埋藏(停留)时间TR之间是否是对等关系?沙丘表面风沙过程是如何影响沙粒的平均埋藏时间?这些问题的解答需要进一步研究来回答。

为了研究沙丘平均埋藏时间,Zhang等[42]将沙丘的颗粒“微团”全部标记为黑色,黑色“微团”将由模型的下游边界流出,然后同等数量的白色“微团”将从沙丘的上游进入沙丘(图 6a),随时间推移,将近一半黑色沙粒离开沙丘主体所需的时间即是沙丘平均埋藏时间的测量值。图 6b显示了被标记沙粒数量随时间推移成指数递减的趋势图,其拟合计算公式如下:

(8)
图 6 新月形沙丘沙平均埋藏时间[42] (a)被标记与未被标记的沙粒分别以黑色和白色来区分;(b)被标记黑色沙颗粒随时间变化图;(c)沙丘平均埋藏时间与其更新时间关系图;(d)沙丘沙平均埋藏时间与V/(QinW)比值的线性关系图;(e)沙丘沙平均埋藏时间与S/Qin比值的线性关系图;(f)新月形沙丘移动过程中,被标记沙粒在垂直、水平空间上的移动路径 Fig. 6 Mean residence time of sedimentary cells in a barchan[42] (a)The tagged and untagged sedimentary cells are shown in black and white, respectively; (b)Number of tagged sedimentary cells decrease with respect to time; (c)The mean residence time with respect to the renewal time of the dune; (d)The mean residence time with respect to V/(QinW); (e)The mean residence time with respect to S/Qin; (f)Horizontal and vertical motions of sedimentary cells during the propagation of the steady state barchan

公式(8) 中,TR为沙丘的平均埋藏时间。为了更好的了解平均埋藏时间与沙丘形态之间的关系,我们模拟计算了沙丘在不同个体大小、不同水平输沙率、不同风力强弱条件下的平均埋藏时间TR。由图 6c所见,平均埋藏时间TR与更新时间Td之间呈线性增长关系,当风力强弱增大时,二者之间数值逐渐靠近,但TR仍远大于Td,表明外界因素改变引起的沙丘表面输沙过程的变化,是导致沙丘平均埋藏时间变化的因素之一,而输沙过程变化引起的沙丘形态的变化或许是导致平均埋藏时间变化的主要因素。从沙丘的形态来看,处于稳定状态的新月形沙丘可被认为是一个与外界环境进行沙物质交换的“容器”,那么沙颗粒的埋藏时间应等于“容器”的体积V与其输沙量QinW的比值(图 6d)。图 6d显示它们之间的线性关系,但其系数值大于1,这说明新月形沙丘并不是一个混合均匀的沙质物体。通过进一步的研究,Zhang等[42]提出平均埋藏时间等于沙丘最大纵向剖面的体积S与单宽输沙率Qin之间的比值(图 6e)。因此,沙丘最大纵向剖面的更新时间可能控制着整个沙丘系统的平均埋藏时间。其关系如公式(9):

(9)

以沙粒“微团”个体为研究对象,研究了沙粒“微团”个体在沙丘内部的停留时间及其空间位置分布状态,追踪了每个沙粒“微团”的运动轨迹(图 6f)。图 6f显示,沙粒“微团”在沙丘落沙坡的埋藏时间是随机分布的。随着时间的推移,沙粒“微团”的水平运动轨迹逐渐向沙丘的两翼角方向偏移,说明沙丘表面的水平输沙过程不仅能“平滑”沙丘脊线的弯曲度和“抑制”次生沙丘的形成[33],而且还影响了沙颗粒在沙丘内部的停留时间。

3.3 复杂风况条件下星状沙丘的形成发育机理

在临界起沙风的作用下,风成沙丘的形态将受到其表面气流的改造[43]。这种持续不断的再改造过程导致星状沙丘(图 7a)的形成。星状沙丘的形成可能是两个或是多个沙丘融合的结果(图 7c)[44~47],也可能是单个沙丘发育出臂膀的结果(图 7b7d)[48]。无论是哪种情况,复杂风况和次生气流在星状沙丘的形成过程中都起到至关重要的作用,但至今难以实现这两种影响因素的量化研究。

图 7 复杂风况下星状沙丘形成发育及其臂膀延伸的动力机制 (a,b)星状沙丘的卫星影像图(中国,44°38′01.72″N,91°36′06.27″E)与数字模拟效果图;(c)星状沙丘形成发育示意图(修改自Lancaster[44]);(d)星状沙丘臂膀延伸动力机制[48];(e)沙丘臂膀形态风况交替周期关系图[48];(f)沙丘臂膀延伸速率与风况交替周期关系图[48] Fig. 7 Star dune morphodynamic and its arms extension mechanism (a, b)Star dune in nature (China, 44°38′01.72″N, 91°36′06.27″E) and model; (c)Formation of star dunes(after Lancaster[44]); (d)Formation and evolution of star dunes using multidirectional wind regimes[48]; (e)Width and height of star-dune arms with respect to Tθ[48]; (f)Arm growth rate with respect to Tθ[48]

Lancaster[44]提出了星状沙丘臂膀形成发育过程(图 7c),翻越沙丘垄脊的初始气流在沙丘的背风侧产生二次流,而在二次流产生的同时,其流动方向也发生了偏移,偏移之后的方向几乎与原来的初始气流方向垂直,二次流可能导致了星状沙丘臂膀的产生,但是它对于沙丘表面风沙活动的贡献是否足以造成沙丘臂膀的产生和延伸,至今没有一个初步的估计。另外,星状沙丘是沙漠环境中形体相对较大的沙丘,它的形态能够记录和综合反映长时间作用于此区域的风况信息。因此,在空间尺度上,星状沙丘往往呈现为一个多空间尺度的沙丘[13, 27],变化范围从其表面上发育的次生沙丘(几十米)到其本身形体的大小(几千米)不等。

在前人研究基础上,Zhang等[48]将星状沙丘的研究重点放在沙丘臂膀的形成发育机理之上(图 7b)。通过对风向变化频率和角度的设定,模拟呈现出具有不同臂膀个数的星状沙丘形态(图 7d)。值得注意的是,不是所有的情况下都能观察到沙丘臂膀的形成发育过程。这种形态上的差异是由不同条件下的风况决定的,但是理论上来说,多风向的风况致使沙颗粒向沙丘的中心位置堆积,造成合成输沙量近似为零,这与星状沙丘臂膀向四周延伸的现象并不吻合。因此,形态对其表面流场的影响是不容忽视的,而前人的研究并未考虑沙丘形态对风速的加速效应,继而造成对输沙量的低估。因此,Zhang等[48]重新考虑了沙丘臂膀表面输沙量的数学公式,将因加速效应产生的输沙量(Q)考虑在内。当任意一方向的风经过沙丘臂膀时,产生的输沙量可分解为垂直于脊线方向的输沙量Q和平行于脊线方向的输沙量QQ是导致臂膀延伸的主要分量,其计算公式如下[48]

(10)

公式(10) 中,Qa是与初始气流相关的输沙量,Qb是与迎风坡坡形比相关的输沙量;HW分别为沙丘臂膀的高度与宽度。为了检验数值模型的准确性,我们计算了不同风况条件下Q在[0°,360°]之间的变化情况,结果显示当风向周期变化的个数为奇数时,Q最大值所对应的角度与我们模型模拟的沙丘臂膀延伸的角度完全一致,而当个数为偶数时,Q的值为0,与我们模拟的结果也是相对应的。在此基础之上,我们进一步改进了沙丘臂膀顶端合成输沙量分量Q的数学计算公式,公式如下[48]

(11)

公式(11) 中,Q1>0是与臂膀坡形比相关的恒定输沙量,Tθ是不同方向的风作用于臂膀的总周期时间。

据公式(11) 计算得知,当风向在一个周期内的个数为奇数时,臂膀会发生延伸,并且臂膀延伸的方向正好与输沙量Q最大值对应的角度一致,角度方程式可表达为:αi= π/n+2/n。当风向个数为偶数时,臂膀不发生延伸,同时Q=0(如图 7d)。通过模拟研究结果显示,因不同来流方向对应的沙丘臂膀坡形比的不同而造成了对不同方向风速的加速效应是有差异的,而这种差异造成了即便是在同样风力大小的条件下,最大合成输沙分量(Q)也不为零,因此沙丘臂膀沿着最大合成输沙分量的方向延伸和发育。进一步的模拟研究显示,沙丘臂膀的形态(图 7e)和延伸速率随着各方向风平均作用时间的增加而系统变化。随着各方向风平均作用时间的拉长,多风向在沙丘表面的合成效果减小,而单个方向风的作用效果逐渐增强。这表示,沙丘臂膀延伸的可能性逐渐减小,甚至消失(图 7f)。

3.4 沙丘形态与近地表风况之间的量化分析

风是塑造沙丘形态的主要因素,换言之,沙丘的形态能够记录基地表风况的变化的历史。因此,量化沙丘形态与近地表风况之间的关系成为研究沙漠地区近地表风场变化历史的关键。Fryberger[49]以沙丘附近气象站点风况资料为基础,计算了合成潜在输沙势(RDP)与总输沙势(DP)之间的比值,并以此为依据局,半定量的确定了风况与沙丘形态之间的对应关系。在此基础之上,Lancaster等[50]利用西沙哈拉沙漠纵向沙丘的脊线走向,并结合沙丘沙的光释光年代数据,重建了该区域的近地表主导风况的变化历史。他们的研究揭示,西沙哈拉沙漠在距今2.5~1.5万年、1.3~1万年和5千年3个不同的时间段,近地表主导风向发生了显著的变化,从而造成了形体大小不一、沙丘脊线走向不同的复合型纵向沙丘的形成。然而,该研究仅利用近地表合成主导风向与纵向沙丘脊线相平行的对应关系,简单推论了地质时期该地区的主导风向。另外,Fryberger[49]的分类方法更多强调的是合成风况与沙丘外部形态之间的对应关系,并未完全实现风况与沙丘脊线走向之间的定量重建,因此并不能真实的反映地质时期近地表主导风况的变化历史。

Hunter等[51]从动力学的角度出发,依据沙丘脊线走向与合成输沙方向之间的关系,重新定义沙丘的分类,将沙丘形态依次分为横向型、斜向型和纵向型(图 8a),逐渐将研究重点由沙丘形态转向沙丘脊线的走向。基于这些认识,Rubin和Hunter[52]对沙丘脊线走向与合成输沙方向之间的关系进行了数学计算,提出了最大总输沙量(Gross Bedform-normal Transport)方法(图 8b8c),半定量的确定了沙丘形态、脊线走向与近地表主导风向之间的关系。随后,Werner和Kocurek[21]进一步优化了Rubin数学模型中斜向沙丘与横向沙丘的过度区域(图 8d)。Fenton等[53]通过美国Great sand dune沙漠中具有不同形体和形态沙丘脊线走向的统计分析,结合Rubin和Hunter[52]的方法,恢复了该区域的近地表风况,同时与现代风况资料计算所得主导风况进行了比较,结果显示二者的方向是近乎一致的。另外,此方法不仅被应用于地球近地表风场的重建,而且还成功应用于其它星体近地表风场的重建[54]

图 8 (a)沙丘形态类型与合成输沙量方向之间的关系[51];(b)沙丘脊线走向与两个输沙量(Q1Q2)之间的关系[52];(c)基于风向夹角(θ)与输沙量比例(N)的沙丘脊线走向预测图,修改自Rubin和Hunter[52];(d)基于风向夹角(θ)与输沙量比例(N)的沙丘脊线走向预测图,修改自Werner和Kocurek[21],其中,θ=α1+α2 Fig. 8 (a)Classification of linear dunes under bimodal wind regimes[51]; (b)The relationship between bedform trend and two incident sand-moving winds with transport Q1 and Q2[52]; (c)Dune orientations within the parameter space[N, θ] of bimodal wind regimes as defined by Rubin and Hunter[52]; (d)Dune orientations within the parameter space[N, θ] of bimodal wind regimes as defined by Werner and Kocurek[21], θ=α1+α2

尽管Rubin和Hunter[52]的数学计算方法得到了很好的检验,然而并未考虑沙丘形成后,其形态对风速的加速效应,继而造成对输沙预测量的低估。最近,Courrech du Pont等[55]等提出的数学模型考虑了沙丘形态对风速的加速效应,改进了沙丘脊线走向的数学模型,并验证了数学模型所预测的沙丘走向与实际沙丘走向一致。沙丘脊线走向数学模型的成熟发展为沙丘模型的应用研究以及沙丘表面主导风况的重建研究提供了理论基础。

在Courrech du Pont数学方法基础上,Gao等[56]利用ReSCAL沙丘模型,不仅模拟总结了风况和沙源双因素影响下的沙丘脊线走向,而且他们的研究结果验证并推进了前人的理论研究成果。结果表明,在双向风况模式下,沙源供应能力的变化不仅能塑造不同的沙丘形态,而且能驱使沙丘脊线按照“Bed instability mode”和“fingering mode”两种不同的模式延伸和发育。另外,Lucas等[57]同样利用ReSCAL模型,模拟验证了障碍物线性沙丘形成发育的机理。早期的研究认为,在线性沙丘的发育过程中,包括延伸和侧向加积两个过程,这也是线性沙丘脊线通常以蜿蜒状向前延伸的原因。然而,在尼日尔的泰内雷沙漠中,发育于山体背风侧的线性沙丘脊线,却以一种近视笔直的状态沿着合成输沙方向延伸。Lucas等[57]的研究成果解释认为这种障碍物沙丘的形成过程中并没有侧向加积作用的参与,障碍物背风侧的阴影区成为了沙颗粒的沉积区域,为线性沙丘的形成和延伸提供了源源不断的沙源。同样在ReSCAL沙丘模型模拟结果的帮助下,Lv等[58]提出了不对称新月形沙丘的形成与发育机理,显示了长期以来造成Bagnold和Tsoar对此类沙丘形成发育产生不同认识的关键点。研究认为准双向风况和有限的沙源条件下发育了不对称新月形沙丘,当双向风况夹角大于临界角度时,将驱使沙丘按照Tsoar所描述的方式发育,反之则按照Bagnold所描述的方式发育。以上研究的显示,风况虽然是控制沙丘形态和走向的关键因素,但沙源供应程度也具有重要意义。

4 总结

沙丘形态的深入研究是对认识沙漠地区环境演化过程具有重要意义。过去几十年在风沙地貌、沙漠环境演化领域有了显著进展,但是对沙丘形成发育机理的认识仍很不足。本文回顾了ReSCAL沙丘模型推动沙丘动力学发展的过程,主要结论如下:

(1) ReSCAL模型较好地耦合了用于模拟沙颗粒运动过程的元胞自动机和用于模拟流体运动的格子气自动机,再现了发育于沙丘表面的次生沙丘形态,并实现了模型的量化分析。

(2) 风况是塑造复合型新月形、横向和星状沙丘形态的主要因素。通过ReSCAL模型研究了风况与沙丘形成发育之间的关系,揭示出因近地表风况与沙丘形态之间的互馈机制引起的近地表风沙流的变化,不仅促使了次生沙丘的形成发育和星状沙丘翼角的延伸,而且控制了新月形沙丘移动过程中,沙粒在其内部的平均埋藏(停留)时间。

(3) 纵向(线性)沙丘普遍存在,但对其形成机理的认识还很有限。在沙丘脊线走向的数学方法发展基础之上,结合ReSCAL沙丘模型的应用研究,本文揭示了沙源供应程度在时空上的变化,不仅能够重塑沙丘的形态,而且能够影响到沙丘的形成发育动力机制,甚至沙丘脊线的走向。

致谢 感谢审稿专家提出了关键性的意见和建议,感谢程宏老师建设性的修改意见,感谢杨美芳编辑对本文的细心审阅。

参考文献(References)
1
朱震达, 吴正, 刘恕等. 中国沙漠概论(修订版). 北京: 科学出版社, 1980, 36-46.
Zhu Zhenda, Wu Zheng, Liu Shu et al. An Outline of Chinese Deserts. Beijing: Science Press, 1980, 36-46.
2
Lancaster N. Geomorphology of Desert Dunes. New York: Routledge, 1995, 190-212.
3
Yang X, Forman S, Hu F et al. Initial insights into the age and origin of the Kubuqi sand sea of Northern China. Geomorphology, 2016, 259: 30-39. DOI:10.1016/j.geomorph.2016.02.004
4
Yang X, Scuderi L, Liu T et al. Formation of the highest sand dunes on Earth. Geomorphology, 2011, 135(1-2): 108-116. DOI:10.1016/j.geomorph.2011.08.008
5
董光荣, 靳鹤龄, 陈惠忠. 末次间冰期以来沙漠-黄土边界带移动与气候变化. 第四纪研究, 1997(2): 158-167.
Dong Guangrong, Jin Heling, Chen Huizhong. Desert-loess boundary belt shift and climatic change since the last interglacial period. Quaternary Sciences, 1997(2): 158-167.
6
杨小平, 师长兴, 李炳元等. 从地球系统科学角度浅析中国地貌若干问题研究的新进展. 第四纪研究, 2008, 28(4): 521-534.
Yang Xiaoping, Shi Changxing, Li Bingyuan et al. Some aspects about Chinese geomorphology: Recent progresses from an earth system science perspective. Quaternary Sciences, 2008, 28(4): 521-534.
7
Yang X, Eitel B. Understanding the interactions between climate change, landscape evolution, surface processes and tectonics in the Earth System: What can the studies of Chinese deserts contribute?. Acta Geological Sinica, 2016, 90(4): 1444-1454. DOI:10.1111/acgs.2016.90.issue-4
8
Yang X, Scuderi L A, Wang X et al. Groundwater sapping as the cause of irreversible desertification of Hunshandake Sandy Lands, Inner Mongolia, Northern China. Proceedings of the National Academy of Sciences of the United States of America, 2015, 112(3): 702-706. DOI:10.1073/pnas.1418090112
9
Werner B T. Eolian dunes: Computer simulations and attractor interpretation. Geology, 1995, 23(23): 1107-1110.
10
王训明, 董治宝, 屈建军等. 塔克拉玛干沙漠简单线性沙丘形态-动力学过程研究. 中国沙漠, 2003, 23(3): 257-262.
Wang Xunming, Dong Zhibao, Qu Jianjun et al. Studies on the morphodynamic processes of simple linear dunes in Taklimakan Desert. Journal of Desert Research, 2003, 23(3): 257-262.
11
Bullard J E, Thomas D S G, Livingstone I et al. Analysis of linear sand dune morphological variability, southwestern Kalahari Desert. Geomorphology, 1995, 11(3): 189-203. DOI:10.1016/0169-555X(94)00061-U
12
凌裕泉, 吴正, 刘绍中等. 新月形沙丘形态的模拟实验研究. 地理科学, 1998, 18(1): 88-93.
Ling Yuquan, Wu Zheng, Liu Shaozhong et al. Simulated study on barchan dune forms. Scientia Geographica Sinica, 1998, 18(1): 88-93.
13
Hersen P, Douady S, Andreotti B. Relevant length scale of barchan dunes. Physical Review Letters, 2002, 89(26): 455-460.
14
Anderson R S, Haff P K. Simulation of eolian saltation. Science, 1988, 241(4867): 820-823. DOI:10.1126/science.241.4867.820
15
Zheng Xiaojing, Bo Tianli, Xie Li. DPTM simulation of aeolian sand ripple. Science China Physics, Mechanics & Astronomy, 2008, 51(3): 328-336.
16
Kroy K, Sauermann G, Herrmann H J. Minimal model for aeolian sand dunes. Physical Review E, 2002, 66(1): 162-173.
17
Parteli E, Duran O, Herrmann H. Minimal size of a barchan dune. Physical Review E, 2007, 75(1): 011301. DOI:10.1103/PhysRevE.75.011301
18
Duran O, Herrmann H. Vegetation against dune mobility. Physical Review Letters, 2006, 97(18): 13831-13840.
19
Hersen P. On the crescentic shape of barchans dunes. The European Physical Journal B, 2004, 37(4): 507-514. DOI:10.1140/epjb/e2004-00087-y
20
Nishimori H, Ouchi N. Computational models for sand ripple and sand dune formation. International Journal of Modern Physics B, 1993, 7(9n10): 2025-2034. DOI:10.1142/S0217979293002742
21
Werner B, Kocurek G. Bedform dynamics: Does the tail wag the dog?. Geology, 1997, 25(9): 727-730.
22
Nishimori H, Yamasaki M, Andersen H. A simple model for the various pattern dynamics of dunes. International Journal Modern Physics B, 1999, 12(12): 257-272.
23
Baas A C W, Nield J M. Modeling vegetated dune landscapes. Geophysical Research Letters, 2007, 34(6): L06405.
24
Narteau C, Zhang D, Rozier O et al. Setting the length and time scales of a cellular automaton dune model from the analysis of superimposed bed forms. Journal of Geophysical Research, 2009, 114(F3): 171-183.
25
Bak P. How Nature Works: The Science of Self-Organised Criticality. New York: Copernicus Press, 1996.
26
Bagnold R A. The Physics of Blown Sand and Desert Dunes. London: Methuen, 1954, 85-94.
27
Elbelrhiti H, Claudin P, Andreotti B. Field evidence for surface-wave-induced instability of sand dunes. Nature, 2005, 437(7059): 720-723. DOI:10.1038/nature04058
28
Iversen J D, Rasmussen K R. Iversen J D, Rasmussen K R. The effect of wind speed and bed slope on sand transport. Sedimentology, 1999, 46(4): 723~731. Sedimentology, 1999, 46(4): 723-731. DOI:10.1046/j.1365-3091.1999.00245.x
29
Iversen J D, Greeley R, Marshall J R et al. Aeolian saltation threshold: The effect of density ratio. Sedimentology, 1987, 34(4): 699-706. DOI:10.1111/sed.1987.34.issue-4
30
Wasson R J, Hyde R. Factors determining desert dune type. Nature, 1983, 304(5924): 337-339. DOI:10.1038/304337a0
31
Lancaster N. Geomorphology of Desert Dunes. New York: Routledge, 1995, 39-40.
32
吴正等. 风沙地貌与治沙工程学. 北京: 科学出版社, 2003, 148.
Wu Zheng et al. Geomorphology of Wind-drift Sands and Their Controlled Engineering. Beijing: Science Press, 2003, 148.
33
Zhang D, Narteau C, Rozier O. Morphodynamics of barchan and transverse dunes using a cellular automaton model. Journal of Geophysical Research, 2010, 115(F3): 1480-1493.
34
Lancaster N. Variations in wind velocity and sand transport on the windward flanks of desert sand dunes. Sedimentology, 1985, 32(3): 581-593.
35
Jackson P S, Hunt J C R. Turbulent wind flow over a low hill. Quarterly Journal of the Royal Meteorological Society, 1975, 101(430): 929-955. DOI:10.1002/(ISSN)1477-870X
36
Vermeesch P, Fenton C R, Kober F et al. One million year residence time of Namib dune sand from cosmogenic nuclides. Nature Geoscience, 2010, 3(12): 862-865. DOI:10.1038/ngeo985
37
Lancaster N. Desert dune dynamics and development: Insights from luminescence dating. Boreas, 2008, 37(4): 559-573. DOI:10.1111/bor.2008.37.issue-4
38
Yang X, Li H, Conacher A. Large-scale controls on the development of sand seas in Northern China. Quaternary International, 2012, 250(2): 74-83.
39
Claudin P, Andreotti B. A scaling law for aeolian dunes on Mars, Venus, Earth, and for subaqueous ripples. Earth Planetary Science Letters, 2006, 252(1-2): 30-44. DOI:10.1016/j.epsl.2006.09.004
40
Andreotti B, Claudin P, Douady S. Selection of dune shapes and velocities. Part 1: Dynamics of sand, wind and barchans. European Physical Journal B, 2002, 28(3): 321-339. DOI:10.1140/epjb/e2002-00236-4
41
Thomas D S G. Arid Zone Geomorphology: Process, Form and Change in Drylands(Third edition). Chichester: John Wiley & Sons, Ltd, 2011, 455-479.
42
Zhang D, Yang X, Rozier O et al. Mean sediment residence time in barchan dunes. Journal of Geophysical Research: Earth Surface, 2014, 119(119): 451-463.
43
Lancaster N. Geomorphology of Desert Dunes. New York: Routledge, 1995, 61.
44
Lancaster N. The dynamics of star dunes: An example from the Gran Desertio, Mexico. Sedimentology, 1989, 36(2): 273-289. DOI:10.1111/sed.1989.36.issue-2
45
屈建军, 凌裕泉, 张伟民等. 金字塔沙丘形成机制的初步观测与研究. 中国沙漠, 1992, 12(4): 20-28.
Qu Jianjun, Lin Yuquan, Zhang Weimin et al. Preliminary observation and study on the formation mechanism of pyramid dune. Journal of Desert Research, 1992, 12(4): 20-28.
46
Zhang W, Qu J, Dong Z et al. The air flow field and dynamic processes of pyramid dunes. Journal of Arid Environment, 2000, 45(4): 357-368. DOI:10.1006/jare.2000.0643
47
Nielson J, Kocurek G. Surface processes, deposits, and development of star dunes: Dumont dune field, California. Geological Society of America Bulletin, 1987, 99(99): 177-186.
48
Zhang D, Narteau C, Rozier O. Morphology and dynamics of star dunes from numerical modeling. Nature Geoscience, 2012, 5(7): 463-467. DOI:10.1038/ngeo1503
49
Fryberger S G. Dune forms and wind regime. In: A Study of Global Sand Seas. Washington: United States Government Printing Office, 1979. 137~170
50
Lancaster N, Kocurek G, Singhvi A et al. Late Pleistocene and Holocene dune activity and wind regimes in the western Sahara Desert of Mauritania. Geology, 2002, 30(11): 991-994. DOI:10.1130/0091-7613(2002)030<0991:LPAHDA>2.0.CO;2
51
Hunter R E, Richmond B M, Alpha T R. Storm-controlled oblique dunes of the Oregon Coast. Geological Society of America Bulletin, 1983, 94(96): 1450-1465.
52
Rubin D, Hunter R E. Bedform alignment in directionally varying flows. Science, 1987, 237(4812): 276-278. DOI:10.1126/science.237.4812.276
53
Fenton L, Michaels T, Beyer R. Inverse maximum gross bedform-normal transport 1: How to determine a dune-constructing wind regime using only imagery. Icarus, 2014, 230(4): 5-14.
54
Fenton L, Michaels T, Chojnacki M et al. Inverse maximum gross bedform-normal transport 2: Application to a dune field in Ganges Chasma, Mars and comparison with HiRISE repeat imagery and MRAMS. Icarus, 2014, 230(4): 47-63.
55
Courrech du Pont S, Narteau C, Gao X. Two modes for dune orientation. Geology, 2014, 42(9): 743-746. DOI:10.1130/G35657.1
56
Gao X, Narteau C, Rozier O et al. Phase diagrams of dune shape and orientation depending on sand availability. Scientific Reports, 2015, 5: 14677. DOI:10.1038/srep14677
57
Lucas A, Narteau C, Rodriguez S et al. Sediment flux from the morphodynamics of elongating linear dunes. Geology, 2015, 43(11): 1027-1030. DOI:10.1130/G37101.1
58
Lv P, Dong Z, Narteau C et al. Morphodynamic mechanisms for the formation of asymmetric barchans: improvement of the Bagnold and Tsoar models. Environmental Earth Sciences, 2016, 75(3): 1-9.
Numerical simulation on sand dune morphodynamic: A case study of ReSCAL model
Zhang Deguo, Yang Xiaoping②,①     
(① Key Laboratory of Cenozoic Geology and Environment, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029;
School of Earth Sciences, Zhejiang University, Hangzhou 310027)

Abstract

In the recent years, the development of sand dune numerical simulation has provided new opportunities for the desert study. Under this background, dune research is no longer limited to the description of a recent depositional state of dunes, but extended to cover investigation on dune morphodynamics. Based on the numerical results of ReSCAL dune model, the following knowledges are obtained:(1) We simulate dunes formed by unidirectional or multi-directional flow regimes, and demonstrate the morphodynamic of superimposed dune emergence and star dune arms extension. Then, we propose that the effect of a complete mechanism between flow and bed form dynamics on dune surface sand flux might be an important reason to shape the superimposed dune and star dune arms. (2) The renewal time of sand particles can be expressed as the ratio of dune length to its migration velocity. However, by tracking individual cells of a 3D cellular automaton dune model, we discover that the depth at which a particle is buried during avalanches is a random variable and the mean residence time is different from the renewal time. In addition, we also calculate that the mean residence time of sediment particles in barchans is equal to the surface of the central longitudinal dune slices divided by the input sand flux. (3) Based on the development of dune orientation's theory and the results of ReSCAL dune model, we elucidate that the changes of sand availability not only shape the dune forms, but also affect the dune growth mechanism and dune arms reorientation in bidirectional wind regimes.
Key words: dune morphodynamic     ReSCAL dune model     superimposed bedforms     star dune     mean residence time     wind regime     sand availability     dune ridge orientation