地壳形变是地壳在力学、热学等复杂环境作用下表现出来的综合响应,研究地壳形变特征及其演化规律是认识地壳运动及其伴生地质灾害的一个重要途径。随着计算机科学的日益发展,数值模拟技术已经成为地壳动力学机理研究的重要手段(Moreno et al,2010; Hergert,Heidbach,2010; Li et al,2014; Zhu,Zhang,2010; He et al,2013; Lu et al,2011; 陈连旺等,1999,2008; 邵志刚等,2008; 李玉江等,2014)。在复杂地质构造作用下,断层作为岩石圈岩石破裂形成的特殊产物,在现今地壳形变特征以及地震的孕育、发生中有重要的作用。在探索地震发生孕育的过程中,国内外科学家从不同角度阐述了断层物理机制(Sibson,1977; 马胜利,本利彦,1995)。但由于地质构造运动的复杂性,现在还没有一种行之有效的办法可以模拟真实地壳环境中断层的形成、发展变化。因此通过数值模拟解释现今地壳形变特征时,往往需要结合野外考察、深地震反射、地震精定位、震源机制解等多学科方法给出的断层产状及断裂带断层泥等做出一定程度的简化,进而构建断层物理模型。

在长期的实践应用过程中,根据物理模型及具体问题的研究,这些简化的断层模型可以归纳为:劈节点模型、块体模型、介质弱化模型以及摩擦接触模型(和平等,2011)。劈节点模型本质上是在预设断层位错的情况下模拟断层,因而在同震或震后效应研究中具有极大的优势,得到了广泛应用; 块体模型本质基础是Shi(2007)提出的DDA(Discontinuous Deformation Analysis)算法,它将地球动力过程看作为一系列块体相互作用的过程,断层其实就是块体与块体之间的边界。介质弱化模型和摩擦接触模型作为2种原理简单且便于在一般软件中设置的模型,在区域地壳动力学形变特征模拟及断层活动状态模拟中得到大量应用。由于断层的介质力学参数对区域应力场以及形变场均有一定程度的影响(李玉江等,2010),因此在模拟过程中,介质弱化模型是将断层视为有一定厚度的弱介质层,通常通过降低该介质层弹性模量来模拟断层对区域地壳动力过程的影响(Lu et al,2011,2012; 祝爱玉等,2015a; 胡勐乾等,2014; 杨兴悦等,2013); 摩擦接触模型则是直接在创建好的模型基础上给断层面赋摩擦系数,进而实现断层两盘的相对运动(Hergert,Heidbach,2010; Zhu,Zhang,2010; He et al,2013; 李玉江等,2014; 祝爱玉等,2015b)。但也存在一些问题:①对于介质弱化模型来说,其本质上是通过降低断层介质弹性模量来近似模拟断层两盘的非连续变形。尽管野外观测到断层带内有糜棱岩或断层泥等弱力学属性材料填充,但这种弱力学属性的填充物已经在长期的构造作用下与周边介质处于力平衡状态。因此,在构建有限元模型时直接赋予该填充物一个较低的弹性模量,以期获得断层两盘的运动差异和断层附近应力状态的同时,往往忽略了降低弹性模量,也使得该填充层内部、垂直于断层走向的力平衡破坏,在断层内部垂直于断层方向产生新的附加力。这可能影响模拟结果的可靠性。此外,实验室关于断层泥的弱力学属性对断层行为的影响也表现在剪切摩擦行为上(马胜利,1986)。在各向同性弹性理论前提下还可能会出现不合理的负泊松比的情况(董培育,石耀霖,2013)。②摩擦接触模型可以很好地刻画断层两盘相对运动的非连续变形,相比介质弱化模型更接近于真实断层状态。但摩擦接触计算是个强非线性问题,涉及到收敛与否以及计算量等问题(李玉江等,2009),特别是在研究区域内断层较多的情况下,如果将每条断层都构建为摩擦基础模型,收敛问题将变得更加紧迫,为此研究中还可以将重点断层考虑为摩擦接触,而将次要断层处理为弱化带模型(Li et al,2015),这一简化处理提高了计算效率,但又增加了弱化介质模型自身的问题。

因此,如何做到既保证计算效率、避免介质弱化模型可能带来的不真实信息,又能快速合理模拟断层两盘相对运动是地壳动力学数值模拟的一个基础性工作。本文在介质弱化带模型基础上引入正交各向异性理论,用正交各向异性本构方程替代传统弱化带中的线弹性本构方程,实现了弹性模量与剪切模量相互独立; 同时考虑到实际建模过程中可能遇到不同断层走向、倾角各不相同的情况,通过坐标转换给出了针对任意走向、倾角的断层本构关系。最后通过数值试验对比分析了介质弱化模型、摩擦接触模型以及正交各向异性断层模型。

