台湾造山带是世界上最为年轻的造山带之一，同时又是正在活动的弧—陆碰撞性大陆边缘。其位于欧亚板块和菲律宾海板块的交界处，欧亚被动大陆边缘与最前端的吕宋火山岛弧在此处发生碰撞。在其西南侧，南海海洋岩石圈沿着马尼拉海沟向西俯冲于菲律宾海板块之下^{[13, 14]}。其东北面一侧为近东西走向的琉球海沟—弧—盆体系，菲律宾海板块向北俯冲于欧亚板块之下(图 1A)。

图  1  (A) 台湾海域构造图(据文献[11]改编); (B)海岸山脉地质图(据文献[12]改编); (C) AA’构造剖面图(据文献[9]改编，α：连续单元；β：破坏单元；γ：分割单元；δ：混杂岩)

Figure 1.  (A)Tectonic map of Taiwan (modified from reference [11]); (B) The geology of Coastal Mountains (modified from reference [12]); (C)AA' structural section (modified from reference [9], α: coherent units, β: broken units, γ: dismembered units, δ: melanges)

由于斜向弧陆碰撞，台湾的造山演化具有由北向南迁移的特点。根据现代20°~25°间活跃弧陆碰撞带的地质及地球物理特征，将其划分为空间上的4个构造带：洋—洋俯冲；初期弧陆碰撞；成熟期弧陆碰撞；岛弧下陷/俯冲^[10]。

台湾造山带主要由3个地质构造单元组成，自西向东依次为：被动大陆边缘逆冲褶皱带(海岸平原—西部麓山带—雪山山脉)；弧陆碰撞前的俯冲增生楔(恒春半岛—中央山脉)；逆冲加附的火山岛弧—弧前盆地(海岸山脉)^{[9, 10]}。

海岸山脉呈NNE—SSW向展布，主要由晚上新世—更新世的残留弧前盆地和火山岛弧组成，残留弧前盆地由北往南依次为水琏、乐合、泰源和台东残留弧前盆地^{[9-12, 15, 16]}(图 2B)。它们与弧上盆地在岩性组成上类似，皆为大港口层深海相浊积岩所充填，整合覆盖于都銮山层之上，以半固结的砂岩—泥岩互层和泥岩为主，偶夹凝灰岩，地层厚度变化并不明显，平均厚度在2 km左右^[17]。火山岛弧在地层上均属于都銮山层，主要为以凝灰岩、火山岩流、岩脉、集块岩、火山角砾岩为主的中新世安山质火山岩^[18](图 2)。

图  2  乐合弧前盆地简易岩性柱状图(据文献[19]改编)

Figure 2.  Simple lithological column of Loho fore-arc basin (modified from reference [19])

海岸山脉在台湾的构造体系中为一独立的单位，是遭受强烈的剪切作用而向西逆冲的吕宋火山岛弧体系，其现今的构造特征受一系列近于NNE走向、向西逆冲的断层所控制^[18]。在构造形态上，海岸山脉被认为是由一系列叠瓦状排列的褶皱和断层组成一个逆冲褶皱带^[20-22]。海岸山脉的褶皱均向NNE延长，较台湾的主构造线更为偏东；主要断层均为向西纵向逆断层，大部分逆断层发生于褶皱形成之时或紧接褶皱之后的褶裂断层，连续的倾向相同的逆断层构成了典型的叠瓦状构造，次级断层有正断层和平移断层，均截切地层走向和纵向断层，明显区别于西部麓山带的逆冲褶皱带^[20]。

由于本文研究区域内的褶皱和断层具有良好的东西走向，因此，为了简化起见，可以用横向的二维模型来模拟整个弧前盆地的形成和演化过程。本文选用的模拟软件是FLAC(Fast Lagrangian Analysis for Continuum)6.0岩土工程软件，最早由美国Itasca Consulting Group. Inc开发应用于岩体力学计算^[23]。这个程序采用拉格朗日分析的二维显示有限差分法，有限差分法是用相邻节点之间的差商代替控制方程的微商，从而求解偏微分方程的问题转化成求解代数方程组的问题。拉格朗日分析法是一种利用拖带坐标系分析大变形问题的数值方法，并利用差分格式按时步积分求解，随着构形的不断变化，不断更新坐标的方法。

