风成沙丘数值模拟研究综述

管 超，陶彬彬，刘 丹，周炎广

（北京师范大学资源学院，北京 100875）

关于风成沙丘的研究，人们从静态或当前状态的沙丘场中的沙丘形态以及分布规律认识到风成沙丘的形成是一个长期的过程，短期内无法解决。早期采用野外观测的方法进行定性研究，这有很多弊端。有些学者采用风洞实验以及水槽试验来模拟沙丘的动力环境，但是其问题在于相似性不够，虽然利用特征尺度的概念可以分别得到大气环境和水环境下的特征尺度值［1］，保证了几何相似性，但是由于尺度缩小后动力发生变化，动力相似并没有很好解决；还有，由于野外观测环境和仪器精度问题的限制，特别是沙丘背风坡流场的湍流度较高，野外观测和风洞试验都很难对流场的详细结构和剪应力分布给出满意的结果。因此既能够保证几何相似和动力相似又能够长时间、大范围的实现风成沙丘的仿真研究，计算机数值模拟的引入就显得极其重要，同时也可以很快地给出不同沙丘的形状参数。就目前的研究现状而言，计算机数值模拟主要从以下两个方面来实现：一方面通过模拟沙丘形成各个阶段表面流场的状况来反演出沙丘各个阶段的形态变化；另一方面通过“沙体元”的迁移计算其沉积概率来实现“沙体元”的堆积与侵蚀，进而得到沙丘的形成过程。这两方面分别是从沙丘的外在流场和内部积蚀的角度来模拟，相关理论已经较为完善，但仍有局限性，与野外实测数据仍存在误差。这主要是由于一些人为规定的参数与实际规律有较大差异。有些学者将前两种方法相结合，即内部积蚀和外在流场是相互影响且共同作用于基底面，进而得到沙丘的形成过程。这种方法已经有了相关研究，但是应用范围有限，并没有覆盖到各个沙丘类型，因此，如何通过较为成熟的模拟方法从不同的初始环境因子得到不同的沙丘类型，是下一阶段研究的重点；同时如何跨尺度的从沙粒传输、风沙流运动模拟到沙丘形成再到沙丘之间的作用，最后到沙丘场的形成，这也是今后研究的重点所在。我们在梳理自20世纪中期至今风成沙丘计算机数值模拟研究的基础上，通过分析不同模拟方法的优劣和不同沙丘类型的模拟结果来浅析未来研究方向，希望为我国风成沙丘计算机数值模拟等领域的研究提供借鉴。

1 数值模拟方法

1.1 非线性离散模型

非线性离散模型主要是利用边界层气象学中低矮地表气流离散模型来反演沙丘流场，从而对沙丘形态进行模拟［2］。近地表流场、床面形态、沙尘沉积和传输是复杂的时空多尺度系统，存在自组织的动态平衡关系［3］。模型中主要参数为原始床面和初始条件，原始床面利用网格数值模型建立，初始条件利用气流离散模型和输沙率模型确定，所以非线性离散模型就是通过计算床面上离散网格的沙粒侵蚀状况来反演沙丘形态变化参数、前移速度和表面风沙流的状况［4］。在模型的3个基本要素中，沙丘流场决定着局地沙尘输运强度，进而影响沙丘形态，因此流场模型的建立尤为重要［3］。早期学者多针对新月形沙丘进行数值模拟［5-8］，模型参数中床面通过实测的雏形新月形沙丘的DEM建立，地表气流离散模型采用阿姆斯特技术大学的气流离散模型［9-10］，气流离散模型可以反演沙丘表面气流、流形和剪切压力的分布状况。输沙率模型从最早Bagnlod的模型开始在不断的改进，包括考虑坡度对床面的影响而提出的坡面输沙率模型［5］。Weng等认为，沙丘近地表风速分布并不遵循对数分布，迎风坡坡角处的输沙量并不仅是取决于摩阻速度［11］。Stam认为气流剪切力与沙丘相对高度变化间的延时性是沙丘发育演化的因素［12］。李振山等认为非线性离散模型中气流模式一般用于平缓沙丘地貌，并不适用于丘体背风侧气流状况模拟［13］。近年来，通过流体力学等涡流理论的建立，开展了新月形沙丘表面全流场数值计算的研究，江丽娟等基于商用CFD流体计算软件FLUENT，采用湍流模型定量给出了不同来流风速及沙丘高度对背风侧回流区长度的影响，认为回流区长度与沙丘高度的比值随着高度的增加而增加［3］；周晓斯等采用Smagorinsky亚格子尺度涡黏模型的大涡模拟（LES）方法，对简化的实验缩比小尺度沙丘模型绕流气相流场进行了数值研究，进一步验证了背风侧回流特征［14］，然而其模拟结果与真实尺度下野外观测结果的验证精度仍有待提高。

