0 引言

中国边缘海大陆架在晚更新世时期曾是海岸平原，在古长江、古黄河泥沙填充下形成了陆架堆积体，并在全新世发育了南黄海辐射沙脊群、废黄河三角洲和长江三角洲^[1-2]。关于古长江三角洲的核心—南黄海辐射沙脊群，早在1957年就有研究发现长江口外水下地形呈扇形分布的古代水下三角洲，前缘水深-50 m，平面中心在32.3º N，面积为7 000 km^{2 [3]}。此后五十多年以来，不断有新的发现和认识。到21世纪，根据区域的地质钻孔，揭示了南黄海辐射沙脊群水深0~25 m，以琼港（32°40′ N，120°54′ E)为中心，褶扇状向海展布，面积22 470 km²。并认为黄洋沙潮流通道是古长江河谷，沉积体形成的时代成早于43 ka B.P. ^{[2, 4]}。然而，根据点状分布的地质钻孔分析和重建，对南黄海—江苏海岸的沉积体系和演变的分布、规模以及机制尚不明了。陆上海侵层难以保存，重建的海面低于现代海面15~50 m，它的海岸线在哪里？大河三角洲沉积位置受沉积基准面控制，由于时间上钻孔地层不连续、空间上海侵层不连续，特别是沉积和侵蚀基准面不清楚，末次冰期古长江的河口沉积在哪里？此外，大河三角洲沉积通量受到陆源泥沙控制，而泥沙受到流域侵蚀环境和搬运水动力制约，在冰期和间冰期中，气候降水变化悬殊，影响着入海泥沙量。晚更新世以来的三角洲建造和泥沙沉积与现代有何差异？

本研究拟进行论证和机制探索。由于沉积体系是多种过程的产物，动力过程与沉积记录的关系十分复杂。采用物理机制进行模拟，已成为过去20多年来包括STRATAFORM在内的多个国际合作项目研究的重要内容^[5-6]。在对现代海洋、海岸沉积模型的基础上，一些学者应用到不同地质时间尺度的边缘海、海岸河口以及三角洲沉积进行数值模拟。例如对晚第三纪^[7]、末次冰盛期^[8]、末次冰消期^[9]以及过去1000年^[10]等不同时期的模拟尝试，获得了对这些重大气候期沉积体系演变的机理认识。由此，本文在前人基于气候—海面—沉积系统的概念模型、动力模型研究基础上，根据气候水文学、沉积学等原理以及泥沙沉积的动态高程计算，构建气候冰川驱动—东黄海地海系统响应—河流沉积建造的数值模式，模拟冰川气候、构造沉降、沉积压实等复杂效应下的海面特征、陆源泥沙沉积和海岸线分布，从动力学机制上揭示南黄海海岸沉积演变。

1 模式架构

首先基于能量守恒定律和质量守恒原理，构建冰川型海面变化模式，模拟晚更新世以来不同气候期南黄海—江苏海岸岸线变化和空间格局，分析制约海面与气候相关点与突变点，以度量和刻画太平洋边缘海海面变化过程和特征。其次，基于水沙沉积动力学和地貌高程学原理，构建沉积物表层高程和物质输移模式，计算潮流海岸和三角洲海岸沉积通量变化，对废黄河三角洲、江苏潮滩和长江三角洲不同类型海岸的沉积过程进行模拟。最后，在地质资料重建的海岸线设置的边界场上，应用于晚更新世以来南黄海—江苏海岸地质时期海面变化和沉积分布模拟。

该系统从全球海面变化模式（模式1）到江苏海岸—长江三角洲沉积模式（模式2）分别进行运算，叠置海岸边界条件后模拟海岸和三角洲的变化分布和过程。多地质钻孔的地层层序提供了沉积速率以及分布特征，被用来验证模拟晚更新世地质海岸边界下的江苏海岸—长江三角洲沉积通量变化。模式流程图见图 1。

模式1的第一模块尽管是以全球尺度进行从辐射到冰量和水量模拟，但以点模式（0维模式）模拟。其结果驱动具有三维空间的第二模块运行，在地形场（东黄海沉积后期海面效应模块）上模拟海面变化，其时间尺度采用了0.25 ka步长模拟。模式2对沉积物输移和累积的模拟设置在东黄海地形场上，采用了0.25 ka时间步长。因此模式1和模式2在空间和时间尺度上能够耦合。

2 模式功能设计 2.1 冰川型海面—后期构造沉积效应模式

该模式基于物理学、气候学和沉积学的基本原理，以辐射—气候驱动热量与冰量热力平衡、全球冰川型海面水量平衡、东黄海海面响应等不同层次，构建三个动力模块。