FLAC软件相比于其他现行的模拟方法有着明显的优势。FLAC计算中使用了“混合离散化”技术，更为精确和有效地模拟计算材料的塑形破坏和塑形流动。同时，显示计算方法更加节约内存空间和运算时间，便于实现非线性大变形的求解。FLAC软件包含了一个强大的内置编程语言——FISH语言(FLACish的缩写)，可以让用户定义新的变量和函数，从而实现对FLAC的扩充。当然，FLAC也有不足之处，当模拟大变形模型时，网格畸变后，重新划分网格易使不同材料边界发生均一化，对后续的变形产生一定的影响。

FLAC软件可以解决多方面的问题，其中包括：平面应力、平面应变几何模型问题、不连续面滑动和分离的问题、蠕变问题、热力学分析和动态分析等问题，如今已经被广泛地应用于褶皱与断层演化^[24-27]、岩石圈结构与形变^[28-31]以及流体-热液成矿^[32-34]的研究。

本文建立的岩石力学模型的尺度、力学参数、边界条件以及加载过程如下：

(1) 模型尺度 根据地震反射剖面解释(图 3)，整个实体模型由弧前盆地、盆地基底以及火山岛组成。为了建立接近真实尺度的数值模型，本文模拟对象水平长75 km，弧前盆地表层高度5 km，火山岛最高点为8 km。整个模型被剖分为180×31个网格，以确保模型不会因为网格过于稀疏而影响正常褶皱的形成和发展。

图  3  台湾弧前盆地的地震反射剖面(据文献[35]改编)

Figure 3.  The seismic reflection profile of Taiwan fore-arc basin (modified from reference [35])

(2) 岩石力学层 根据前人对于北吕宋海槽弧前盆地的地层分析^[12]，并考虑到计算存储空间和时间的关系，根据不同年代地层单元岩性之间的差异，将弧前盆地划分为5个力学层。为了简化起见，把弧前盆地恢复到最初始状态，将弧前盆地的地层平行于滑脱带分布。由于弧前盆地和火山岛弧具有不同的岩性组成，我们需要区别对待，把火山岛弧填充为安山质火山岩。模型中各层的岩石力学行为遵循摩尔—库伦弹塑性理论，变形主要是由弹性变形和塑性变形组成。需要注意的是，地质构造的模拟一般只有少量的限制数据，而且地质体不均一特性可能有很大的差异，不可能得到实际地质体的全部资料。所以，我们在设计过程中应该进行适当的运算，进行反复的参数调整。模型中不同的岩石力学层的具体参数见表 1，表 1中各层岩石力学参数的选取参照了前人的实测结果^[36]。

(3) 滑脱带 根据地震反射剖面解释(图 3)，滑脱带位于弧前盆地与火山弧基底之间。在FLAC中可以用于模拟断层、节理等不连续界面的方法主要为interface单元法，此功能可以模拟2个界面之间的滑动、分开等相互作用^[37]。在本文中将用interface单元法来模拟滑脱带，滑脱带的法向刚度Kn和剪切刚度Ks取值见表 2。

(4) 网格剖分 在模拟过程中，滑脱带附近的网格往往由于较大的应变量而变成畸形，从而导致计算终止。因此，在本文中将采用Auto-rezoning函数，其目的是当网格在模拟中变成畸形时，系统将自动分配现存的压力、速度、位移到一个新的、规律的网格中，以确保计算在大变形中能顺利进行。因为在重新分配中用到插值法，所以一些不确定性是必然存在的，这些不确定性会随着网格密度的增大，而逐渐减小。一般来说，对于模拟结果的影响可以忽略不计。

(5) 加载过程 在模型中，固定其右边界和底部边界，并在左侧滑脱带以上施加水平向右、0.4 m/时步的速度。模型的加载过程到总缩短率不超过50%或者因网格畸变而终止计算为止。