图1 Werner理论模型

2.2 网格动力模型

元胞自动机（Cellular automotion，CA）是网格动力模型的代表，它具有强大的空间运算能力，常用于自组织系统演变过程的研究。它是一种时间、空间、状态都离散，空间相互作用和时间因果关系都为局部的网格动力学模型，具有模拟复杂系统时空演化过程的能力［15］。Hagerstrand在空间扩散模型的研究中首次采用了类似于元胞自动机的思想［16］。20世纪70年代，Tobler将元胞自动机的概念引入地理学研究［17］，他用CA模型来模拟当时美国五大湖边底特律地区城市的迅速扩展，他认为类似元胞自动机的地理模型的采用是分析模拟地理动态现象的一次方法革命。近年来，国内冯永玖、黎夏和孙战利等人在元胞自动机与地理时空演化特性的交叉研究方向上都做了深入的研究［18-20］，探讨了地理元胞自动机在空间复杂现象中的知识发现、数据挖掘和时空过程建模等问题［21-23］，更加完善了CA应用在地学方面的理论。

在沙丘地貌模拟方面，沙波纹、流动沙丘和沙丘场已应用CA模型来模拟［24-26］。最早提出较为完善的沙丘CA模型的是Werner［27］，其理论主要涉及到沉积概率和休止角判断的问题。Bagnold曾提出由于地面的性质不同，如果没有沙物质（石质地面）则跃移中沙粒的回弹程度较大，其沉积概率较低，如果有沙物质则沉积概率较高［28］。在Dune-CA中，沙丘是由“沙体元”形成的，它们被限制在一个无侵蚀能力的床面网格内（如图1），一个沙体元是在所有的沙体元中被随机选择用于传输的，沙体元移动的网格数记为，移动之后沉积的位置取决于沉积概率，它和该位置的沙体元数量有关；如果移动后该位置没有沙体元，那么其沉积概率值一定小于该位置有一个沙体元的沉积概率值；如果该沙体元没有沉积，那么它将再次重复移动距离网格直到沉积发生，之后下一个随机选取的沙体元才开始移动。当沙体元移动到阴影区的时候，需要判断该位置的坡度是否大于休止角，从而判断该沙体元是否会滑落。

从沙丘形态上来看，Werner得到的模拟结果与自然界实际的沙丘形态十分相似，即在单一风向作用下，沙源供给量较多则会形成横向沙丘，较少则会形成新月形沙丘；在多风向作用下，2个风向作用会形成线形沙丘，3个风向作用会形成星形沙丘。近年来，Thomas和Chris基于Werner沙丘模型，采用C++和Python开发出的GUI界面沙丘模拟工具，并集成在ArcGIS中，提出崩塌算法、边界周期循环算法等内容［29］，使得模拟结果更接近真实（如图2所示）。

