近年来, 人们已经开发了许多具有不同驱动机制和不同地质时间尺度的数值模拟软件, 以模拟地表的演化过程(Whipple and Tucker, 2002; Tucker and Hancock, 2010; Salles and Hardiman, 2016; Adams et al., 2017; Campforts et al., 2017)。这些模拟软件均可表达为一组数学方程, 将经验数据和概念方法有效结合, 用于重建地貌演化和相关的沉积物通量(Howard et al., 1994; Hobley et al., 2011)。它们的应用范围很广, 除一般的地貌演化研究外, 还可延伸到许多研究领域, 如水文学、土壤侵蚀、山坡稳定性和泥石流、珊瑚礁演变与全球变化、海底滑坡灾害、四维层序地层学等。

但也有许多模拟软件仅专注于侵蚀过程, 特别是有河流系统演化的情况下(Sklar and Dietrich, 2001; Attal et al., 2008; Cowie et al., 2008; Hobley et al., 2011)。这类模拟软件或仅限于河道水系系统的沉积物运移、沉积、固结等过程, 或建立在基于简单扩散方程之上的地表系统的演化过程(Whipple and Tucker, 1999; Willgoose, 2005; Turowski et al., 2007), 导致它们难以将特定地点的观测数据与数值模型联系起来。此外, 涉及到源-汇系统的数值模拟时, 其与相关地层资料的结合也非常少(Howard et al., 1994; Salles et al., 2017), 从而限制了对沉积物深时时期从源到汇过程的全面理解。

Badlands软件的开发旨在解决上述这些缺点, 其重点和核心主要是模型化描述大陆尺度和地质时间(数千年到数百万年)尺度上的地表地貌演化和区域沉积成盆过程(Tucker and Slingerland, 1997; Salles et al., 2011, 2017; Hobley et al., 2017)。Badlands软件被视为一个综合框架, 能提供一个简单、适应性强的数值工具, 以探索地球表面的动态演化过程, 并量化气候、构造作用、侵蚀作用和沉积作用之间的反馈机制。该软件可通过改变边界条件和初始古地形来设置不同演化模型, 并将不同模拟结果与现今或特定地质历史时期的地貌形态、沉积层序和海进海退情况进行对比。

2.2.1 初始模型的建立

Badlands软件使用有限体积法定义续性方程, 应用三角形不规则网格(TIN)的方式来求解地貌方程。初始模型设定时, 需要确定研究区内某一特定地质时间点的古地形并作为模拟演化的初始条件。初始古地形的恢复, 可以结合板块重建、平衡剖面法等进行。之后, 将获得的古地形三维曲面进行网格化处理, 并加载到初始模型中。

2.2.2 可选模块的设置

该软件也可以加入部分可选模块, 包括构造运动、海平面波动、降水分布、波浪条件和挠曲均衡等。

1)构造运动具有可变的三维空间累积位移, 从而可以模拟复杂的空间构造演化, 包括垂直方向(隆起和沉降)和水平方向的运动。同时, 还可以结合其他地球动力学软件或算法(Gerya and Burg, 2007)计算出动力地形, 并将该结果作为深部动力(如地幔流动产生的深部驱动力)的地貌贡献。当施加三维空间位移量时, 该模型使用Thieulot等(2014)提出的节点加密技术。在自然的地表演化之前, 其表面几何形状首先被构造形变所改变。此后, 节点密度随着时间推移而演变, 这可能导致局部产生不均匀的网格分辨率。为了避免这一效应, 需通过添加或删除节点来修改几何表面, 以确保节点的均匀分布。由于将深部动力过程和浅部地形变化相结合, 因此它有助于研究其两者之间的相互关系, 同时也可以更好地对古地貌和古地理演化进行模拟再现。

2)海平面变化一直被认为是影响地球表层系统演化的基本因素之一。在进行长时间的模拟分析时, 由于其长周期的全球海平面波动很可能会受到区域性大陆架和边缘海盆地的扰动(Hallam, 2003; Haq and Al-Qahtani, 2005; Miller et al., 2005; Mü ller et al., 2016), 那么恢复出的区域性海平面变化数据可以作为模型演化的约束条件。为此, 软件中的海平面变化曲线既可以直接应用一些公认的全球海平面曲线数据库, 或由研究人员直接定义。