模拟结果显示(图 4)，随着变形的增加，褶皱向右侧发展，当缩短率为10%时，先在弧前盆地的断坡处形成了一个断弯褶皱，随后在弧前盆地的左侧形成一连串的完整的褶皱雏形。从模型中可以明显的看出，左侧的褶皱①波长较长，测得为10.4 km，右侧的断弯褶皱的波长较小，为4.5 km。为了方便研究，统一测量模型第2层砂岩的波幅，取第2层背斜最高点和向斜最低点的一半值，褶皱①、②的波幅分别为0.9和0.8 km，波长振幅比分别为12:1和6:1。

图  4  台湾弧前盆地逆冲构造的模型演化

Figure 4.  The evolution of thrust structure in Taiwan fore-arc basin

当缩短率达20%时，褶皱逐渐延伸到断弯褶皱②的右侧，同时，褶皱②的轴面发生向右的倾斜，倾角为50°，并且在模型的左侧具有了明显的褶皱形态特征，为向斜宽，背斜窄，宽度比值为1.4:1。褶皱①、②的波长分别缩短为7.8 km和3.9 km，波幅分别增大为1.25和1.5 km，波长振幅比为6:1和3:1。

当缩短量达到30%的时候，褶皱继续向右侧延伸，并沿着火山弧表面上冲，同时在褶皱②的左翼上形成了一个以向西倾斜的次级褶皱。我们发现此时的缩短变形主要集中在模型的右侧，而左侧的褶皱形态并未发生明显的变化。褶皱③的形态与缩短量为10%时褶皱②的形态相似，倾角为50°，波长为3.3 km，波幅为0.9 km，波长振幅比为4:1。在这之后的持续变形中，主要为褶皱的紧闭与隆起。

当缩短量达到40%的时候，在褶皱②和褶皱③之间存在一个明显西倾的次级褶皱，同时褶皱③的倾角增大到80°，近垂直状态，波长缩短为2.6 km，波幅增大为1.6 km，波长振幅比为1.6:1。

在当缩短量达到最大值50%的时候，各个褶皱的波长进一步减小，而波幅进一步增大，褶皱③的倾向逐渐由东倾转变为西倾，倾角为85°。

研究剪切应变增量(质点速度)的变化对于研究褶皱的形成演化同样具有重要的意义。正如人们所预料的那样，依次记录有限剪切应变增量可以显示岩层逐渐遭受的变形和破坏，这显示了地层发生局部变形的位置。

图 5描绘了剪切应变增量随着缩短量增加而发生的变化。我们发现，随着时间步的推移，增量剪切应变带的累计产生了有限应变，当缩短量为10%的时候，有限剪切应变主要集中在弧前盆地第一个向东倾斜的剪切带中和断弯褶皱前方的剪切带中；当缩短量增大到15%时，有限剪切增量在向东倾斜的剪切带中更为集中，从而形成了向东的逆冲推覆构造。当缩短量不断增大到20%和25%的时候，向西倾斜的剪切带中的剪切应变增量也在不断的累积，这就导致了在先前形成的叠瓦状褶皱的左翼形成了向西倾斜的次级褶皱，同时在弧前盆地最东侧褶皱的倾向逐渐发生反转，并最终演化为向西的逆冲构造。

图  5  有限剪切增量变化

Figure 5.  The change in finite shear strain increment

对比图 6、7模型演化过程中东西两侧褶皱的形态发现，弧前盆地的地层厚度在东西两侧存在明显的差异，呈现西厚东薄的形态，弧前盆地西侧的地层厚度平均为3.5 km，能干层为中间的砂岩，平均厚度为1.8 km，在砂岩上下均是较为软弱的泥岩；而在弧前盆地的东侧地层较薄，平均厚度只有1.7 km，能干层的平均厚度为0.6 km，两侧能干层厚度之比为3:1。

图  6  不同坡角的模型演化((A)坡角为0°，(B)坡角为3°)

Figure 6.  The evolution of different slope angle in the model((A) the slope angle is 0°, (B) the slope angle is 3°)

图  7  不同坡角的模型演化((A)坡角为5°；(B)坡角为10°；地层岩性同图 5)