在Werner的基础上，Anderson提出，在沙丘场中，随着沙源的提供、风向的改变，所有沙丘通过沙丘间的融合、分离和链接等最后均可形成横向沙垄和纵向沙垄，即沙丘场中存在这两类吸引子，解决了风速的变化决定了最终吸引子走向的问题［30］。然而Anderson并未解决Werner对于单一风速的限制，在模型中由于风速决定了“沙体元”的传输距离，即决定了每次迭代后床面沙体元的重新分布，自然也决定了下一次迭代的不同位置的沉积概率。显然，风速设为定值是与现实不符的，同时由于Werner和Anderson的CA模型均默认模拟中砂粒粒径为定值，导致侵蚀沙量也为常数，这也是与现实不符的。由于植被在沙丘形态变化过程中起着重要作用，而植被与沙丘的相互作用机理和反馈关系尚不清楚［31］，故限制其发展。Werner在CA模型的基础上加入植被因子，得到了加入植被因子CA模型的雏形［27］，随后Baas建立CA扩展模型，以加入植被因子的新月形沙丘和抛物线形沙丘［33］为例，模拟了沙丘的侵蚀堆积和演变过程［32-33］。图3为新月形沙丘结合植被因子的CA扩展模型结果［32］。近年来，Mark A.Fonstad提出“CAs”的概念来丰富CA模型，将生态因子加入CA模型中来模拟地貌—生态耦合系统，并以风沙地貌的CAs模型为例，模拟出更符合真实生态环境下的风成沙丘地貌和沙丘场［34］。Niedl和 Baas利用 DECAL（discrete ecogeomophic Aeolian landscape model）模型，模拟了新墨西哥州抛物线形沙丘的形态演变过程，认为较高植被增长速度利于沙丘的形成，较低植被增长速度将导致横向沙丘的形成［24］。Baas和Nield又将面向对象方法引入到DECAL模型中，对模型的应用范围、环境参数输入和状态转换规则等进行了改进，使其模拟预测的结果与现实场景更加接近［35］。

图2 Werner模型的计算机模拟结果

2.3 系统动力学模型

系统动力学模型是通过分析系统结构选取适当因素，建立它们之间的关系，并在此基础上建立微分方程，进而模拟和预测不同参数输入下沙丘演变过程和趋势［36］。该方法多使用于封闭系统单体沙丘的模拟和计算，可将野外观测的植被盖度和侵蚀堆积状况等数据作为参数输入。Herrmann和Sauermann建立了沙丘演变偏微分方程［37］，并对结果进行验证。Duran等将野外观测的沙丘植被盖度和侵蚀堆积数据输入砂质搬运模型中，模拟巴西东北部海岸抛物线形沙丘的发育过程，并分析新月形沙丘向抛物线形沙丘演变模式及风蚀坑坑后集沙区发育成抛物线形沙丘的过程［38］。

图3 新月形沙丘CA扩展模型结果

3 结束语

就3种模型而言，非线性离散模型主要应用于流场的模拟，在新月形沙丘上应用较多，新月形沙丘中最重要的就是背风坡分离形成水平涡流，因此早期学者利用背风坡的无分离湍流封闭模型数值求解N-S方程，得到其表面剪应力分布以及风场扰动特征，但对于实际情况中的背风坡强湍流，这种数学模型是不适用的，也就不能得出沙丘背风坡的剪应力分布。江丽娟等近年来引入CFD流体计算软件FLUENT后得到了一些结论，指出了背风坡回流区与沙丘高度的关系［3］，但与实际沙丘状况相比模拟精度仍有待验证。非线性离散模型只适用于局部气流模拟，对于大范围环境影响而言模拟精确性不高。

网格动力模型主要应用于环境变量对沙丘的作用，主要体现在植被的固定对抛物线沙丘的演变过程的研究，得出的结论主要涉及沙丘的蚀积和演变过程［33，39］，但由于基础数据为栅格数据，模型精度常常不高，对于高精度的流畅研究不适用。

非线性模型在金字塔沙丘的流场研究应用较多，但一般需要结合其他模型分析。如柳本立等利用剪切压力传输模型（SST）证明了金字塔沙丘是外来风力和自身二次流共同作用的产物［40］；张德国等提出一种由元胞自动机模型（用来模拟颗粒运动）和栅格气体模型（用来模拟流体运动）组合而成的混合模型，通过对金字塔形沙丘三维形态分析，得出沙丘臂膀的形态和延伸速度依赖于主风向的变换周期的结论等［41］。

尽管目前风成沙丘计算机模拟已经取得一定的进展，但是仍然存在许多问题。如图4所示，

图4 风成沙丘数值模拟存在的问题