3)影响地貌变化的因素很多, 而水是极其重要的因素之一。水包括降水和河水, 这两者均可影响相关的沉积过程。该软件不仅可以设置相应的降水量、可蚀性系数、岩石类型, 还可设定通道宽度、泥沙压实率、孔隙度、海洋波浪和潮流对沉积物进行再搬运的相关模拟参数。陆地地貌演化也同时受山体沉积物滑移和沟谷通道运移2种方式影响, 这可以根据研究区的具体情况进行设置。

在大陆上, 降水可控制岩石风化以及河流下游沉积物的分布, 从而驱动河流侵蚀和影响流域内水系网络的形成。该软件中, 降水量的时间变化可以作为恒定值或一组代表空间变化的降雨状况来进行设置。此外, 为了了解降雨量与地形之间的相互作用, 可以选择使用Smith和Barstad(2004)的地形降水线性模型。例如, 降水模式和地形的耦合演化可用于量化山地地貌中气候、侵蚀和构造三者之间的相对重要性。

4)波浪的状态可以影响数千年来海洋沉积物的运移, 这里采用的方法依赖于当前区域内多年波浪状态参数的平均值。波变换模型通常设计在5~50年的时间间隔内进行, 目的是通过施加一系列波浪强迫条件来模拟真实波场。在任何给定的时间间隔内, 定义每个深水波条件的活动百分含量以及波高。然后, 再根据测深法计算相关波参数。上述强迫机制将直接控制沉积物运移、相关的地层结构以及碳酸盐岩的发育。

5)在模型的演化过程中, 可能存在古山脉的剥蚀和沉积物负荷再分配。一般认为, 地表地形的起伏造成的载荷差异将在地壳深部乃至更深的部位得到充分补偿, 过多的地表负荷会导致在补偿界面之上要有等量的质量亏缺才能达到静态平衡, 反之亦然。因此, 需要考虑重力均衡作用对地表产生的影响。该软件也同时考虑到岩石圈表层沉积物负荷再分配对地壳形变的影响, 采用挠曲均衡模块进行计算, 并可根据具体的地质背景, 自由选择挠曲或非挠曲的地壳均衡模块进行模拟计算。

东海陆架盆地位于华南大陆东部的邻海区域, 西接浙闽隆起区, 东邻钓鱼岛隆褶带, 东西宽250~300 km, 整体呈NNE向展布, 为中国东部海域中、新生代叠合含油气盆地(江东辉等, 2017; 刘泽, 2018)。东海陆架盆地南部区域主要受NE— NNE向断裂控制, 区域上表现为NNE向隆拗相间的特征, 自西向东包括多个构造单元: 瓯江凹陷、雁荡低凸起、闽江凹陷、台北低凸起、基隆凹陷等(图 1)。

图 1

Fig.1

	Figure Option ViewDownloadNew Window
	图 1 东海陆架南部盆地构造单元划分(据刘泽, 2018)Fig.1 Structural units of the southern East China Sea Continental Shelf Basin(after Liu, 2018)

东海陆架盆地在中生代曾广泛接受沉积, 厚度较大, 总体具有“ 东厚西薄、南厚北薄” 的特征。沉积中心分别位于东部的基隆凹陷、闽江凹陷一带和西部的瓯江凹陷带内。区域内沉积层一般厚2000~5000 m, 最厚可达6000 m。东部的基隆凹陷、闽江凹陷带的平均沉积厚度为5000 m左右, 而瓯江凹陷平均沉积厚度较小, 约为3000 m。东海陆架盆地南部多发育滨浅海相沉积物, 由南向北逐渐变薄。

垂向上, 研究区则表现为双层盆地结构。雁荡低凸起和闽浙隆起带在中生代期间对研究区的盆地构造格架有重要影响。通过对东海陆架盆地西部地震和重磁资料的综合地球物理解释, 认为雁荡低凸起在中生代晚期发育于瓯江凹陷和闽江凹陷之间, 长约170 km、宽15~50 km, 呈NE向展布。其可能为低山或高地形地貌形态, 分割瓯江凹陷和闽江凹陷(杨传胜, 2014), 形成了2个低凸起和3个凹陷相间分布的盆地构造格局(图 1)。

为了获取中生代早期的古地形, 笔者通过质量平衡法和平衡剖面法, 并结合研究区域内现有的地震剖面、测井、平衡剖面等资料对其进行恢复。具体来说, 中生代研究区的陆地区域由于一直处于剥蚀状态, 故应用质量平衡法对该期间的陆地区域古地形进行恢复。而在海域范围内, 由于存在大量的地震剖面资料, 基于封闭体系中体积守恒、面积守恒和线长守恒3项基本原则(Cristallini and Allmendinger, 2001; N’ stor et al., 2003), 先进行二维平衡剖面的恢复, 再对大量恢复后的平衡剖面海域区域的盆地基底数据进行分段三次Hermite多项式插值, 进而拟合出更为精细的中生代海域古地形。之后, 将获得的海、陆古地形结合起来进行三维曲面网格化处理, 并加载到初始模型中。文中设置的三维模型分辨率为1 km。