（1） 辐射—冰量热量平衡模块：根据全球太阳辐射温度的能量平衡的Planck辐射定律，计算温度变化：

其中T是理论计算地表温度，α是地表反射率，δ是Stefan-Boltzmann常数（5.67×10^-8 W/m² K⁴），A是表面积（计算全球理论温度采用5.1×10¹⁴ m²），S是太阳辐射总量（1.74×10¹⁷ W）。对晚更新世不同气候期（冰期和间冰期），根据陆面和海面面积变化，海陆采用不同的反射率计算。

其次，根据气候状况（温度变化）和冰量（第四纪冰盖变化），采用冰融化能量平衡模式(能量=比热×质量×温度变化)，计算温度升高时吸收的热量：

其中，c是热容量或比热，水体比热 C_w=4 100J/(kg·℃），冰比热 C_i=4 100J/(kg·℃)。m是冰质量(kg)，ΔT是冰体融化所吸收的温度：ΔT=T_i-T_m，其中T_i冰温度，T_m冰融化温度(0℃)。根据温度升高/降低时冰体吸收/释放的热量变化，针对晚更新世不同地质时期期相对现代的变化，模拟从状态1(地质时期) 到状态2(现代)在热量恒定下冰体积的变化：

（2）冰量—海水物质平衡模块

从冰量变化模拟海洋水体体积变化，拟采用全球固态水与液态水的水量物质平衡：V=V_i+V_w，其中V是全球水体总量，V_i为固态水（冰量），V_w为液态水。地质时期可采用变化量计算，水量变化（%）和冰量变化（%）的总和为100%（V=1）。由于冰密度ρ_i与水密度差异（冰融化水后体积减小1/11），固态水转化液态水的增量ρ_i V_i。获得高程动态变化的微分方程：

其中 h海面变化量(m)，V_i受温度控制的冰量（%）。ρ_i和ρ_w分别是冰密度与水密度。A(h)为模拟的海域面积，它是h的状态函数。由此对微分方程4求数值解，获得晚更新世水量变化以及不同冰量体积变化下海水体积变化。

（3）东黄海沉积后期海面效应模块

依据海面与地面系统变化地质学原理，分别对海侵层后期的海、陆升降效应进行模拟。地面系统(V_L)与海面(V_S)之间的相对海面变化（H_R），定义为古海面在沉积层中的位置^[11]。它以高程计量，是冰川型海面(H_S)与后期海、陆效应(H_L、H_M)的平衡：H_R=H_S- (H_L + H_M)。

H_S 是冰量引起的海水体系变化，通过冰期—间冰期辐射与冰量能量转换—冰量与水量的物质平衡—地形响应多层模式的模拟获得。晚更新世不同气候期（冰期和间冰期）的冰量模拟结果与现代相比为减少。

H_L是陆面系统变化，主要考虑构造沉降和沉积压实两个效应：H_L=H_G+H_P。其中H_G 是构造沉降，其效应使海侵层降低，相对海面增高；以现代百年级的沉降速率估计。H_P为上覆沉积压实作用，其效应使海侵层降低，相对海面增高；以钻孔岩芯层的孔隙率计算。

H_M是海面系统变化，主要考虑海盆体积和陆源泥沙充填二个效应：H_M=H_F+H_B。其中H_F为大陆河流泥沙充填边缘海，其效应使海洋容积减少、海面相对增高。H_B是海盆构造沉降，使海盆体积增大，引起海面相对降低。

综合上述5个分量和变化方向，可计算古海面在沉积层中的位置：

2.2 晚更新世江苏海岸—长江三角洲沉积模式

当沉积速率大于可容空间增速下，形成加积层序，依次向盆地方向进积，形成高位体系域和低位前积楔状体的沉积特征。理想三角洲沉积是指由陆向海的加积沉积(包括水下沉积和水上沉积)，一般是从三角洲和潮滩海岸顶点，沉积物在河口水流、泥沙流和重力作用下向前推移和累积^[12-13]。根据泥沙运动公式^[14]，其垂向沉降的微分方程有：

其中Q_sed是输沙量，_∇是泥沙扩散散度，在垂直方向受到沉积速率（r）与水深（H）控制，沉积速率与泥沙重力惯性（ω）、泥沙质地半径（D）与水流阻力（ν）有关。设α₁是与泥沙质地属性有关的经验系数：。由此6-1式有：

在河流泥沙向海洋输运的加积作用下，泥沙沉积的分布过程可采用逻辑斯谛方程（Logistic）表达^[8]：

其中K是Logistic方程中半饱和系数，它在沉积过程中是与河口形态高差有关的参数，可采用参数率定。输沙量受到流量和含沙量控制，而流量受到水动力变化控制。在地质时期，主要考虑受到流域降水变化的控制。因此，流量变率(Q)采用降水变率(P)相关：Q=α₂ P；泥沙变率采用流量变率相关Q_sed=α₃ Q，进一步用降水变率相关：Q_sed =α₃ Q = α₃ α₂ P，其中α₂与流量、α₃与含沙量相关的经验系数。上式成为：

上式中降水变率（P）采用为状态变量，重力惯性（ω）、泥沙质地半径（D）与水流阻力（ν）根据河口海岸物理属性设置。在6-4式中有3个参数α₁ 、α₂α₃、K，分别采用参数率定。根据常微分方程Logistic式的积分式为：

对6-4式的积分式而言，。因此，可运用MATLAB程式中方程求解函数（Solve）和非线性最小二乘法拟合函数（nlinfit），对6-5式的3个参数进行拟合。

3 边界场资料处理和设置 3.1 地形和海陆岸线

地形边界条件根据中国东部大陆—东黄海等高/等深地形图，GIS划分格点（5×5经纬格点精度的数值化）以计算逐米等深线面积，以多项式拟合变量为深度(m) 和面积(km²)函数，用于模式的状态变量。在地形海拔高度-20~+10 m范围内，以占东中国海总面积542 000 km² 的46.8%面积^[15]设置微分范围和步长。通过各个格点微分方程数值求解，获得不同海拔高度对应的水域面积(图 2)。

全新世中期以来海岸线作为模拟的动态边界条件，采用据江苏沿海和南黄海记录地质钻孔和海岸线分布的重建资料^[1-4]（图 3）。根据6 ka B.P.、4.5 ka B.P.、2 ka B.P.、1.2 ka B.P.以及历史时期海岸线和河口不同时期的位置，进行2D的同化处理，物理量纲为岸线的高程[m]。

3.2 边界场和初始场设置

模式1的边界场，采用了矢量形式的变量输入，包括相对于现代的太阳辐射变率、陆地冰量变率以及研究区域的降水变率。根据文献综述^[16-17]和本文分析，设置了晚更新世（约140 ka B.P.）以来的主要气候期的全球辐射、相对现代的大陆冰量变化、研究区降水变化等边界阈值。表 1列出140 ka B.P.以来主要气候期边界场阈值和相应的参考文献。模式2边界场主要是各个时期的海岸线分布，采用了地质和沉积资料重建的全新世以来江苏海岸线分布，以岸线的经度、纬度相对0 ka B.P.的位置表示。各个时期的岸线分布位置见图 3。

模拟试验的初始场和参数采用设置和率定获得。模式1采用了平均区域构造垂直升降速率、沉积物压实率作为沉积后期高程变化的主要参数，系根据前人研究的成果设置^[30-33]。模式2采用的海陆边界高程差（H₀）、与海岸和河口沉积高差、坡降有关的半饱和系数（K），以及计算沉积通量的入海河流径流量（Q）、入海河流输沙量(Q_sed)。由于区域地貌差异，划分研究区为三段不同类型海岸：北部废黄河三角洲、中部海岸潮滩和南部长江三角洲，因此初始场分三段海岸设置。对6-4式中的参数，系根据多年的江苏长江河口段、废黄河口段水沙资料拟合获得。采用6-5式，当t=1,2,…n年和年降水P=P₁,P₂ …P_n时，设置了河口初始高程（H₀）后，通过非线性最小二乘法拟合函数拟合获得3个参数解。由于江苏海岸中段水沙资料有限^[1]，仅采用类比性设置。表 2列出设置的各项参数和相应的参考文献。

3.3 沉积物埋深

多地质钻孔的地层层序^{[1-4, 10, 33-35]}提供了沉积物埋深的实际情况，被用来验证模拟晚更新世地质海岸边界下的江苏海岸—长江三角洲沉积通量变化。根据钻孔沉积物测定的年龄，划分晚更新世以来的不同气候特征期（35±3 ka B.P.、21±3 ka B.P.、15±1 ka B.P.、12±1 ka B.P.、10±1 ka B.P.、8±1 ka B.P.、6±1 ka B.P.、4±0.7 ka B.P.、2±0.5 ka B.P.、1±0.3 ka B.P.）为标志性年代。采用钻孔样品年代和层序计算沉积速率，利用Kringe几何数学方法做空间内插，确定区域空间不同时代的沉积层厚度和埋藏（图 4）。这些钻孔沉积物显示了自4万年以来的分布特征，特别是在几个变化显著阶段。在35 ka B.P.时长江三角洲北部是一个大海湾，在20 ka B.P.时沉积重心自西向东迁移，到了10 ka B.P.时琼港外形成面状沉积体，在6 ka B.P.时候海岸沉积带向西推进。

4 模拟和验证

地质时期的沉积物分布特征和过程，采用上述二个模式和一个动态边界场，时间总长度140 ka，步长为0.25 ka。采用模式1模拟冰期—间冰期海面，计算冰川型海面，并对后期构造沉降、沉积物压实等复杂效应下的海面高度的时间序列。输出物理量纲为海面垂直高度[m]。同时，采用模式2中模拟河流三角洲泥沙，系根据晚更新世长江和黄河流域气候经历了冰期和间冰期变化，导致降雨干湿并引起入海流量和泥沙巨大变化，驱动陆源泥沙入海和海岸与河口沉积通量变化。输出物理量纲是沉积物埋深[m]。

根据黄海、东海的大陆与海底地形，沉积和构造分量，模拟了140 ka B.P. 以来区域冰川—沉积海面垂直变化，海面变化2D再现。模拟表明，当海面在东黄海下降到-55 m时，渤海海峡关闭。台湾海峡在海面下降到-80 m时关闭。当海面达到-90 m时，黄海成陆。在末次盛冰期的低海面(-145 m)阶段，对马海峡关闭，整个东亚大陆的大陆架暴露成陆，海岸线向东推进了1 000 km (图 5a,b)。这个模拟与全球海面研究成果^{[20, 36]}对比，冰期与间冰期海面变化一致，而区域性的不同变化部分占10%~18%。

根据降雨—流量—泥沙运动特征和参数，模拟不同地质时期大陆入海泥沙的沉积通量、沉积层高程以及区域空间分布。由于沉积物的堆积，岸线和河口向海推进，改变了地形深度和分布。把黄河三角洲海岸、江苏潮汐海岸和长江三角洲海岸不同地貌类型的沉积侵蚀与堆积模拟结果，与海面变化叠置模拟了距今140 ka B.P. 以来江苏海岸线和长江三角洲地貌与沉积长期变化过程。2D模拟输出显示，大陆泥沙输移和堆积造就了海岸带和三角洲的依次变化。在末次冰盛期ca. 19 ka B.P.时海面达到最低，与琉球群岛链呈现内陆湖泊。海面在晚冰期ca. 14 ka B.P.时，从低海面不断上升，形成北黄海内陆湖泊。从全新世早期到中期ca. 6 ka B.P.，海侵到江苏低平原地区，苏北达到范公堤以西(图 5c，d）。

根据江苏沿海径向8个钻孔点的35±3 ka B.P.海侵层记录，对比模拟35 ka B.P.沉积物分布 (图 6)。对比表明，模拟的东海和黄海海面高度在-30.2~-23.7 m。沉积物经3.5万年以来海陆构造沉降、上覆地层压实、河流泥沙充填后，分布高度在-33.5~-22.8 m，不同地点模拟误差在±2.5~4.5 m。对比地质记录，模拟误差小于10%，结果可以接受。从区域上来看，模拟结果能够与钻孔层位（图 4）对比，模拟的空间分布与地质记录基本一致。

5 结语

本文介绍地质时期海面和沉积分布模式设计与实现，分别从气候—泥沙水动力驱动、下垫面海面—地形两个方面构建了模式系统的主要功能，并把两个动力系统模式进行了耦合。模式1的辐射—冰量—水量的0维模式叠加在东黄海地形场进行3D模拟，输出海面高度变化作为模式2边界条件，继续在3D区域地形场上模拟沉积物分布，较好地解决了在两个模式在空间和时间尺度的耦合。

本研究模拟了距今14万年、4万年和1万年不同时间尺度江苏海岸线和长江三角洲地貌与沉积长期变化过程。模拟发现，当海面在东黄海下降到-55 m时，渤海海峡关闭。台湾海峡在海面下降到-80 m时关闭。当海面达到-90 m时，黄海成陆。在末次盛冰期的低海面(-145 m)阶段，对马海峡关闭，整个东亚大陆的大陆架暴露成陆，海岸线向东推进了1 000 km。这个模拟与全球海面研究成果对比，冰期与间冰期海面变化一致，而区域性的不同变化部分占10%~18%。该模拟与数个地质钻孔资料进行的对比验证，表明在晚更新世以来冰期—间冰期—晚冰期—冰后期的海岸分布和沉积记录的变化趋势基本一致。

本研究尝试的物理机制模式对地质时间尺度的模拟，结果表明在晚更新世以来气候环境下，能够模拟由地质推断的沉积体系，表明在河海交互下复合动力能够建造大扬子沉积系统。在今后的工作中，笔者将进一步对模式进行完善和扩展，使其成为一个诊断和预测海面变化和我国海岸响应研究的有效工具。