高精度参数目前多为基于野外样方的调查经验参数，尚需要结合3S技术，从宏观、全区的角度提取信息；而对于模型验证标准，多数学者模拟出的结果仍采用与风洞试验数据相比较的方法进行验证，由于尺度不同，难以保证精度，建议通过实测沙丘形态数据建立沙丘Kappa系数评价体系，实测数据可以选择某一典型沙丘的连续时序数据，通过高精度RTK-GPS获得空间离散点坐标，经过插值得到多期连续DEM数据，最后对比模拟网格数据和真实DEM数据得到Kappa系数，通过评价体系来得到验证精度。单一的模拟沙丘的形成过程而不涉及小尺度的沙粒运动、风沙流运动以及大尺度的沙丘场形成是不能够完全模拟出精度较高的沙丘，因此如何跨尺度的研究仍是进一步的重点。最后复合型沙丘动态变化及沙丘之间相互转化目前也没有实现，如抛物线形沙丘与新月形沙丘和横向沙丘之间的相互转化［42-46］，如何仅仅通过环境变量达到这一目的，仍然需要在沙丘形成机理与环境变量的相互关系上进一步研究。

［1］ SAUERMANN G，KROY K，HERRMANN H J.Continuum saltation model for sand dunes［J］.Physical Review E，2001， 64（3）：301-305.

［2］ ZEMAN O，JENSEN N O.Modification of turbulence characteristics in flow over hills［J］.Quarterly Journal of the Royal Meteorological Society， 1987， 113（475）：55-80.

［3］ 江丽娟，马高生，郑晓静.新月形沙丘表面流场的数值模拟［J］. 兰州大学学报，2010，46（1）：48-52.

［4］ 李振山，董治宝，陈广庭.风沙运动数值模拟研究的进展［J］. 干旱区研究，1997，14（1）：63-68.

［5］ HOWARD A D，MORTON J B.Simulation model of erosion and deposition on a barchan dune［J］.NASA CR-2838，1977：1-73.

［6］ HOWARD A D，MORTON JB，GAD-EL-HAK M，et al.Sand transport model of barchan dune equilibrium［J］.Sedimentology，1978，25（3）：307-338.

［7］ HOWARD A D，WALMSLEY J L.Simulation model of isolated dune sculpture by wind［C］//Proceedings of international workshop on the physics of blown sand，University of Aarhus，Aarhus，Denmark.1985：377-391.

［8］ FISHER P F，GALDIES P.A computer model for barchan-dune movement［J］.Computers&Geosciences，1988，14（2）：229-253.

［9］ WALMSLEY J L， HOWARD A D.Application of a boundary layer model to flow over an eolian dune［J］.Journal of Geophysical Research：Atmospheres（1984—2012）， 1985，90（D6）：10631-10640.

［10］ WIPPERMANN F K，GROSS G.The wind-induced shaping and migration of an isolated dune：a numerical experiment［J］.Boundary-Layer Meteorology，1986，36（4）：319-334.

［11］ WENG W S，HUNT JC R，Carruthers D J，et al.Air flow and sand transport over sand-dunes［M］.Aeolian Grain Transport∶Springer Vienna，1991：1-22.

［12］ STAM JM T.On the modelling of two-dimensional aeolian dunes［J］.Sedimentology，1997，44（1）：127-141.

［13］ 李振山，倪晋仁.国外沙丘研究综述［J］.泥沙研究，2000（5）：73-80.

［14］ 周晓斯，王 元，李志强.小尺度新月形沙丘背风侧流场特性的大涡模拟分析［J］.西安交通大学学报，2013，47（7）：114-119.

［15］ 周成虎，孙战利，谢一春.地理元胞自动机研究［M］.北京：科学出版社，2001.

［16］ HAGERSTRAND T.A monte-carlo approach to diffusion［J］.Archives Europeennes de Sociologie，1965，6（1）： 43-67.

［17］ TOBLER W R.A computer movie simulating urban growth in the Detroit region［J］.Economic Geography，1970， 46：234-240.

［18］ 冯永玖，童小华，刘妙龙，等.基于GIS的地理元胞自动机模拟框架及其应用［J］.地理与地理信息科学，2010，26（1）：41-43.

［19］ 黎 夏，叶嘉安.知识发现及地理元胞自动机［J］.中国科学（D辑），2004，34（9）： 865-872.

［20］ 孙战利.空间复杂性与地理元胞自动机模拟研究［J］. 地球信息科学，1999，1（2）：32-37.

［21］ 黎 夏，叶嘉安，刘小平，等.地理模拟系统：元胞自动机与多智能体［M］.北京：科学出版社，2007.

［22］ 全 泉，田光进，沙默泉.基于多智能体与元胞自动机的上海城市扩展动态模拟［J］.生态学报，2011，31（10）：2 875-2 887.

