数值模拟方法
如前所述，基于镜质组化学反应动力学的煤热演化史数值模拟方法经历了简单函数关系模式、受热时间－经验法模式、反应活化能－温度函数模式和平行反应化学动力学模式四个发展阶段（秦勇等，1995）。其中，目前在化石能源地质界最为常用的是TTI法、LOM法和Tissot法，近年建立起来的EASY%Ro方法也已成为数值模拟的重要发展方向。（一）时间－温度指数（TTI）法这种方法的基本思想由Lopatin（1971）首先提出，后经Waples（1980）根据31个盆地402个样品的实测资料加以完善，建立起了其数值模拟的方法体系，故又称为Lopatin-Waples法。TTI法遵循化学反应动力学的基本法则来衡量温度和时间两个因素对煤化作用的贡献，即煤化作用速度随受热温度的增高而增大，温度每增加10℃煤化反应速度增加一倍。由此，可将TTI数学模式用连续函数加以表示：山西南部煤化作用及其古地热系统：兼论煤化作用的控气地质机理T=Ts＋G·D式中：T——煤层受热的古地温温度，℃；Ts—－古地表温度，℃；G——古地温梯度，℃/m;D—煤层埋藏深度，m。上述连续函数模式不便于常规运算，故在实际工作中通常采用分段积分方式，根据如下方程进行计算：TTI＝∑γn·△t式中：γ——温度效应因子，考虑到温度每增高10℃煤级量值增加一倍的关系，取r=2;n——温度指数，其大小取决于温度（表5-1）；△t——煤层在某一温度段的受热时间，Ma。表5-1 TTI数学模式中温度指数（n）与煤层受热温度段的对应关系注：温度每增（减）10℃，温度指数增（减）1。TTI值本身可以作为煤级或有机质成熟度指标，Waples（1980）也提供了TTI值与镜质组最大油浸反射率的关系（表5-2）。此外，Wood（1988）通过回归拟合，将二者之间关系采用下式表示：lg%Ro＝—0.006（lgTTI）3+0.042（lgTTI）2+0.162（lgTTI）—0.397表5-2 时间温度指数（TTI）与镜质组反射率（Ro）的相关关系应予指出，TTI法尽管几十年来被化石能源地质界所广泛采用，但仍存在某些明显的不足有待于改进，国内外学者为此做过一定努力。例如，温度每增加10℃煤化速率加快一倍的假设有其局限性，只能在活化能10～25kcal/mol（相当于20～160℃的温度范围，1kcal=4186.8J，下同）区间是可靠的（Magoon,1983）。再如，沉积有机质热降解的反应活化能是受热温度的函数，随成熟度的增加而逐步加大，不同类型有机质达到相同成熟度所需的活化能也有差异，而TTI模式将整个煤化过程中的反应活化能作为常量，显然与实际情况有出入。（二）有机成熟度水平（LOM）法该方法由Hood（1975）建立，后经Bostick等（1979）完善而成为一种国内外广泛采用的数值模拟方法，故又被称为Hood-Bostick方法。Hood等考虑到总体反应活化能随受热温度而增高的变化趋势（18～33kcal/mo1），采用镜质组油浸反射率标定有机质成熟度水平（LOM）（图5-1）。在模式中采用了有效受热时间的概念，即温度不低于最高受热温度15℃范围内的受热时间，建立起有机质成熟度、温度及有效受热时间之间的相互关系（图5-2）。图5-1 有机质成熟度水平（LOM）与其它有机质成熟度指标的关系（引自周中毅，1990）Ⅰ—褐煤；Ⅱ1—亚烟煤；Ⅱ2—高挥发性烟煤；Ⅱ3—中挥发性烟煤；Ⅱ4—低挥发性烟煤；Ⅲ1—半无烟煤；Ⅲ2—无烟煤1Btu=1055.06J在实际工作中，LOM法采用图解方法求取有机质成熟度或镜质组反射率（图5-2）。国内外应用效果表明，利用这种方法推算出的盆地古地温温度是比较可靠的，与由综合研究所得出的结果以及盆地实际情况较为吻合（Vote,1981；周中毅等，1983,1984,1985）。
数值模拟方面
虽然该区的研究很丰富，但是在动力学数值模拟等方面，大部分都是基于对地震的关注，着眼于地壳等浅部的应力分布、变化以及地壳水平运动、形变等问题，并未从地质的角度研究构造形成的机制、深层对浅部形态、应力分布等的控制，以及时空演化、作用过程。在国内，数值模拟的研究工作在西部，尤其是青藏高原处做的比较多（张贵宾等，2001；傅容珊等，2000），但是涉及我们要模拟的鄂尔多斯和华北的区域的数值模拟相对较少（陈连旺，2001；傅容珊，2000；张东宁，1999；等）。但是陈连旺只涉及到地壳，给出了鄂尔多斯隆起东缘以东的华北地区三维构造应力场的模型，探讨的是与地震的关系。他所设计的三维模型思路是很值得借鉴的。张东宁涉及的是动力边界模拟。傅容珊等用数值模拟的手段研究了印度板块和欧亚板块碰撞后40 Ma年以来青藏高原的挤压、隆升、东亚大陆形变及应力场的演化过程。通过应力场格局与地震、地质观测相验证，现代水平形变速度与GPS观测相结合，证明印度板块和欧亚板块的碰撞、挤压是构成岩石层水平形变和应力场格局的主要驱动力。从目前国际岩石圈研究发展不难看出，作为一种现代方法技术，动力学数值模拟正在并将继续起到重大作用，带动了岩石圈研究的发展。数值模拟技术本身的不断发展，各种新技术的运用如GPS，使得动力学、数值模拟的研究空间和效率得以拓展。Martin H.P.Bott教授是一直从事有限元构造模拟的国际知名构造地球物理学家之一。他与他的科研组的研究领域涉及地球动力学模拟的多个方面。他利用有限元分析模拟了大陆裂谷的应力分布，研究大陆裂谷的形成机制（Bott,1992）。1984年研究了岩石圈中构造应力的起源及其与构造活动的关系（Bott and Kusznir,1984）。1989年利用粘弹有限元模拟技术研究了俯冲带处仰冲板块和俯冲板块中的应力状态以及俯冲带板块边界力（Bott,1989）。1990年模拟了与碰撞山脉相关的应力分布和板块边界力，研究了造山带和边缘盆地的产生机制（Bott,1990）。1991年模拟了由洋脊推动引起的在被动大陆边缘两侧海洋岩石圈和大陆中的应力分布及板块边界力（Bott,1991）。1992年模拟了把南极洲分成东西两部分的隆升山脉和与之相邻的盆地沉降（Bott and Stern,1992）。1993年Bott总结自己的一系列有关论文，研究了上地幔密度异常如何在岩石圈中产生构造应力从而驱动包括板块运动在内的绝大多数构造活动（Bott,1993）。Bott结合重力异常研究，利用粘弹塑性有限元方法模拟了这种小尺度局部地壳均衡隆升机制以及这其中断裂的作用，认为这种均衡隆升是下地壳牛顿流体或幂律流体的蠕变引起的（Bott,1999）。Zhang和Bott（2000）模拟了上地壳中由平面高角度逆断层界定的非对称盆地的演化机制，为解释既窄又深的欧洲型前陆盆地（不同于较宽的西加拿大型）是如何形成的提供了一种可能的理论。张贵宾、高锐和肖序常（2001）模拟研究了青藏高原横过西昆仑和塔里木结合带的山隆盆降动力学机制。通过这一系列的研究、模拟，形成了比较系统的具有针对性的粘弹塑性模拟技术FEVPLIB。本文将利用PEVPPLIB技术，基于上地壳、下地壳、岩石圈上地幔和软流圈具有不同的流变性的粘弹模型，模拟研究华北地区新生代以来的地质、地貌、岩石圈三维结构图像的形成机制和演化问题。