Figure 7.  The evolution of different slope angle in the model ((A) the slope angle is 5°, (B) the slope angle is 10°, formation lithology with Fig. 5)

通过上述的模型结果分析发现，弧前盆地西侧的褶皱相比于东侧的褶皱具有更大的波长和更小的振幅。当缩短量为10%时，对比一下弧前盆地东西两侧的褶皱形态特征，以褶皱①、②为例，2个褶皱的轴面近似垂直，均为直立褶皱，同时能干层砂岩在发生正弦状弯曲的弯曲，而上下较为软弱的泥岩会被动的调整和适应这种弯曲形态，可以看到在模型图中并无较大的地形起伏。测得其波长分别为10.4和4.5 km，波幅分别为0.9和0.8 km，波长振幅比分别为12:1和6:1。当缩短量为20%时，褶皱①、②的轴面发生向东的倾斜，由直立褶皱转变为斜歪褶皱，同样泥岩会被动的调整和适应砂岩引起的弯曲，在向斜处具有更厚的泥岩层，相反在背斜处泥岩较薄。同时2个褶皱的波长分别缩短为7.8和3.9 km，波幅分别增大为1.25和1.5 km，波长振幅比分别为6：1和3：1，同样为西侧褶皱具有更大的波长和更小的振幅。

在图 6、7中，随着基底断坡倾角的不断增大，断弯褶皱的几何形态越发明显，其隆起幅度越大。当不存在断坡的时候可以发现，褶皱的演化是从西往东逐渐发育的，每一个褶皱的形态特征相似，并且为向斜宽，背斜窄。当坡角增大到3°的时候，整个褶皱的演化过程同样是自西向东逐渐发育的，与没有断坡时的演化方式类似，并无明显的差别。

当基底坡角继续增大到5°的时候，在断坡处形成了一个断弯褶皱，此处为最大起伏的隆起，而在断弯褶皱的左侧只是形成了一连串的褶皱雏形，波长较大，而振幅相对较小。当坡角增大为10°的时候，我们对比模拟图发现，缩短量为10%时，在缩短变形的过程中优先形成的是断坡处的断弯褶皱，相比较于5°坡角的时候，断弯褶皱明显变得更加紧闭，其隆起幅度也相应增大。而在断弯褶皱的左侧并没有明显的褶皱形成，相反在其右侧发育一个小褶皱；当缩短量达15%的时候，随着挤压变形的继续，褶皱从断弯褶皱处向右逐渐延伸，而在断弯褶皱的左侧有稍微的隆起，但并无明显的褶皱形态；当缩短量达到20%的时候，褶皱从断弯褶皱处逐渐延伸至弧前盆地的最右侧，断弯褶皱的倾角变小，同时变得更加紧闭。

根据数值模拟的演化结果可以发现，台湾弧前盆地的构造反转是一个剪切增量累积的过程，构造反转开始于成熟的弧陆碰撞时期，即更新世早期，由局部的构造反转最终演化为整体的构造反转，这与Jacques等的沙盘模型结果相一致^[38]。并且通过观察图 2中的构造剖面图可以发现，图中同时存在2个方向的大角度逆断层，这也与我们的模拟结果相吻合，也可以通过模拟结果对Huang提出的弧前盆地的演化模型进行补充^[12](图 8)：

图  8  台湾弧前盆地的构造演化模型(据文献[12]改编)

Figure 8.  Tectonic evolution model of Taiwan fore-arc basin (modified from reference [12])

(1) 在初期弧陆碰撞早期(3.7~3.2 Ma)，即晚中新世时期，菲律宾海板块沿着马尼拉海沟俯冲于欧亚板块之下，此时由于强烈的弧陆碰撞，台湾中央山脉逐渐抬升隆起。

(2) 在初期弧陆碰撞的后期(3~2 Ma)，北吕宋海槽弧前盆地受到向东的顺层挤压力，形成了一系列向东倾斜的逆冲叠瓦状褶皱和断层，并逐渐向东延伸至最东侧的北吕宋火山岛弧，演化形成了花东海脊，并在其东侧发育年轻的残留弧前盆地。