由于海平面升降也会显著地影响到大陆边缘盆地沉积演化的整个过程, 所以它是模型设置中非常重要的参数(Watts and Thorne, 1984; Haq et al., 1987; Hallam, 2003; Haq and Al-Qahtani, 2005; Miller et al., 2005)。首先, 笔者使用Didger软件对Haq等(1987)恢复的全球海平面变化曲线进行数字化, 之后对数据采用初始化函数得到标准化曲线(图 2中绿色曲线), 并截取中生代从200— 100 Ma的海岸线变化数据(图 2中右侧曲线部分)存储为csv文件格式以便输入模型时进行加载。

图 2

Fig.2

	Figure Option ViewDownloadNew Window
	图 2 东海陆架盆地南部中生代古地貌模型中输入的全球海平面变化曲线(据Haq et al., 1987)Fig.2 Global eustatic curve input in the Mesozoic paleogeography model of southern East China Sea Continental shelf sea(after Haq et al., 1987)

区域降水量是判断区域气候类型的重要参数, 并控制岩石风化速率以及河流的径流量, 从而驱动河流侵蚀和河流体系的演化(Willgoose, 2005; Whipple, 2009)。文中模型的降水参数设置为每年1000 mm(陈云华, 2008; Whipple, 2009)。地表地貌的塑造不仅和降水量有关, 同时也和地表的风化剥蚀率有着密切的关系。Э 是模型中的侵蚀率系数, 用来衡量冲蚀速率的大小。这个剥蚀系数, 被假定为在整个区域中是均匀的, 其数值为每年0.00005 m(1-2n)。m和n的值分布表示在恒定的泥沙通量和输沙能力条件下的河床底部剪切应力和河床下切率大小之间的关系, 也就是指示了冲蚀速率随河床底部剪切应力呈比例地变化。m和n都为正值, 在文中m和n的值设置为0.5。

另外, 还要获得一个能够较好反映研究区深部动力学过程的区域动力数据。笔者将结合其他地球动力学软件和有限差分算法(Gerya and Burg, 2007)计算出的动力地形模拟结果作为深部动力的来源(Rubey et al., 2017; Liu et al., 2018)。计算过程中的动量守恒和连续性方程解的边界条件是非滑动边界条件。

文中设定研究区域地壳的有效厚度为50 km(蒋玉波等, 2013), 平均沉积物和地幔密度设定分别为2700 kg/m³和3500 kg/m³。在每一个模拟时间步结果保存后, 可以通过加载这个模块, 计算重力均衡作用对地表产生的影响。

从模拟结果来看, 在深部地幔物质运动和岩石圈尺度的构造运动共同作用下, 浙闽隆起发生裂陷, 在其东侧形成了雁荡低凸起。这点在模拟结果剖面和地震剖面上有着很好的一致性, 具体表现为: 雁荡低凸起的西侧为大断距的瓯江凹陷, 其东侧为具凹陷特征的闽江凹陷, 而两者之间的雁荡低凸起为相对高地形, 并将2个凹陷分开。闽浙和粤东沿海地区一直处于山地地形, 几乎没有沉积盆地形成, 而该区域实际地层中也缺乏这个时期的地层记录, 表明该区基本没有这一时期的沉积盆地(蒋玉波等, 2013), 当时地势高, 长期处于剥蚀状态(图 3)。另外, 从集中在闽江凹陷、瓯江凹陷内部的钻井揭示的岩性来看, 白垩系是以陆相碎屑岩为主, 这也说明了瓯江凹陷和闽江凹陷的沉积物可能是来自于雁荡低凸起和闽浙隆起带的陆相地层风化剥蚀的产物(杨传胜, 2014)。

图 3

Fig.3

	Figure Option ViewDownloadNew Window
	图 3 东海陆架盆地南部中生代古地貌模拟结果(据刘泽, 2018) 海拔高程是根据Haq等(1987)的海平面曲线进行计算的Fig.3 Simulated Mesozoic paleogeomorphy in the southern East China Sea Continental Shelf Basin(after Liu, 2018)

