地壳流变结构控制作用下的龙门山断裂带地震发生机理

第５７卷第２期地　球　物　理　学　报Ｖｏｌ．５７，Ｎｏ．２柳畅，石耀霖，朱伯靖等．２０１４．地壳流变结构控制作用下的龙门山断裂带地震发生机理．地球物理学报，５７（２）：４０４－４１８，ｄｏｉ：１０．６０３８／ｃｊｇ２０１４０２０７．Ｌｉｕ　Ｃ，Ｓｈｉ　Ｙ　Ｌ，Ｚｈｕ　Ｂ　Ｊ，ｅｔ　ａｌ．２０１４．Ｃｒｕｓｔａｌ　ｒｈｅｏｌｏｇｙ　ｃｏｎｔｒｏｌ　ｏｎ　ｔｈｅ　ｍｅｃｈａｎｉｓｍ　ｏｆ　ｔｈｅ　ｅａｒｔｈｑｕａｋｅ　ｇｅｎｅｒａｔｉｏｎ　ａｔ　ｔｈｅ　ＬｏｎｇｍｅｎＳｈａｎ　ｆａｕｌｔ．Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．（ｉｎ　Ｃｈｉｎｅｓｅ），５７（２）：４０４－４１８，ｄｏｉ：１０．６０３８／ｃｊｇ２０１４０２０７．地壳流变结构控制作用下的龙门山断裂带地震发生机理柳畅１，２，３，石耀霖３＊，朱伯靖３，程惠红３，杨小林４１华中科技大学，物理学院，地球物理研究所，武汉　４３００７４２Ｌａｂｏｒａｔｏｉｒｅ　ｄｅ　Ｇéｏｌｏｇｉｅ，ＣＮＲＳ－ＵＭＲ８５３８，éｃｏｌｅ　Ｎｏｒｍａｌｅ　Ｓｕｐéｒｉｅｕｒｅ，Ｐａｒｉｓ　７５２３１３中国科学院大学，中国科学院计算地球动力学重点实验室，北京　１０００４９４陕西省地震局，西安　７１００６８摘要　青藏高原东缘低地形变速率的龙门山断裂带上相继发生了２００８汶川Ｍｗ７．９级地震和２０１３芦山Ｍｗ６．６级地震．地震勘探与震源定位结果揭示了龙门山区域地震空间分布特征：纵向上，龙门山断裂带这两次地震主震均发生在龙门山断裂带上地壳的底部（１４～１９ｋｍ），绝大部分余震均发生在上地壳范围（５～２５ｋｍ），而在其中、下地壳深度范围内鲜见余震发生；横向上，地震（Ｍｗ＞３）在龙门山断裂带青藏高原一侧密集分布且曾有大震发生，而四川盆地地震稀少（Ｍｗ＞３）．为探讨龙门山断裂带地震发生机理，并解释以上龙门山区域地震空间分布特征，本文建立了龙门山断裂带西南段跨芦山地震震中区域的四种不同流变结构的龙门山断裂带三维岩石圈模型，以地表ＧＰＳ观测资料为约束边界条件，数值模拟龙门山断裂带岩石圈在数千年以上长期匀速构造挤压作用下的应力积累特征，探讨了地壳分层流变性质对地壳应力积累的影响，分析了该区域地震空间分布与构造应力积累速率的关系．计算结果表明：该区域在数千年的应力积累过程中，脆性上地壳中应力表现近于恒定值的线性增长趋势，龙门山断裂带上地壳底部出现应力集中积累现象，这一应力集中现象可以解释龙门山断裂带汶川地震与芦山地震主震的发生，及其大部分余震在脆性上地壳中的触发；青藏高原一侧上地壳应力积累速率远远高于四川盆地的应力积累速率，这一应力积累分布现象可以解释龙门山区域青藏高原一侧地震密集而四川盆地地震稀少的地震空间分布特征；通过比较不同流变结构模型中的应力积累状态，认为导致这一应力积累空间分布状态的重要控制因素在于青藏高原中、下地壳较低的黏滞系数与四川盆地中、下地壳较高的黏滞系数的差异．在柔性的中、下地壳内，应力增长近于指数形式，稳定状态之后其应力增长速率近于零，构造应力积累难以达到岩石破裂强度，因而鲜见地震发生．地壳各层位的应力增长率差异与地震成层分布的现象共同揭示了龙门山区域岩石圈分层流变结构：脆性上地壳、韧性中、下地壳（青藏高原一侧较弱，四川盆地一侧较强）、韧性岩石圈上地幔．关键词　龙门山断裂带，汶川地震，芦山地震，应力集中，黏性差异，Ｍｏｈｏ面突变ｄｏｉ：１０．６０３８／ｃｊｇ２０１４０２０７中图分类号　Ｐ３１５，Ｐ３１３收稿日期２０１３－１１－２１，２０１４－０１－１５收修定稿Ｃｒｕｓｔａｌ　ｒｈｅｏｌｏｇｙ　ｃｏｎｔｒｏｌ　ｏｎ　ｔｈｅ　ｍｅｃｈａｎｉｓｍ　ｏｆ　ｔｈｅ　ｅａｒｔｈｑｕａｋｅ　ｇｅｎｅｒａｔｉｏｎａｔ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔＬＩＵ　Ｃｈａｎｇ　１，２，３，ＳＨＩ　Ｙａｏ－Ｌｉｎ３＊，ＺＨＵ　Ｂｏ－Ｊｉｎｇ３，ＣＨＥＮＧ　Ｈｕｉ－Ｈｏｎｇ３，ＹＡＮＧ　Ｘｉａｏ－Ｌｉｎ４１　Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｇｅｏｐｈｙｓｉｃｓ，Ｓｃｈｏｏｌ　ｏｆ　Ｐｈｙｓｉｃｓ，Ｈｕａｚｈｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｓｃｉｅｎｃｅ　ａｎｄ　Ｔｅｃｈｎｏｌｏｇｙ，Ｗｕｈａｎ　４３００７４，Ｃｈｉｎａ２　Ｌａｂｏｒａｔｏｉｒｅ　ｄｅ　Ｇéｏｌｏｇｉｅ，ＣＮＲＳ－ＵＭＲ８５３８，éｃｏｌｅ　Ｎｏｒｍａｌｅ　Ｓｕｐéｒｉｅｕｒｅ，Ｐａｒｉｓ　７５２３１，Ｆｒａｎｃｅ３　Ｌａｂｏｒａｔｏｒｙ　ｏｆ　Ｃｏｍｐｕｔａｔｉｏｎａｌ　Ｇｅｏｄｙｎａｍｉｃｓ，Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　ｔｈｅ　Ｃｈｉｎｅｓｅ　Ａｃａｄｅｍｙ　ｏｆ　Ｓｃｉｅｎｃｅｓ，Ｂｅｉｊｉｎｇ１０００４９，Ｃｈｉｎａ４　Ｅａｒｔｈｑｕａｋｅ　Ａｄｍｉｎｉｓｔｒａｔｉｏｎ　ｏｆ　Ｓｈａａｎｘｉ　Ｐｒｏｖｉｎｃｅ，Ｘｉ′ａｎ　７１００６８，Ｃｈｉｎａ基金项目　国家科技支撑计划《地震分析预测若干实用技术研究》项目（２０１２ＢＡＫ１９Ｂ０３－５）和国土资源部深部探测项目（ＳｉｎｏＰｒｏｂｅ－０７）资助．作者简介　柳畅，博士，主要从事岩石圈动力学，地震力学方面的研究．Ｅ－ｍａｉｌ：ｌｉｕ＠ｇｅｏｌｏｇｉｅ．ｅｎｓ．ｆｒ＊通讯作者　石耀霖，教授，主要从事地球动力学方面的研究．Ｅ－ｍａｉｌ：ｓｈｉｙｌ＠ｕｃａｓ．ａｃ．ｃｎ　２期柳畅等：地壳流变结构控制作用下的龙门山断裂带地震发生机理Ａｂｓｔｒａｃｔ　Ｉｎ　ｔｈｅ　ｒｅｃｅｎｔ　５ｙｅａｒｓ　ｔｈｅ　２００８　Ｍｗ７．９Ｗｅｎｃｈｕａｎ　ｅａｒｔｈｑｕａｋｅ　ａｎｄ　２０１３　Ｍｗ６．６Ｌｕｓｈａｎｅａｒｔｈｑｕａｋｅ　ｏｃｃｕｒｒｅｄ　ｉｎ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ　ａｔ　ｔｈｅ　ｅａｓｔｅｒｎ　ｍａｒｇｉｎ　ｏｆ　ｔｈｅ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ，ｗｈｅｒｅ　ｔｈｅ　ｃｒｕｓｔ　ｃｏｎｖｅｒｇｅｎｃｅ　ｒａｔｅ　ｉｓ　ｌｏｗ（～３ｍｍ／ａ）．Ｉｎ　ｐｒｅｖｉｏｕｓ　ｒｅｓｅａｒｃｈｅｓ　ｉｔ　ｉｓ　ｓｈｏｗｎ　ｔｈｅ　ｍａｉｎｓｈｏｃｋｓ　ｏｆ　ｔｈｅｓｅ　ｔｗｏ　ｅａｒｔｈｑｕａｋｅｓ　ａｒｅ　ｌｏｃａｔｅｄ　ａｔ　ｔｈｅ　ｂｏｔｔｏｍ　ｏｆ　ｔｈｅ　ｕｐｐｅｒ　ｃｒｕｓｔ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎｆａｕｌｔ，ａｎｄ　ａｆｔｅｒｓｈｏｃｋｓ　ｍａｉｎｌｙ　ａｒｅ　ｌｏｃａｔｅｄ　ａｔ　ｔｈｅ　ｄｅｐｔｈ　ｒａｎｇｉｎｇ　ｆｒｏｍ　５ｔｏ　２５ｋｍ　ｉｎ　ｔｈｅ　ｕｐｐｅｒ　ｃｒｕｓｔ．Ｉｔ　ｉｓ　ａｌｓｏ　ｎｏｔｅｄ　ｔｈａｔ　ｅａｒｔｈｑｕａｋｅｓ（Ｍｗ＞３）ａｒｅ　ｄｅｎｓｅｌｙ　ｒｅｃｏｒｄｅｄ　ｉｎ　ｔｈｅ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ　ｓｉｄｅ，ｗｈｉｌｅｒａｒｅｌｙ　ｉｎ　ｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ．Ｔｈｉｓ　ｓｔｕｄｙ　ａｉｍｓ　ｔｏ　ｅｘｐｌａｉｎ　ｔｈｅ　ｍｅｃｈａｎｉｓｍ　ｏｆ　ｔｈｅ　ｅａｒｔｈｑｕａｋｅｇｅｎｅｒａｔｉｏｎ　ｉｎ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ，ａｎｄ　ｔｈｅ　ｅａｒｔｈｑｕａｋｅｓ　ｓｐａｔｉａｌ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｉｎ　ｔｈｅ　ＬｏｎｇｍｅｎＳｈａｎ　ａｒｅａ　ｔｈｒｏｕｇｈ　３－Ｄ　ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌｌｉｎｇ．Ｓｅｖｅｒａｌ　ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｍｏｄｅｌｓ　ｂａｓｅｄ　ｏｎ　ｄｉｆｆｅｒｅｎｔｃｒｕｓｔａｌ　ｒｈｅｏｌｏｇｙ　ａｒｅ　ｓｅｔ　ｕｐ　ｗｉｔｈ　ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎｓ　ｏｆ　ｓｔｅａｄｙ　ｃｏｍｐｒｅｓｓｉｏｎａｌ　ｄｅｆｏｒｍａｔｉｏｎ　ｒａｔｅｃｏｎｓｔｒａｉｎｅｄ　ｂｙ　ＧＰＳ　ｏｂｓｅｒｖａｔｉｏｎｓ．Ｗｅ　ｃａｌｃｕｌａｔｅ　ｔｈｅ　ｌｏｎｇ　ｔｅｒｍ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎ　ｐｒｏｃｅｓｓ　ｉｎ　ｔｈｅｄｉｆｆｅｒｅｎｔ　ｌｉｔｈｏｓｐｈｅｒｅ　ｍｏｄｅｌｓ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ａｒｅａ．Ｏｕｒ　ｒｅｓｕｌｔｓ　ｓｈｏｗ　ｔｈａｔ　ｓｔｒｅｓｓ　ｉｎｃｒｅａｓｅｓａｌｍｏｓｔ　ｌｉｎｅａｒｌｙ　ｗｉｔｈ　ｔｉｍｅ　ｉｎ　ｔｈｅ　ｂｒｉｔｔｌｅ　ｕｐｐｅｒ　ｃｒｕｓｔ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ　ｄｕｒｉｎｇ　ｔｈｅ　ｉｎｔｅｒ－ｓｅｉｓｍｉｃ　ｐｅｒｉｏｄ　ｏｆ　ｓｅｖｅｒａｌ　ｔｈｏｕｓａｎｄ　ｙｅａｒｓ　ｄｕｅ　ｔｏ　ｔｈｅ　ｃｒｕｓｔ　ｓｈｏｒｔｅｎｉｎｇ．Ｓｔｒｅｓｓ　ｃｏｎｃｅｎｔｒａｔｅｓ　ａｔ　ｔｈｅｂｏｔｔｏｍ　ｏｆ　ｔｈｅ　ｂｒｉｔｔｌｅ　ｕｐｐｅｒ　ｃｒｕｓｔ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．Ｔｈｉｓ　ｓｔｒｅｓｓ　ｃｏｎｃｅｎｔｒａｔｉｏｎ　ｐｒｏｃｅｓｓ　ｉｓｒｅｓｐｏｎｓｉｂｌｅ　ｆｏｒ　ｔｈｅ　ｅａｒｔｈｑｕａｋｅ　ｇｅｎｅｒａｔｉｏｎ　ａｔ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．Ｔｈｅ　ｐｒｉｍａｒｙ　ｒｅａｓｏｎ　ｆｏｒ　ｔｈｅｓｔｒｅｓｓ　ｃｏｎｃｅｎｔｒａｔｉｏｎ　ｉｓ　ｔｈｅ　ｌａｒｇｅ　ｖｉｓｃｏｓｉｔｙ　ｄｉｆｆｅｒｅｎｃｅ　ｏｆ　ｔｈｅ　ｍｉｄｄｌｅ　ａｎｄ　ｌｏｗｅｒ　ｃｒｕｓｔｓ　ｂｅｔｗｅｅｎ　ｔｈｅＴｉｂｅｔａｎ　Ｐｌａｔｅａｕ　ａｎｄ　ｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ．Ｓｐａｔｉａｌｌｙ　ｔｈｅ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎ　ｒａｔｅ　ｉｎ　ｔｈｅ　ｕｐｐｅｒ　ｃｒｕｓｔｏｆ　ｔｈｅ　ｅａｓｔｅｒｎ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ　ｉｓ　ｍｕｃｈ　ｈｉｇｈｅｒ　ｔｈａｎ　ｔｈａｔ　ｏｆ　ｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ．Ｔｈｉｓ　ｓｔｒｅｓｓ　ｒａｔｅｓｐａｔｉａｌ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｈｅｌｐｓ　ｔｏ　ｅｘｐｌａｉｎ　ｔｈｅ　ｅａｒｔｈｑｕａｋｅ　ｓｐａｔｉａｌ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｔｈａｔ　ｅａｒｔｈｑｕａｋｅｓ（Ｍｗ＞３）ａｒｅ　ｄｅｎｓｅｌｙ　ｒｅｃｏｒｄｅｄ　ｉｎ　ｔｈｅ　ｅａｓｔｅｒｎ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ，ｗｈｉｌｅ　ｒａｒｅｌｙ　ｉｎ　ｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ．Ｓｔｒｅｓｓｉｎｃｒｅａｓｅｓ　ｅｘｐｏｎｅｎｔｉａｌｌｙ　ｔｏ　ａ　ｓｔｅａｄｙ　ｌｅｖｅｌ　ｉｎ　ｔｈｅ　ｄｕｃｔｉｌｅ　ｍｉｄｄｌｅ　ａｎｄ　ｌｏｗｅｒ　ｃｒｕｓｔｓ，ｗｈｅｒｅ　ｌｅｓｓａｆｔｅｒｓｈｏｃｋｓ　ａｒｅ　ｒｅｃｏｒｄｅｄ　ｄｕｒｉｎｇ　ｔｈｅ　Ｗｅｎｃｈｕａｎ　ｅａｒｔｈｑｕａｋｅ　ａｎｄ　Ｌｕｓｈａｎ　ｅａｒｔｈｑｕａｋｅ．Ｏｕｒ　ｒｅｓｕｌｔｓｓｕｐｐｏｒｔ　ｔｈａｔ　ｔｈｅ　ｒｈｅｏｌｏｇｉｃａｌ　ｓｔｒｕｃｔｕｒｅ　ｏｆ　ｔｈｅ　ｌｉｔｈｏｓｐｈｅｒｅ　ｏｆ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ａｒｅａ　ｉｓ　ａｓ　ｆｏｌｌｏｗｉｎｇ：ｂｒｉｔｔｌｅ　ｕｐｐｅｒ　ｃｒｕｓｔ－ｄｕｃｔｉｌｅ　ｍｉｄｄｌｅ　ａｎｄ　ｌｏｗｅｒ　ｃｒｕｓｔ（ｍｏｒｅ　ｄｕｃｔｉｌｅ　ｉｎ　ｔｈｅ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ　ｔｈａｎ　ｔｈａｔｉｎ　ｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ）－ｄｕｃｔｉｌｅ　ｍａｎｔｌｅ　ｌｉｔｈｏｓｐｈｅｒｅ．Ｋｅｙｗｏｒｄｓ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ，Ｗｅｎｃｈｕａｎ　ｅａｒｔｈｑｕａｋｅ，Ｌｕｓｈａｎ　ｅａｒｔｈｑｕａｋｅ，Ｓｔｒｅｓｓ　ｃｏｎｃｅｎｔｒａｔｉｏｎ，Ｖｉｓｃｏｓｉｔｙ　ｄｉｆｆｅｒｅｎｃｅ，Ｍｏｈｏ　ｓｕｒｆａｃｅ　ｊｕｍｐ１　引言２００８年５月１２日１４时２８分在青藏高原东缘与四川盆地交汇处的龙门山断裂带上发生了汶川Ｍｗ７．９级强烈地震．地震造成的西南向北东方向发展的破裂带长度约３５０ｋｍ，破裂持续时间长达９０ｓ，整个断层面上的平均滑动量约２．４ｍ，最大滑动量达７．３ｍ（张勇等，２００８；张培震等，２００８）；并且沿北川—映秀断裂和彭县—灌县断裂分别形成了长达２４０ｋｍ和７２ｋｍ的地表破裂带，最大垂直位移量约为６．５ｍ，右旋走滑位移量约为４．９ｍ（徐锡伟等，２０１３）．震源机制解显示，汶川地震主震以逆冲为主兼有少量右旋走滑．汶川地震较大余震的“缺失”分析认为汶川地震发生在龙门山断裂带北东段从映秀到青川之间大约３５０ｋｍ的断裂上，而留下了西南段从汶川西南到芦定之间大约１２０ｋｍ的地震亏空区（陈运泰等，２０１３）．５年后，２０１３年４月２０日８时２分在汶川地震的亏空区发生了芦山地震．地震破裂过程研究表明，断层滑动面长宽约为３０ｋｍ×３０ｋｍ，最大滑动量为１．６ｍ；破裂起始点接近于地震滑动量集中区，破裂面从震源处向下延伸至２０ｋｍ左右深的地方，破裂面并未到达地表，破裂过程没有明显的方向性；震源机制解显示，芦山地震主要以逆冲为主兼有非常少量的右旋走滑分量（陈运泰等，２０１３）．结合地震勘探成果所揭示的龙门山断裂带地壳结构，余震震源定位工作表明龙门山区域地震空间分布特征为：纵向上，龙门山断裂带这两次地震主震均发生在上地壳的底部，汶川地震的震源深度为１９ｋｍ（ＵＳＧＳ），芦山地震的震源深度为５０４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５７卷　１４ｋｍ（ＵＳＧＳ），全球其他各个研究机构对芦山地震的震源深度定位结果在１２～１９ｋｍ深度范围之间（杜方等，２０１３）；两次地震的绝大部分余震均发生在上地壳范围（５～２５ｋｍ）（张瑞青等，２００８；黄媛等，２００８；刘巧霞等，２００８；陈九辉等，２００９；张广伟等，２０１３；吕坚等，２０１３；陈运泰等，２０１３），而在中、下地壳深度范围内鲜见余震发生；横向上，地震（Ｍｗ＞３）在龙门山断裂带青藏高原一侧密集分布，且有历史大震发生，如：１６３０年Ｍ６．５级虎牙地震、１９１３年Ｍ７．０级叠溪地震、１９３２年Ｍ７．５级叠溪地震、１９６０年Ｍ６．７级松潘地震与１９７６年３次震级为６．６＜Ｍ＜７．３级的松潘地震群；而四川盆地地震（Ｍｗ＞３）稀少（李勇等，２００６；滕吉文等，２００８；Ｒｏｂｅｒｔ　ｅｔ　ａｌ．，２０１０；雷兴林等，２０１３）．在低地形变速率（约３ｍｍ／ａ）的龙门山断裂带相隔５年发生两次强震：为探讨龙门山断裂带地震孕震成因，不同学者就青藏高原东缘的构造运动模式，提出了两种不同的概念性地震地质孕震模型．Ｔａｐｐｏｎｎｉｅｒ等（２００１）主张龙门山断裂带可能贯穿青藏高原东缘整个脆性地壳；Ｈｕｂｂａｒｄ等（２００９）认为汶川地震是龙门山断裂带在构造挤压环境下的再一次活动结果，并强调地壳缩短是地震发生的首要机制．更多的研究（张培震等，２００８；滕吉文等，２００８；Ｒｏｙｄｅｎ　ｅｔ　ａｌ．，２００８；Ｂｕｒｃｈｆｉｅｌ　ｅｔ　ａｌ．，２００８；Ｓｔｏｎｅｅｔ　ａｌ．，２００９；Ｒｏｂｅｒｔ　ｅｔ　ａｌ．，２０１０；柳畅等，２０１２ｂ）则支持如下的观点：认为印度板块对欧亚板块的推挤作用造成了青藏高原物质的东向运动，高原柔性的中、下地壳物质在龙门山断裂带处遭到相对坚硬的四川盆地的阻挡之后，部分中、下地壳物质在龙门山断裂带下堆积、产生应力集中而导致龙门山断裂带地震的发生．这一观点所描述的青藏高原脆性上地壳与韧性的中、下地壳组合的地壳分层流变结构，区别于Ｔａｐｐｏｎｎｉｅｒ等（２００１）所认为的全脆性地壳流变结构．总结以上的观点，其共同点在于认同青藏高原的东向挤压作用，但又强调了完全不同的地壳流变结构．如何对这两种概念性地震地质孕震模型从数值模拟的角度作出验证，进一步对龙门山断裂带地震发生机理加以更深层次的认识，是一个尚未解决且需要探讨的问题．为此，本文将建立四种不同流变结构的龙门山断裂带三维岩石圈模型，数值模拟不同岩石圈模型中的应力积累过程与分布；并探讨青藏高原中、下地壳不同的黏滞系数，以及青藏高原下地壳流Ｃｈａｎｎｅｌ　ｆｌｏｗ（Ｒｏｙｄｅｎ　ｅｔ　ａｌ．，１９９７；Ｃｌａｒｋｅｔ　ａｌ．，２０００，２００５；Ｂｅａｕｍｏｎｔ　ｅｔ　ａｌ．，２００１，２００４；石耀霖等，２００８）对应力积累过程的影响．２　地质构造背景龙门山断裂带处于年轻活跃的青藏高原与古老稳定的扬子克拉通的交汇地带，为ＮＥ－ＳＷ走向上长约４７０ｋｍ、宽约５０ｋｍ的活动断裂带．该断裂带在空间上呈ＮＥ－ＮＮＥ向展布，且以北ＮＷ－ＳＥ方向逆冲为主而兼具少量右旋走滑分量；是由汶川—茂汶逆断裂、映秀—北川逆断裂、彭县—灌县断裂、龙门山山前断裂、山前隐伏断裂段和相应的推覆体组成的一组断裂系．２００８年汶川地震发生在映秀—北川断裂带上，２０１３年芦山地震发生在山前断裂南段大川—双石断裂上．为揭示龙门山断裂带的岩石圈结构，该地域进行过大量的地震探测工作：如层析成像工作（Ｈｕａｎｇｅｔ　ａｌ．，２００７；郭飚等，２００９），接收函数工作（刘启元等，２００９；Ｒｏｂｅｒｔ　ｅｔ　ａｌ．，２０１０；郑勇等，２０１３），人工地震勘探工作（Ｗａｎｇ　ｅｔ　ａｌ．，２００５；Ｗａｎｇ　ｅｔ　ａｌ．，２０１０）和一些地质勘查工作（Ｂｕｒｃｈｆｉｅｌ，１９９５，２００８）．以上研究结果表明从青藏高原东缘到四川盆地地壳结构变化强烈，青藏高原一侧地壳厚度在６３ｋｍ左右，而四川盆地一侧地壳厚度在４５ｋｍ左右．在龙门山区域横向５０ｋｍ范围内地壳厚度的垂直变化幅度可达１８ｋｍ左右，地形高差在龙门山陡降近５ｋｍ；因此，结合Ｍｏｈｏ面形状与地表地形看，龙门山区域地壳结构呈瓶颈状（图２）．大陆岩石圈一般存在脆性的上地壳、柔性的下地壳和较强上地幔这种三明治式的分层流变结构，尤其在青藏高原这种地壳厚、地热流量密度高的地区会更为显著（Ｒｏｙｄｅｎ　ｅｔ　ａｌ．，１９９７）．岩石流变性质与温度密切相关，安美建利用地震波速估算了上地幔５０～２００ｋｍ的温度（安美建等，２００７）；石耀霖据此温度模型给出了中国大陆岩石圈不同层位的等效黏滞系数（石耀霖等，２００８），得到了相应区域的三维流变结构．胡圣标研究了川滇地区相对偏高的平均大地热流值（胡圣标等，２００１），认为川滇地区的中、下地壳较热，介质强度可能相对较软．ＧＰＳ观测结果也显示，川滇地区的现今运动模式支持较硬的脆性上地壳和软弱的柔性中、下地壳的分层流变结构（Ｓｈｅｎ　ｅｔ　ａｌ．，２００５）．Ｐ波层析成像工作所揭示的四川盆地下方地壳速度结构，认为扬子克拉通的稳定速度结构可能延伸至２５０ｋｍ深度（Ｌｉ　ｅｔ　ａｌ．，２００６，２００８）．并且，中６０４　２期柳畅等：地壳流变结构控制作用下的龙门山断裂带地震发生机理生代与新生代的构造活动对四川盆地周边构造单元产生了显著的变形作用，而四川盆地并未受到明显的影响（Ｌｉ　ｅｔ　ａｌ．，２００６，２００８），从而表明四川盆地下方岩石圈的力学强度较周边构造单元较强．３　三维有限元模型３．１　模型参数本文建立如图１中黑色虚线方框所示区域的三维黏弹性Ｍａｘｗｅｌｌ体岩石圈模型（见图２）．坐标系Ｘ轴（ＳＥ４１°）垂直于龙门山断裂带走向，长５００ｋｍ；Ｙ轴（ＮＥ４９°）平行于龙门山断裂带走向，长１００ｋｍ，图１　青藏高原东缘与四川盆地区域构造与地表ＧＰＳ水平速度场红色圆点代表２００８汶川地震Ｍｗ５．０级以上余震，粉色圆点代表２０１３芦山地震Ｍｗ５．０级以上的余震．黄色箭头代表地表ＧＰＳ水平速度（Ｓｈｅｎ　ｅｔ　ａｌ．，２００５）．黑色小三角代表地名，分别为：１，汶川；２，成都；３，叠溪；４，松潘；５，卢霍；６，康定；７，芦山．黑色虚线框表示数值模拟的模型区域．Ｄ－Ｓ　ｆａｕｌｔ代表大川—双石断裂．Ｆｉｇ．１　Ｍａｐ　ｏｆ　ｔｈｅ　ｅａｓｔｅｒｎ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ　ａｎｄ　ｔｈｅＳｉｃｈｕａｎ　Ｂａｓｉｎ　ｓｈｏｗｉｎｇ　ｔｈｅ　Ｌｏｎｇｍｅｎｓｈａｎ　ｆａｕｌｔ　ａｎｄ　ｔｈｅｈｏｒｉｚｏｎｔａｌ　ＧＰＳ　ｖｅｌｏｃｉｔｙＴｈｅ　ｅｐｉｃｅｎｔｒｅｓ　ｏｆ　ｔｈｅ　ｍａｉｎ　ｓｈｏｃｋ　ｏｆ　２００８Ｗｅｎｃｈｕａｎ　ｅａｒｔｈｑｕａｋｅａｎｄ　２０１３Ｌｕｓｈａｎ　ｅａｒｔｈｑｕａｋｅ，ａｎｄ　ａｆｔｅｒｓｈｏｃｋｓ（Ｍｗ＞５．０）ｏｆｔｈｅ　Ｗｅｎｃｈｕａｎ　ｅａｒｔｈｑｕａｋｅ（ｒｅｄ　ｄｏｔｓ）ａｎｄ　Ｌｕｓｈａｎ　ｅａｒｔｈｑｕａｋｅ（ｐｉｎｋ　ｄｏｔｓ）．Ｔｈｅ　ｙｅｌｌｏｗ　ａｒｒｏｗ　ｉｎｄｉｃａｔｅｓ　ｔｈｅ　ＧＰＳ　ｖｅｌｏｃｉｔｙ（Ｓｈｅｎｅｔ　ａｌ．，２００５）．Ｔｈｅ　ｂｌａｃｋ　ｔｒｉａｎｇｌｅｓ　ｉｎｄｉｃａｔｅ　ｔｈｅ　ｌｏｃａｔｉｏｎｓ　ｏｆｃｉｔｉｅｓ．Ｔｈｅ　ｎｕｍｂｅｒｓ　ｍａｒｋ　ｔｈｅ　ｃｉｔｙ　ｎａｍｅｓ：１，Ｗｅｎｃｈｕａｎ；２，Ｃｈｅｎｇｄｕ；３，Ｄｉｅｘｉ；４，Ｓｏｎｇｐａｎ；５，Ｌｕｈｕｏ；６，Ｋａｎｇｄｉｎｇ；７，Ｌｕｓｈａｎ．Ｔｈｅ　ｂｌａｃｋ　ｓｏｌｉｄ　ｌｉｎｅｓ　ｉｎｄｉｃａｔｅ　ｔｈｅ　ｆａｕｌｔｓ．Ｄ－Ｓ　ｆａｕｌｔｉｎｄｉｃａｔｅｓ　ｔｈｅ　Ｄａｃｈｕａｎ－Ｓｈｕａｎｇｓｈｉ　ｆａｕｌｔ．Ｔｈｅ　ｂｌａｃｋ　ｄａｓｈｅｄ　ｌｉｎｅｓｆｏｒｍｉｎｇ　ａ　ｒｅｃｔａｎｇｌｅ　ｍａｒｋ　ｔｈｅ　ａｒｅａ　ｗｈｅｒｅ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎ　ｉｓｓｉｍｕｌａｔｅｄ．包含了龙门山断裂带在２００８汶川地震未发生破裂的西南段部分；Ｚ轴向上（岩石圈的深度在Ｚ轴负向），模型中深度为从地面至地下１００ｋｍ；坐标原点位于模型上表面的西南端．模型具体考虑龙门山断裂带青藏高原与四川盆地的地形高差、地壳厚度在龙门山断裂带下方的突变与中、下地壳的黏滞系数在高原与盆地下方的差异．本文计算所涉及的黏弹性介质的３个物质参数分别为：杨氏模量Ｅ，泊松比ν和黏滞系数η．岩石杨氏模量Ｅ与泊松比ν依据波速计算．一般随深度增加，地震波速增加，岩石的杨氏模量亦随深度增加．依据所选择物质参数，龙门山断裂带岩石圈１００ｋｍ深度范围内介质杨氏模量随深度的变化曲线见图３．岩石圈模型中青藏高原与四川盆地的上地壳及岩石圈上地幔部分具有相同的流变性质与参数，而中、下地壳流变性质有非常明显的差异．对于黏滞系数而言，青藏高原中、下地壳较软；而四川盆地下方中、下地壳相对坚硬．参照石耀霖等（２００８）的中国大陆岩石圈等效黏滞系数的研究结果，相应的层位黏滞系数取为：高原下方中地壳为３．０×１０２０～５．０×１０２０　Ｐａ·ｓ，下地壳为３．０×１０１９～６．０×１０１９Ｐａ·ｓ，四川盆地中地壳为１．０×１０２２～１．２×１０２２Ｐａ·ｓ，下地壳为７．０×１０２１～８．０×１０２１Ｐａ·ｓ．相应岩石圈各层位黏弹性物质的松驰时间（介质黏滞系数与杨氏模量之比的二分之一）分别为：青藏高原与四川盆地上地壳为３００００ａ左右；青藏高原中、下地壳分别为２００ａ、５０ａ左右，四川盆地中、下地壳分别为１０００ａ、５００ａ左右；青藏高原与四川盆地岩石圈上地幔为３００ａ左右．模型岩石圈物质黏滞系数随深度分布的变化曲线见图４．针对以上黏弹性Ｍａｘｗｅｌｌ体模型，在数千年的时间尺度内考察岩石圈的应力积累过程中，不论是青藏高原还是四川盆地，在地壳的缩短过程中，上地壳将表现为脆性，而中、下地壳介质的流变性质则可以充分得到体现（这取决于黏弹性介质的松驰时间）．由于Ｍａｘｗｅｌｌ体具有这种可以自然处理脆－韧性转变的优点，本文有限元模型实际可以有效地描述弹性上地壳覆盖在黏弹性中、下地壳和岩石圈上地幔的应力积累过程．我们将该岩石圈模型命名为参考模型０，以方便与下文中另外的几种岩石圈模型加以区别．模型０从上到下一共分为１３层，其中１—７层为上地层，８—９层为中地壳，１０—１１层为下地壳，１２—１３层为岩石圈上地幔，图２中不颜色代表不同分层的物质７０４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５７卷　图２　龙门山断裂带三维岩石圈模型模型尺寸为５００ｋｍ×１００ｋｍ×１００ｋｍ（Ｘ×Ｙ×Ｚ）．Ｙ轴方向平行于龙门山断裂带走向．不同颜色代表不同岩石圈层位的不同物质，依据对应的黏滞系数，青藏高原中、下地壳物质较软而四川盆地中、下地壳物质相对较硬．黑色星星代表２０１３芦山Ｍｗ６．６级地震的震中位置．Ｆｉｇ．２　Ｌｉｔｈｏｓｐｈｅｒｉｃ　ｍｏｄｅｌ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　ＳｈａｎＴｈｅ　ｓｉｚｅ　ｏｆ　ｔｈｉｓ　ｍｏｄｅｌ　ｉｓ　５００ｋｍ×１００ｋｍ×１００ｋｍ（Ｘ×Ｙ×Ｚ）．Ｔｈｅ　Ｙａｘｉｓ　ｉｓ　ｐａｒａｌｌｅｌ　ｔｏ　ｔｈｅ　ｓｔｒｉｋｅ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．Ｔｈｅｄｉｆｆｅｒｅｎｔ　ｃｏｌｏｕｒｓ　ｒｅｆｅｒ　ｔｏ　ｔｈｅ　ｄｉｆｆｅｒｅｎｔ　ｍａｔｅｒｉａｌ　ｐａｒａｍｅｔｅｒｓ　ｏｆ　ｌａｙｅｒｓ　ｉｎ　ｔｈｅ　ｌｉｔｈｏｓｐｈｅｒｅ．Ｔｈｅ　ｍａｔｅｒｉａｌｓ　ｏｆ　ｔｈｅ　ｍｉｄｄｌｅ　ａｎｄ　ｌｏｗｅｒ　ｃｒｕｓｔｓ　ｏｆｔｈｅ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ　ａｒｅ　ｓｏｆｔｅｒ　ｔｈａｎ　ｔｈｏｓｅ　ｏｆ　ｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ，ｃｏｒｒｅｓｐｏｎｄｉｎｇ　ｔｏ　ｔｈｅ　ｌｅｖｅｌ　ｏｆ　ｖｉｓｃｏｓｉｔｙ．Ｂｌａｃｋ　ｓｔａｒ　ｉｎｄｉｃａｔｅｓ　ｔｈｅｅｐｉｃｅｎｔｒｅ　ｏｆ　ｔｈｅ　２０１３Ｌｕｓｈａｎ　ｅａｒｔｈｑｕａｋｅ．图３　龙门山断裂带岩石圈物质杨氏模量随深度分布图（Ｗａｎｇ　ｅｔ　ａｌ．，２００５；Ｈｕａｎｇ　ｅｔ　ａｌ．，２００７；石耀霖等，２００８；Ｒｏｂｅｒｔ　ｅｔ　ａｌ．，２００９；Ｗａｎｇ　ｅｔ　ａｌ．，２０１０；Ｚｈａｎｇ　ｅｔ　ａｌ．，２０１１）Ｆｉｇ．３　Ｒｅｌａｔｉｏｎｓｈｉｐ　ｏｆ　Ｙｏｕｎｇ′ｓ　ｍｏｄｕｌｕｓ　ｗｉｔｈ　ｄｅｐｔｈｉｎ　ｔｈｅ　ｌｉｔｈｏｓｐｈｅｒｅ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ａｒｅａ图４　龙门山断裂带岩石圈各层位岩石介质黏滞系数在青藏高原一侧与四川盆地一侧随深度分布图（石耀霖等，２００８）Ｆｉｇ．４　Ｒｅｌａｔｉｏｎｓｈｉｐ　ｏｆ　ｖｉｓｃｏｓｉｔｙ　ｗｉｔｈ　ｄｅｐｔｈ　ｉｎｔｈｅ　ｌｉｔｈｏｓｐｈｅｒｅ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ａｒｅａ参数，其具体介质参数见表１与表２．用六面体单元对模型进行网格划分，单元总数为１８２４００个，节点总数为２７９６００个．３．２　边界条件和初始条件大量的研究对青藏高原的流变结构进行过探讨（Ｒｏｙｄｅｎ　ｅｔ　ａｌ．，１９９７，２００８；Ｃｌａｒｋ　ｅｔ　ａｌ．，２０００，２００５；Ｂｅａｕｍｏｎｔ　ｅｔ　ａｌ．，２００１，２００４），并对青藏高原下地壳的运动模式持有不同的见解（Ｒｏｙｄｅｎ　ｅｔ　ａｌ．，２００８）．Ｃｌａｒｋ等（２０００，２００５）认为青藏高原重力驱动作用下的低黏度（黏滞系数为２．０×１０１８Ｐａ·ｓ）的下地壳流动速度大约在８０ｍｍ／ａ左右（为地表运动速度的８倍左右）．Ｃａｏ等（２００９）通过建立青藏高原的三维岩石圈模型拟合了青藏高原地表的水平运动速度，其结果表明青藏高原下地壳的流动速度仅比地表运动速度快８ｍｍ／ａ左右，并且所估计的下地壳黏滞系数也比Ｒｏｙｄｅｎ与Ｃｌａｒｋ的估计结果高出一个量级，为３．０×１０１９Ｐａ·ｓ．而Ｂｅｎｄｉｃｋ等（２００７）和Ｗａｎｇ等（２００８）则支持下地壳与上、中地壳和岩石圈上地幔之间无差异运动的观点．因此，我们在有限元模型中对青藏高原一侧的物质运动采用了三种８０４　２期柳畅等：地壳流变结构控制作用下的龙门山断裂带地震发生机理表１　青藏高原东缘岩石圈物质参数Ｔａｂｌｅ　１　Ｍａｔｅｒｉａｌ　ｐａｒａｍｅｔｅｒｓ　ｏｆ　ｒｏｃｋｓ　ｉｎ　ｔｈｅ　ｌｉｔｈｏｓｐｈｅｒｅｏｆ　ｔｈｅ　ｅａｓｔｅｒｎ　ｍａｒｇｉｎ　ｏｆ　ｔｈｅ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ层Ｈ／ｋｍ　Ｅ／１０１０　Ｐａνρ／（ｋｇ·ｍ－３）η／（１０２３Ｐａ·ｓ）上地壳３～－２　３．７５　０．２１　２５００　１．０－２～－５　５．１５　０．２６　２５６９　１．０－５～－１０　６．０　０．２４　２６８５　１．０－１０～－１３　８．１　０．２５　２６８４　１．０－１３～－１６　８．６　０．２４　２６８４　１．０－１６～－１９　９．１　０．２６　２６８３　１．０－１９～－２２　９．７　０．２６　２６８３　１．０中地壳－２２～－３２　１１　０．２４　２７５４　０．００５－３２～－４２　１２　０．２６　２７７３　０．００３下地壳－４２～－５１　１３　０．２５　２８３５　０．０００６－５１～－６０　１４　０．２６　２９５９　０．０００３岩石圈地幔－６０～－８０　１７　０．２５　３２７０　０．００６－８０～－１００　１８　０．２８　３２７０　０．００４ 注：Ｈ为深度，Ｅ为岩石的杨氏模量，ν为泊松比，ρ为密度，η为黏滞系数（Ｗａｎｇ　ｅｔ　ａｌ．，２００５；Ｈｕａｎｇ　ｅｔ　ａｌ．，２００７；石耀霖等，２００８；Ｒｏｂｅｒｔ　ｅｔ　ａｌ．，２００９；Ｗａｎｇ　ｅｔ　ａｌ．，２０１０；Ｚｈａｎｇ　ｅｔ　ａｌ．，２０１１）．表２　四川盆地岩石圈物质参数Ｔａｂｌｅ　２　Ｍａｔｅｒｉａｌ　ｐａｒａｍｅｔｅｒｓ　ｏｆ　ｒｏｃｋｓ　ｉｎ　ｔｈｅｌｉｔｈｏｓｐｈｅｒｅ　ｏｆ　ｔｈｅ　Ｙａｎｇｔｚｅ　ｃｒａｔｏｎ层Ｈ／ｋｍ　Ｅ／１０１０Ｐａνρ／（ｋｇ·ｍ－３）η／（１０２３Ｐａ·ｓ）上地壳０～－３＊＊＊＊－３～－５　３．７５　０．２１　２５００　１．０－５～－７　５．１５　０．２６　２５６９　１．０－７～－９　６．０　０．２４　２６８５　１．０－９～－１１　８．１　０．２５　２６８４　１．０－１１～－１３　８．６　０．２４　２６８４　１．０－１３～－１５　９．１　０．２６　２６８３　１．０－１５～－１９　９．７　０．２６　２６８３　１．０中地壳－１９～－２６　１１　０．２４　２７５４　０．１２－２６～－３２　１２　０．２６　２７７３　０．１下地壳－３２～－３７　１３　０．２５　２８３５　０．０８－３７～－４５　１４　０．２６　２９５９　０．０７岩石圈地幔－４５～－８０　１７　０．２５　３２７０　０．００６－８０～－１００　１８　０．２８　３２７０　０．００４ 注：＊表示该深度段没有物质．不同的位移边界条件，以满足以上几种不同的观点，见图５．边界条件１（ＢＣ１）：地壳运动观测网络给出了研究地区１９９８—２００４年的地表ＧＰＳ速度场结果，见图１．将ＧＰＳ在研究区域边界附近实测速度值插值到四个侧边界，作为水平速度约束条件，且假定从地表到１００ｋｍ深度保持一致，以拟合藏青高原下地壳的无差异运动的情况（Ｂｅｎｄｉｃｋ　ｅｔ　ａｌ．，２００７；Ｗａｎｇｅｔ　ａｌ．，２００８）；垂直方向位移可以保持自由．上表面为自由边界，即法向应力和剪应力均为零．对于底部边界，鉴于目前对该区域岩石圈上地幔运动状态的未知，暂且将底面垂直方向速度约束为０，而水平方向自由．在黏弹性问题中，边界条件随时间的变化也是重要的问题，在目前的模拟中，我们初步假定边界位移速率不随时间变化，见图５．边界条件２（ＢＣ２）：我们采用Ｃａｏ等（２００９）的观点，认为青藏高原下地壳以Ｃｈａｎｎｅｌ　ｆｌｏｗ形式的运动速度Ｖ１比地表运动速度Ｖ快８ｍｍ／ａ．其他边界条件同边界条件１，见图５．边界条件３（ＢＣ３）：青藏高原下地壳以Ｃｈａｎｎｅｌｆｌｏｗ形式的运动速度Ｖ２为地表运动速度Ｖ的３倍．需要指出的是，Ｖ２的取值并没有理论依据，目的是为了测试下地壳不同的流动速度对模型构造应力积累的影响．其他边界条件同边界条件１，见图５．初始条件是构造应力场模拟中最困难的问题，尽管现今在龙门山断裂带进行了不同程度的应力测量结果，但基本都限于沉积层深度的钻孔应力测量，而对深部三维应力分布和应力演变历史几乎仍然没有定量的资料．在这种情况下，我们只能先假定初始应力为０，然后计算在定长的边界位移速率条件下应力的演变．虽然我们不可能模拟现今真实的应力状态，但我们可以了解在定长边界位移速率下的应力增长率．应力增长率高的地方，未必一定是目前应力最高的地方，但如果初始应力类似，则较高应力增长率的地方则更有可能是现今应力绝对值较大的区域．本文主要就应力增长率与地震活动性的关系进行讨论．４　计算结果计算所使用的程序是利用“飞箭有限元程序自动生成系统（ＦＥＰＧ）”生成的Ｍａｘｗｅｌｌ体三维有限元计算程序，程序计算的可靠性已经在大量事例中通过验证．在保证程序可靠性的前提下，对模型进行计算是合理的．本文采用１ａ为一个时间步长，根据时间步长逐步加载边界位移约束．本文计算主要的目标是计算龙门山断裂带岩石圈各层位在数千年时间尺度以上的应力积累过程，在计算过程中只考虑模型在边界条件作用下产生的构造应力的积累变化，并不考虑重力的因素．９０４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５７卷　图５　模型边界条件剖面１垂直于龙门山断裂带走向，剖面２平行于龙门山断裂带走向．ＢＣ１，ＢＣ２和ＢＣ３分别为边界条件１，２和３的缩写．Ｆｉｇ．５　Ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎｓＴｒａｎｓｅｃｔｉｏｎ　１ｉｓ　ｐｅｒｐｅｎｄｉｃｕｌａｒ　ｔｏ　ｔｈｅ　ｓｔｒｉｋｅ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ，ｗｈｉｌｅ　ｔｒａｎｓｅｃｔｉｏｎ　２ｉｓ　ｐａｒａｌｌｅｌ　ｔｏ　ｔｈｅ　ｓｔｒｉｋｅ　ｏｆｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．ＢＣ１，ＢＣ２ａｎｄ　ＢＣ３ａｒｅ　ａｂｂｒｅｖｉａｔｉｏｎｓ　ｏｆ　ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎｓ　１，２ａｎｄ　３，ｒｅｓｐｅｃｔｉｖｅｌｙ． 黏弹性Ｍａｘｗｅｌｌ体在外部载荷下的变形，不仅与边界条件及其随时间的变化有关，而且与初始条件和以前的应力演变历史有关．鉴于初始条件缺乏测量资料，只能假定为零应力、初始应变速率为０．边界条件也假定了位移速率为常量、不随时间变化．在这种条件下，开始的数百年内，下地壳、岩石圈上地幔等黏滞系数较小、弛豫时间较短（数十到数百年）的柔性可以占主导的层位，在压缩位移下的应力增长与柔性介质内的应力松弛将达到平衡，应力维持在一个较低水平而不再增加．相反，上地壳黏滞系数高（弛豫时间达数万年）的层位，在几百到几万年的期间内，弹性仍然占主导，应力随压缩几乎可以接近线性的速率增长．该结论在我们的另一研究中已经得到证明（柳畅等，２０１２ｂ）．因此，在模型中各层位应力积累达到稳定速率增长后，我们分别取经过芦山地震震源的纵剖面（图５中剖面１）和芦山地震震源深度的横切面上的应力积累速率分布，以分析龙门山断裂带数千年以上的应力积累及其与区域地震活动性之间的关系．４．１　应力积累结果在边界条件１（下地壳无差异运动）作用下参考模型０，垂直于龙门断裂带方向的压应力积累分布情况如下．我们取经过芦山地震震中、且垂直于龙门山断裂带走向的纵剖面（图５中剖面１）上的应力积累速率分布，见图６ａ，以揭示应力积累与龙门山区域地震纵向空间分布特征之间的关系．结果显示在脆性的上地壳应力以近乎定值的速率持续增长，在龙门山断裂带上地壳底部出现应力集中现象，最大应力积累速率为－３．８６ｋＰａ／ａ．而中、下地壳及岩石圈上地幔的应力积累速率因物质的黏性应力松驰而为０．这种应力积累状态表明，在构造挤压作用下龙门山断裂带的上地壳内应力可以持续增长；在应力积累速率最大的龙门山断裂带上地壳底部，应力可以优先增长至该部分岩石破裂强度而导致主震的发生（汶川地震主震震源深度１９ｋｍ，芦山地震主震震源深度１４ｋｍ），进而触发上地壳内部的余震（５～２５ｋｍ）；而在中、下地壳内，当构造挤压作用下的应力增长与介质黏性应力松驰效应达到动态平衡后，其应力趋于稳定，不再增长，因此该层位应力积累很难达到岩石的破裂强度．汶川地震和芦山地震的余震震源定位结果均表明，在龙门山断裂带中、下地壳０１４　２期柳畅等：地壳流变结构控制作用下的龙门山断裂带地震发生机理图６　（ａ）边界条件１（ＢＣ１）作用下模型０中经过芦山地震震中且垂直于龙门山断裂带的剖面上的应力积累速率分布状况；（ｂ）边界条件２（ＢＣ２）作用下该剖面上的应力积累速率分布状况．应力在龙门山断裂带下集中积累，在龙门山断裂带上地壳底部有最大增长率．黑色星星表示汶川地震和芦山地震主震震源在该剖面上的投影ＴＰ表示青藏高原，ＬＭＳ表示龙门山，ＳＢ表示四川盆地Ｆｉｇ．６　（ａ）Ｎｏｒｍａｌ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎ　ｒａｔｅ　ｄｉｓｔｒｉｂｕｔｉｏｎｉｎ　ｔｈｅ　ｃｒｏｓｓ　ｓｅｃｔｉｏｎ　ｏｆ　ｍｏｄｅｌ　０ｉｎ　ｔｈｅ　ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎ１（ＢＣ１）；（ｂ）Ｎｏｒｍａｌ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎ　ｒａｔｅ　ｄｉｓｔｒｉｂｕｔｉｏｎｉｎ　ｔｈｅ　ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎ　２（ＢＣ２）；Ｔｈｅ　ｃｒｏｓｓ　ｓｅｃｔｉｏｎ　ｉｓｐｅｒｐｅｎｄｉｃｕｌａｒ　ｔｏ　ｔｈｅ　ｓｔｒｉｋｅ　ｏｆ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．Ｓｔｒｅｓｓ　ｃｏｎｃｅｎｔｒａｔｅｓ　ａｔ　ｔｈｅ　ｂｏｔｔｏｍ　ｏｆ　ｔｈｅ　ｕｐｐｅｒ　ｃｒｕｓｔ　ｏｆｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．Ｔｈｅ　ｌａｒｇｅｓｔ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎｒａｔｅ　ｉｓ　ｌｏｃａｔｅｄ　ａｔ　ｔｈｅ　ｂｏｔｔｏｍ　ｏｆ　ｔｈｅ　ｕｐｐｅｒ　ｃｒｕｓｔ　ｏｆ　ｔｈｅＬｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．Ｔｈｅ　ｂｌａｃｋ　ｓｔａｒｓ　ｉｎｄｉｃａｔｅ　ｔｈｅｈｙｐｏｃｅｎｔｅｒｓ　ｏｆ　ｔｈｅ　２００８Ｗｅｎｃｈｕａｎ　ｅａｒｔｈｑｕａｋｅ　ａｎｄ　ｔｈｅ２０１３Ｌｕｓｈａｎ　ｅａｒｔｈｑｕａｋｅ．ＴＰ　ｉｎｄｉｃａｔｅｓ　ｔｈｅ　ＴｉｂｅｔａｎＰｌａｔｅａｕ，ＬＭＳ　ｉｎｄｉｃａｔｅｓ　ｔｈｅ　Ｌｏｎｇｍｅｎ　Ｓｈａｎ，ａｎｄ　ＳＢｉｎｄｉｃａｔｅｓ　ｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ深度范围内鲜见余震发生．我们取上地壳深度１４ｋｍ处（芦山地震源深度）横切面上的应力积累速率分布，见图７．结果显示青藏高原上地壳应力增长率（约－２．５ｋＰａ／ａ）远远高于四川盆地的（约－１．０ｋＰａ／ａ），而最大应力增长率（－３．４ｋＰａ／ａ）则位于龙门山断裂带．这种应力积累分布现象有利于解释龙门山区域地震横向空间分布特征：龙门山断裂带青藏高原一侧地震（Ｍｗ＞３）密集分布（李勇等，２００６；滕吉文等，２００８；雷兴林图７　边界条件１作用下模型０中深度１４ｋｍ处横切面上的应力积累速率分布状况；青藏高原应力增长速率远远高于四川盆地，在龙门山断裂带集中积累有最大值．黑色星星表示Ｍｗ６．６级芦山地震主震．黑色三角标志城市位置．白色线表示断裂Ｆｉｇ．７　Ｎｏｒｍａｌ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎ　ｒａｔｅ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｉｎｔｈｅ　ｍａｐ　ａｔ　ｔｈｅ　ｄｅｐｔｈ　ｏｆ　１４ｋｍ　ｉｎ　ｍｏｄｅｌ　０ｉｎ　ｔｈｅｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎ　１（ＢＣ１）．Ｔｈｅ　ｓｔｒｅｓｓ　ａｃｃｕｍｕｌａｔｉｏｎｒａｔｅ　ｉｎ　ｔｈｅ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ　ｉｓ　ｍｕｃｈ　ｈｉｇｈｅｒ　ｔｈａｎ　ｔｈａｔ　ｏｆｔｈｅ　Ｓｉｃｈｕａｎ　Ｂａｓｉｎ．Ｔｈｅ　ｓｔｒｅｓｓ　ｃｏｎｃｅｎｔｒａｔｅｓ　ａｔ　ｔｈｅＬｏｎｇｍｅｎ　Ｓｈａｎ　ｆａｕｌｔ．Ｔｈｅ　ｂｌａｃｋ　ｓｔａｒ　ｉｎｄｉｃａｔｅｓ　ｔｈｅｈｙｐｏｃｅｎｔｒｅ　ｏｆ　ｔｈｅ　２０１３Ｌｕｓｈａｎ　ｅａｒｔｈｑｕａｋｅ．Ｂｌａｃｋｔｒｉａｎｇｌｅｓ　ｉｎｄｉｃａｔｅ　ｔｈｅ　ｌｏｃａｔｉｏｎ　ｏｆ　ｃｉｔｉｅｓ等，２０１３），且有历史大震发生，如：１６３０年Ｍ６．５级虎牙地震、１９１３年Ｍ７．０级叠溪地震、１９３２年Ｍ７．５级叠溪地震、１９６０年Ｍ６．７级松潘地震与１９７６年３次震级为６．６＜Ｍ＜７．３级的松潘地震群；最大级别地震———汶川Ｍｗ７．９级地震，发生在龙门山断裂带；而四川盆地地震稀少（Ｍｗ＞３）（李勇等，２００６；滕吉文等，２００８；雷兴林等，２０１３）．４．２　青藏高原下地壳流的影响我们计算了边界条件２（下地壳流动速度快于地表８ｍｍ／ａ）与边界条件３（下地壳流速度为地表速度的３倍）作用下模型０中的应力积累速率分布状态，以探讨青藏高原下地壳流对龙门山断裂带应力积累过程的影响．计算结果表明，除了在下地壳的上、下底界面处（相应于图６ｂ中Ｌ层位的上、下界面）因下地壳的快速流动而引起的拖曳力（张应力）外，边界条件２和３作用下模型的应力积累速率分布与边界条件１作用下的结果并无太大区别．因此，我们只展示边界条件２作用下，模型０中的应力积累速率分布状况（图６ｂ）．结果显示龙门山断裂带上地壳底部同样出现应力集中现象，且该部位应力增１１４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５７卷　长速率有所增大．取芦山地震震源处的应力增长率分析，边界条件１，边界条件２和边界条件３作用下该处的应力增长率分别为－３．６０ｋＰａ／ａ，－３．７５ｋＰａ／ａ和－３．９ｋＰａ／ａ，见表５．可见，下地壳流动有助于龙门山断裂带的应力积累，并且随流动速度增大而增大，但增大量并不显著．４．３　地壳流变结构的影响为探讨青藏高原地壳流变结构对龙门山断裂带应力积累的影响，我们在参考模型０的基础上另外引入了四个不同流变结构的岩石圈模型（模型１、模型２、模型３和模型４）．在这四个模型中，我们仅在参考模型０的基础上改变了岩石圈不同层位的黏滞系数，见表３和表４，而其他物质参数保持不变，且均加载边界条件１．我们复述模型０的参数取值以方便与其他四个模型做比较．各个模型的具体黏滞系数随深度分布见图８．模型０：上地壳表现为脆性；中、下地壳表现为韧性，且横向上青藏高原与四川盆地的黏滞系数有差异；岩石圈上地幔表现为韧性．黏滞系数的取值在上地壳为１．０×１０２３　Ｐａ·ｓ；高原下方中地壳为３．０×１０２０～５．０×１０２０　Ｐａ·ｓ，下地壳为３．０×１０１９～６．０×１０１９Ｐａ·ｓ，四川盆地中地壳为１．０×１０２２～１．２×１０２２Ｐａ·ｓ，下地壳为７．０×１０２１～８．０×１０２１　Ｐａ·ｓ；岩石圈上地幔为４．０×１０２０～６．０×１０２０Ｐａ·ｓ．模型１：整个岩石圈表现为脆性，岩石圈各层位有统一的较高的黏滞系数，为１．０×１０２３Ｐａ·ｓ，且横向上从青藏高原到四川盆地黏性无差异，具体黏滞系数随深度分布见图８ｅ．模型２：整个地壳表现为脆性，地壳各层位有统一的较高的黏滞系数，为１．０×１０２３Ｐａ·ｓ；岩石圈上地幔表现为韧性且黏滞系数与模型０中相同．具体黏滞系数随深度分布见图８ｆ．该模型中所描述的岩石圈流变结构符合Ｔａｐｐｏｎｎｉｅｒ等（２００１）支持的青藏高原全脆性地壳模型特征．模型３：上地壳表现为脆性；中、下地壳和岩石圈上地幔表现在为韧性，且横向上青藏高原与四川盆地的黏滞系数无差异（这是该模型与模型０的唯一区别），中、下地壳黏滞系数分别为３．０×１０２０～５．０×１０２０Ｐａ·ｓ和３．０×１０１９～６．０×１０１９　Ｐａ·ｓ；其他层位黏滞系数均与模型０相同．具体黏滞系数随深度分布见图８ｇ．模型４：仅将模型０中青藏高原中、下地壳的黏滞系数均减小为原来的２０倍，则其中、下地壳黏滞系数分别为１．５×１０１９～２．５×１０１９　Ｐａ·ｓ和１．５×１０１８～３．０×１０１８　Ｐａ·ｓ．这一取值仍在Ｃｌａｒｋ等（２００５）、石耀霖等（２００８）和Ｇｏｄａｒｄ等（２００９）对青藏高原中、下壳的黏滞系数的估计范围之内，为合理取值．引入模型３、模型０和模型４其目的是为了研究青藏高原与四川盆地中、下地壳黏滞系数差异从无到有、到增大差异这三种情况下，相应岩石圈模型应力积累分布的区别．由于结果显示模型４与模型３的应力积累结果分布十分相似，因此我们只给出边界条件１作用下模型１、模型２、模型３和模型０中经过芦山地震震中且垂直于龙门山断裂带的纵剖面上的应力积累速率分布状况，见图８．结果具体如下：模型１中，如图８ａ所示，整个地壳内部横向上应力积累速率均匀分布，没有应力集中的现象．纵向上应力积累速率随深度递增大，不同层位的应力增长率随层位加深、杨氏模量的增大而增加．芦山地震的震源处应力速率为－２．１ｋＰａ／ａ．模型２中，如图８ｂ所示，在龙门山断裂带下地壳的底部出现应力集中现象，且模型１中Ｍｏｈｏ面拐点处Ａ应力积累速率从－２．８ｋＰａ／ａ增大到模型２中的－４．３ｋＰａ／ａ．可见，在Ｍｏｈｏ面起伏的岩石圈结构中，构造应力积累在很大程度上受控于岩石圈的流变结构，这一点在我们的另外一个研究中也同样得到验证（柳畅等，２０１２ａ）．地壳其他部分应力积累分布与模型１中相似．韧性的岩石圈上地幔部分应力积累速率为０．模型３中，如图８ｃ所示，上地壳内龙门山断裂带下方并无应力集中现象出现，并且在上地壳内青藏高原的应力积累速率小于四川盆地的．韧性的中、下地壳和岩石圈上地幔部分应力积累速率为０．以上结果表明，模型１、模型２与模型３的应力积累结果均不能与龙门山地区的地震空间分布特征相对应（包括纵向与横向）．从应力积累的角度而言，这三种模型既无法解释汶川地震与芦山地震的发生，也无法解释如上文所述的龙门山断裂带地震空间分布特征（李勇等，２００６；滕吉文等，２００８；雷兴林等，２０１３）．只有模型０中的应力积累分布特征与龙门山地区的地震空间分布特征有较好的吻合，这一点在上文中已得到详细分析，这里不再赘述．模型０与其他几个模型相比较，其最重要的特征就在于模型中青藏高原与四川盆地的中、下壳的黏性差异．模型４中各层位的应力积累分布与模型０的结果相似．但是由于模型４相对于模型０增大了青藏２１４。