假定断层带为夹在块体之间的类似于层合板的一类特殊介质:其在垂直于断层方向力学性质与周边块体相近,不易产生大挤压或拉张变形; 而在平行于断层运动方向上的剪切模量较小,能够产生较大的剪切变形。则根据复合材料力学(沈观林,胡更开,2006)基础知识,断层内应力应变关系在其3个主方向上满足如下关系:

ε_6×1=K_6×6*σ_6×1(1)

式中:ε=[e₁₁ e₂₂ e₃₃ e₁₂ e₂₃ e₁₃]^T,σ=[s₁₁ s₂₂ s₃₃ s₁₂ s₂₃ s₁₃]^T分别为3个主方向的应变和应力分量; K为柔度矩阵; E₁,E₂,E₃,V₁₂,V₂₃,V₁₃,G₁₂,G₂₃,G₁₃分别为其3个主方向的工程弹性常数。K可表示为:

K=[1/(E₁)-(v₁₂)/(E₂)-(v₁₃)/(E₃)0 0 0

-(v₂₁)/(E₁)1/(E₂)-(v₂₃)/(E₃)0 0 0

-(v₃₁)/(E₁)-(v₃₂)/(E₂)1/(E₃)0 0 0

0 0 0 1/(G₁₂)0 0

0 0 0 0 1/(G₂₃)0

0 0 0 0 0 1/(G₁₃)](2)

且:(v_ij)/(E_j)=(v_ji)/(E_i)(i,j=1,2,3,且i≠j)。

由于本文重点考虑断层介质在走向和倾向的剪切行为,其他方向上的介质力学性质与周边介质一致。为此将式(2)进一步简化为:

K=[1/E -v/E -v/E 0 0 0

-v/E 1/E -v/E 0 0 0

-v/E -v/E 1/E 0 0 0

0 0 0(2(1+v))/E 0 0

0 0 0 0 1/(G_slip)0

0 0 0 0 0 1/(G_slip)](3)

式中:E,v分别为周边介质弹性模量和泊松比; G_slip为断层在剪切方向上的剪切模量,独立于周边介质参数E,v。

在实际操作过程中,由于区域断层较多且其走向和倾角各不相同,因此需要将其主方向下的应力应变矩阵旋转变换到地理坐标(x^-,y^-,z^-)下。图1为断层面主方向(1^-,2^-,3^-)以及地理坐标系(x^-,y^-,z^-)相对关系,其中:1^-,2^-,3^-分别为断层走向(面向走向,右手在上盘)、倾向(逆冲为正)以及张拉方向(张为正); x^-,y^-,z^-分别为地理东、北向以及垂直地表向上的方向。α为断层走向,β为断层倾角。令断层面应变和应力分量在地理坐标系(x^-,y^-,z^-)下分别为:

e=[e_xx e_yy e_zz e_xy e_yz e_xz]^T

s=[s_xx s_yy s_zz s_xy s_yz s_xz]^T(4)

由图1b可知,将地理坐标系下的应力张量绕z轴逆时针旋转π/2-α,然后再绕主轴1逆时针旋转β便可得到主方向下的应力张量。则:

[σ₁₁ σ₁₂ σ₁₃

σ₁₂ σ₂₂ σ₂₃

σ₁₃ σ₂₃ σ₃₃]=T·[σ_xx σ_xy σ_xz

σ_xy σ_yy σ_yz

σ_xz σ_yz σ_zz]·T'(5)

式中:T=A*B*[1 0 0

0 1 0

0 0 1]=[T₁₁ T₁₂ T₁₃

T₁₂ T₂₂ T₂₃

T₁₃ T₂₃ T₃₃]

A=[1 0 0

0 cos(β)sin(β)

0 -sin(β)cos(β)]

B=[cos(90°-α)sin(90°-α)0

-sin(90°-α)cos(90°-α)0

0 0 1]

结合式(1)和(4),式(5)可以转化为:

σ_6×1=F_6×6·s_6×1(6)

式中:F为坐标转换矩阵,可表示为:

F=[T₁₁T₁₁ T₁₂T₁₂ T₁₃T₁₃ 2T₁₁T₁₂ 2T₁₂T₁₃ 2T₁₁T₁₃

T₂₁T₂₁ T₂₂T₂₂ T₂₃T₂₃ 2T₂₁T₂₂ 2T₂₂T₂₃ 2T₂₁T₂₃

T₃₁T₃₁ T₃₂T₃₂ T₃₃T₃₃ 2T₃₁T₃₂ 2T₃₂T₃₃ 2T₃₁T₃₃

T₁₁T₂₁ T₁₂T₂₂ T₁₃T₂₃ T₂₁T₁₂₊T₂₂T₁₁ T₂₂T₁₃₊T₂₃T₁₂ T₂₁T₁₃₊T₂₃T₁₁

T₂₁T₃₁ T₂₂T₃₂ T₂₃T₃₃ T₃₁T₂₂₊T₃₂T₂₁ T₃₃T₂₂₊T₃₂T₂₃ T₃₁T₂₃₊T₃₃T₂₁

T₃₁T₁₁ T₃₂T₁₂ T₃₃T₁₃ T₃₁T₁₂₊T₃₂T₁₁ T₃₃T₁₂₊T₃₂T₁₃ T₃₃T₁₁₊T₃₁T₁₃]