随后, 进一步将沉积体系模拟结果与中生代残留沉积地层及地貌演化等进行对比。东海陆架盆地的中生界分布特征为南厚北薄(图 4-b), 其南部中生界整体呈北东向残留分布, 东西宽90~110 km, 南北长约 550 km, 面积约4× 10⁴ km², 沉积层厚度一般在2000~4000 m之间(毛建仁, 1994; 杨传胜等, 2012; 李三忠等, 2013)。从模拟结果来看, 中生代残留沉积与模拟沉积区平面分布具有比较好的可对比性: 两者均为NE向条带状展布, 自西向东依次为闽浙隆起带、瓯江凹陷、雁荡低凸起和闽江斜坡带; 区域性盆地形态也极为相似; 中生界沉积厚度一般为2500~5000 m, 最大厚度为5500 m左右。上述这些特征与中生界残留沉积特征较为相符(图 4-a)。

通过对比分析后发现, 模拟结果与已知的中生代的盆地构造格架特征和沉积分布规律具有较好的一致性, 这也说明了文中的模拟结果可能较好地反映了研究区中生代的地质演化过程。

图 4

Fig.4

	Figure Option ViewDownloadNew Window
	图 4 东海陆架盆地南部中生界残留沉积厚度分布与模拟结果对比(据刘泽, 2018)Fig.4 Comparison of thickness of Mesozoic residual sedimentary thickness with simulation results in the southern East China Sea Continental Shelf Basin(after Liu, 2018)

在200— 175 Ma的侏罗纪早期(图 5-a), 闽浙隆起带逐渐在沿海边缘开始隆起, 但隆升不明显, 在北部沿海区域出现一系列丘陵地带。闽浙隆起区为主要物源区, 在东南部地区主要发育大陆边缘沉积盆地, 总体处于滨海— 浅海陆架沉积环境。盆地形态北高南低, 研究区的西北侧(对应现今中国东南沿海)发育NW向河流, 在河流下游存在山间盆地, 在盆地南部发育由北向南流向的河流, 在河流下游发育沉积盆地。这一模拟结果可能正好对应于浙闽陆域分布较广、西北部以陆相为主、多呈NNE向条带状分布的特征(江东辉等, 2017)。

图 5

Fig.5

	Figure Option ViewDownloadNew Window
	图 5 东海陆架盆地南部中生代地貌演化模拟结果Fig.5 Simulated evolution of Mesozoic geomorphology in the southern East China Sea Continental Shelf Basin

在175— 150 Ma的侏罗纪晚期(图 5-b), 闽浙隆起带北部隆升明显, 仍为主要物源区, 东侧近海沉积范围在之前的基础上有所扩大, 近海河流三角洲向东拓展。在整个侏罗纪, 未见雁荡低凸起在隆升, 仍为一系列丘陵地貌特征。此阶段的雁荡低凸起与浙闽隆起区可能连为一体。前人观点认为, 从早侏罗世早期开始, 陆域西南部海侵面积扩大, 但多数地区仍处于隆起状态。早侏罗世末期海侵结束, 海水退出, 至中侏罗世全境上升, 气候也由早期的温暖潮湿逐渐转为炎热干燥, 仅有陆相沉积, 多为山间盆地沉积(江东辉等, 2017)。这与笔者的模拟结果相对应。东海陆架盆地南部发育滨浅湖相, 东侧的大陆边缘前缘区域有多个河流如海, 形成多个河流三角洲。这很可能是模型东北部分区域地层中多为以紫红色为主的陆相细碎屑岩建造的原因(蒋玉波等, 2013)。

大约在150— 125 Ma(图 5-c, 5-d), 闽浙隆起带和雁荡低凸起开始逐渐隆升, 成为其东侧的闽江凹陷和斜坡带的物源供应区。南部的部分沉积区物质抬升到地表, 遭到风化剥蚀。盆地构造格局也发生了变化, 大陆边缘沉积盆地以西, 古海岸山脉的两侧是山间盆地。瓯江凹陷和闽江斜坡带为主要沉积区域。在此阶段, 研究区的西部主要沿着河流发育冲积平原、三角洲沉积, 东北部沿海地区沉积区扩大显著, 雁荡低凸起东侧沿海区域则以冲积平原— 三角洲发育为主。至此, 瓯江凹陷和闽江凹陷的中生代盆地形态基本形成。

以上应用实例表明, 在强烈的深部构造运动背景下, Badlands模拟软件可能是古地形恢复的一种有效方法。该方法可以获得更为连续的古地理、古地貌信息, 也可为其他盆地的沉积过程、能源矿产勘探的研究提供有益的思路。