(3) 在成熟的弧陆碰撞时期(< 2 Ma)，变形的弧前盆地继续遭受着强烈的挤压作用，随着弧前盆地中有限剪切应变增量不断的累计，在先前形成的叠瓦状褶皱的左翼发育了向西倾斜的次级褶皱，同时最东侧部分褶皱的倾向发生反转，最终演化为向西的逆冲构造，形成了如今的海岸山脉。

海岸山脉作为一个独立的地质单元，是在弧陆碰撞后期距今2 Ma才向西逆冲的火山岛弧体系，其构造特征受到弧陆碰撞的控制。在整个弧陆碰撞过程中，弧前盆地在受到向东水平挤压的同时，西侧的北吕宋火山岛弧提供了强烈的支撑作用，才会导致构造反转的发生。

根据北吕宋海槽北部的剖面图可以看出(图 9)，弧前盆地在不同的位置，其地层的厚度是有较大差异的，而地层厚度与褶皱的形成演化有着密切的联系。在叠瓦状褶皱形成的过程中，褶皱的初始波长受褶皱地层厚度决定，表现为褶皱波长与地层厚度成正比，这与其他学者利用不同的研究方法总结归纳获得的结论相同^{[39, 40]}。

图  9  北吕宋海槽北部南北向构造剖面(据文献[41]改编)

Figure 9.  A N—S line near the northern half of the North Luzon Trough (modified from reference [41])

而且不同能干层的褶皱形态是不同的，台湾弧前盆地的宽度是远远大于地层厚度的，同时地层厚度呈现东厚西薄的形态。所以，在弧前盆地西侧的褶皱波长较大，地层向东逐渐减薄，褶皱的波长逐渐减小。

同时，岩层间的力学性质差异在褶皱形成中也起着主导作用。弧前盆地在岩性组成上，以砂岩—泥岩互层和泥岩为主，岩系中砂岩和泥岩的力学性质不一致，则在顺层挤压作用下，强硬层即砂岩层就会失稳而发生正弦曲线状弯曲；相对软的层即泥岩层作为介质，在均匀压扁的同时被动地调整和适应由强硬层引起的弯曲形态。进一步挤压下，强硬层的褶皱变得越紧闭，可使翼部被压扁而成IC型褶皱。

由于岛弧基底的起伏不平，同样导致了不同位置上坡角的差异。对比不同坡角的演化模型(图 6，7)可以看出，当坡角为0°和2°的时候，褶皱的演化为由西向东规律性发育；当坡角增加到5°的时候，由于坡角的影响，优先在断坡处发育断弯褶皱，随着缩短量的增加，其他位置的褶皱才逐步发育，同时也是由西向东逐渐演化；当坡角达到10°的时候，这一现象更加明显，优先形成的断弯的波长更小，而振幅更大，由于东侧的地层更薄，所以褶皱向东延伸的速度比坡角为5°时明显要快。

通过4个模型的对比可以看出，坡角对褶皱发育的次序具有明显的控制作用，坡角达到一定值时，在断坡处优先发育断弯褶皱。而临界坡角的大小与地层的厚度成正比关系，即地层越厚，临界坡角越大。

(1) 在台湾弧前盆地中，逆冲构造的反转是由于向西的剪切带中应变增量的增加形成了向西倾斜的剪切带，同时使得在先前形成的褶皱的左翼上发育了向西倾斜的次级褶皱，发生了局部的构造反转，弧前盆地最东侧的部分褶皱的倾向也逐渐由东向西，最终演化成了以向西为主的逆冲构造引起的。

(2) 地层的厚度对于褶皱的形态具有控制作用，褶皱波长与地层厚度成正比，地层越厚，波长越大，同时软弱层会被动地调整和适应由强硬层引起的弯曲形态。

(3) 岛弧基底的坡角控制了弧前盆地中褶皱发育演化的次序，当坡角达到一定值之后，在断坡处的褶皱优先发育，并向地层薄的一侧转播地更快，而这一临界坡角与地层的厚度成正比。