［23］ 杨青生，黎 夏.基于动态约束的元胞自动机与复杂城市系统的模拟［J］.地理与地理信息科学，2006，22（5）：10-15.

［24］ NIELD J M，BAASA C W.The influence of different environmental and climatic conditions on vegetated aeolian dune landscape development and response［J］.Global and Planetary Change，2008，64（1）：76-92.

［25］ BISHOPSR，MOMIJI H，CARRETERO-GONZÁLEZ R，et al.Modelling desert dune fields based on discrete dynamics［J］.Discrete Dynamics in Nature and Society，2002，7（1）：7-17.

［26］ 袁 鸿，王志伟.工程力学科学与实践［M］.广州：暨南大学出版社，2010：21-39.

［27］ WERNERB T.Eolian dunes：Computer simulations and attractor interpretation［J］.Geology， 1995， 23（12）：1 107-1 110.

［28］ BAGNOLDR A.The physics of blown sand and desert dune［M］.London：Methuen，1941.

［29］ BARCHYN T E，HUGENHOLTZ C H.A new tool for modeling dune field evolution based on an accessible，GUI version of the Werner dune model［J］.Geomorphology，2012，138（1）： 415-419.

［30］ ANDERSON R S.The attraction of sand dunes［J］.Nature，1996，379：24-25.

［31］ HESP P.Foredunes and blowouts：initiation，geomorphology and dynamics［J］.Geomorphology，2002，48（1）： 245-268.

［32］ BAASA CW.Chaos，fractals and self-organization in coastal geomorphology：simulating dune landscapes in vegetated environments［J］.Geomorphology，2002，48（1）：309-328.

［33］ BAASA CW.Complex systems in aeolian geomorphol－ogy［J］.Geomorphology，2007，91（3）：311-331.

［34］FONSTAD M A.Cellular automata as analysis and synthesis engines at the geomorphology-ecology interface［J］.Geomorphology，2006，77（3）：217-234.

［35］ BAASA C W，NIELD JM.Ecogeomorphic state variables and phase space construction for quantifying the evolution of vegetated aeolian landscapes［J］.Earth Surface Processes and Landforms， 2010， 35（6）：717-731.

［36］ FORRESTERJ W.Urban Dynamics［M］.Cambridge，Massachusetts，USA：The MITPress，1969：10-33.

［37］ HERRMANNH J，SAUERMANNG.The shape of dunes［J］.Physica A：Statistical Mechanics and its Applications，2000，283（1）：24-30.

［38］ DURÁN O，SILVA M V N，BEZERRA L J C，et al.Measurements and numerical simulations of the degree of activity and vegetation cover on parabolic dunes in north-eastern Brazil［J］.Geomorphology， 2008， 102（3）：460-471.

［39］ DURÁN O，HERRMANN H J.Vegetation against dune mobility［J］.Physical Review Letters，2006，97（18）：188 001.

［40］ 柳本立，张伟民，彭 飞，等.金字塔沙丘流场的三维数值模拟.中国沙漠［J］.2011，31（2）：386-392.

［41］ ZHANG D，NARTEAU C，ROZIER O，et al.Morphology and dynamics of star dunes from numerical modelling［J］.Nature Geoscience，2012，5（7）：463-467.

［42］ HACK J T.Dunes of the western Navajo country［J］.Geographical Review，1941，31（2）：240-263.

［43］ ANTON D，VINCENT P.Parabolic dunes of the Jafurah desert，Eastern Province，Saudi Arabia［J］.Journal of Arid Environments，1986，11：187-198.

［44］ GAYLORD D R，STETLER L D.Aeolian-climatic thresholds and sand dunes at the Hanford Site，southcentral Washington，USA［J］.Journal of Arid Environments，1994，28（2）：95-116.

［45］ MUCKERSIE C，SHEPHERD M J.Dune phases as time-transgressivephenomena，Manawatu，New Zealand［J］.Quaternary International，1995，26：61-67.

［46］ ANTHONSEN K L，CLEMMENSEN L B，JENSEN J H.Evolution of a dune from crescentic to parabolic form in response to short-term climatic changes：Råbjerg Mile，Skagen Odde，Denmark［J］.Geomorphology，1996，17（1）：63-77.