同理可得:

ε_6×1=F_6×6·e_6×1(7)

将式(6)和式(7)带入式(1),可得:

F·e=K·F·s(8)

即:

e=D·s=F^-1·K·F·s(9)

式中:D为正交各向异性材料在地理坐标系(x^-、y^-、z^-)下的柔度矩阵,且:

D=F^-1·K·F=[D₁₁

D₂₁ D₂₂ 对

D₃₁ D₃₂ D₃₃ 称

D₄₁ D₄₂ D₄₃ D₄₄

D₅₁ D₅₂ D₅₃ D₅₄ D₅₅

D₆₁ D₆₂ D₆₃ D₆₄ D₆₅ D₆₆](10)

式中:柔度矩阵D为21个独立参数的对称矩阵,因此在实际模拟计算过程中,只需修改正交各向异性本够关系柔度矩阵K在其相应主方向下对应滑动方向的剪切模量(式(3)中的G_slip),结合坐标转换矩阵F,代入式(10)便可生成断层有限元模型的柔度矩阵,进而可以构建有限元计算的整体刚度矩阵,实现垂直断层方向无大变形,平行断层滑动方向有大变形滑动的断层运动状态。为方便控制,本文令断层走向和倾向剪切模量与周边弹性介质剪切模量之比为λ∈[0,1],即G_slip=λ·E/(2(1+v))。

图1 断层变形模式示意图(a)以及断层面主方向与地理坐标的关系(b)
Fig.1 Fault deformation mode(a)and the relation-ship between geographic coordinate system and principal direction of fault deformation(b)

近年来,随着GPS地壳速度场观测数据的丰富,基于GPS约束的地壳动力学数值模拟研究日渐深入(Hergert,Heidbach,2010; He et al,2013; Lei et al,2010; Cho,Kuwahara,2013; 龙小刚,朱守彪,2015; 刘峡等,2010)。 一般情况下,研究人员通常都会通过一些途径选定有限元模型的边界条件,模拟计算给出模型内部的速度场模拟结果,然后通过对比模拟结果与观测结果之间的残差来确定模型是否可靠,参数选择是否合理,在大量比较之后给出合理模型的基础上分析区域的动力学机制及应力状态等,这就要求我们首先要保证所采取的模型是合理的。介质弱化模型作为最为简单的断层模型,其建模方便、不涉及非线性问题且不会存在难于收敛的情况,在一些研究中得到大量的应用。因此在应用介质弱化模型的情况下,在对比模拟地壳运动速度场与GPS观测速度场残差的同时,如何评价类似于图4a、d以及图5d的形变场畸变对模拟结果的影响也是不容忽视的问题。在不可避免使用介质弱化断层模型的时候,如能进一步修改其有限元刚度矩阵,使其剪切模量与弹性模量相互独立,对于提高模拟结果的可靠性具有一定的意义。

本文在介质弱化带模型基础上,引入正交各向异性理论,同时考虑到实际建模过程中可能遇到不同断层的走向、倾角各不相同的情况,本文通过坐标转换给出了针对任意走向、倾角的正交各向异性断层本构关系。最后通过数值试验对比分析了介质弱化模型、摩擦接触模型以及正交各向异性断层模型在震间地壳速度长图像模拟中的差异,结果显示:(1)无论断层力学性质强弱,传统的介质弱化模型都会对地表地壳运动位移场产生一定程度的影响,断层力学性质(对应于图2不同断层模型中断层的力学参数,如弱化介质模型的弹性模量、各向异性模型的剪切模量以及摩擦接触模型的摩擦系数)越弱,影响越显著,这主要由于是降低断层弹性模量,使得断层内垂直于断层方向的力平衡被打破,进而影响了断层附近位移场图像的空间分布,这在模拟中同时也会对断层附近的应力场和应变场产生一定程度的影响;(2)从位移场图像以及应力场图像的空间分布来看,正交各向异性本构关系与摩擦接触模型给出的模拟结果基本一致,因此在区域构造复杂的情况下,可以考虑使用各向异性本构关系来表征断层的变形机制。

非常感谢陈连旺研究员的悉心指导,同时也感谢中国地震局访问学者计划的支持和地壳应力研究所提供良好的办公环境和软件支持!部分图件采用GMT(Wessel et al,2013)绘制,特此表示感谢!