地震学报  2007, Vol. 29 Issue (3): 265-273
唐山井水温的同震变化及其物理解释
石耀霖1, 曹建玲1, 马丽2, 尹宝军2,3,4    
1.中国北京100049中国科学院研究生院;
2.中国北京100036中国地震局地震预测研究所;
摘要:唐山市一口水井观测到多次同震水位振荡, 以及伴之的高精度测量到的深部层水温变化. 水温变化幅度与水位振荡幅度有关, 水位幅度变化数厘米到一米, 水温变化幅度为0.001℃~0.01℃, 而且总是温度下降, 震后一到数小时内水温恢复正常. 我们进行了有限单元法模型计算, 认为井水垂直振荡时搅动井水引起的弥散效应, 是造成同震水温变化的主要原因, 后续的热传导作用可以解释水温的复原过程. 今后如果在不同深度同时进行高精度井水温度观测可以检验该模型.
关键词地震    井水位    井水温    同震效应    
引言

文献中经常见到关于井水与地震关系的研讨(King et al,1999Yechieli,Bein,2002陆明勇等,2005a杨竹转等,2005). 关于井水位、 水温、 水化学成分等的变化是否能反映地震前兆问题,仍是人们感兴趣和继续探讨的课题(Igarashi et al,1992; Chadha et al, 2003; 曹新来,边庆凯,2004陆明勇等,2005b); 而某些井水水位可以像地震仪一样记录到地震波动、 像重力仪一样记录到固体潮汐,已经是公认的事实(Cooper et al, 1965; 王道,2002Wang et al, 2003; 庄光国,1994). 虽然对水井温度观测和分析已有报道和研究(车用太等,2003; 刘伟等,2000; 毛德培等,2002),特别是苏门答腊大地震后对我国大量观测井资料的整理和机制讨论(中国地震局监测预报司,2005),然而对单井水温的高精度测量和同震变化机制进行定量模型分析,文献中讨论仍然少见. 本文将以唐山一口水井观测为例,进行描述和模型分析.

本文讨论的井水现象简单,即每次地震波到达时,井水位发生振荡,而井深部水温降低. 不论地震的方位、 远近、 大小,记录的水温变化无一例外,总是发生降低,且较大的水位振荡导致的水温降低幅度较大,恢复原正常值所用的时间也较长. 谷元珠等(2003)曾经注意到潮汐、 地震及提探头等微扰动对水温的影响,并从水的动态扰动方面进行了定性的探讨和解释. 中国地震局监测预报司(2005)在编辑的《2004年印度尼西亚苏门答腊8.7级大地震及其对中国大陆地区的影响》一书中,也介绍了有关专家的一些定性讨论. 我们认为,我国为地震预报目的建立了许多地下水观测井台站,取得了世界不多见的丰富记录,我们应该对这些宝贵的数据进行认真的分析: 对于曲线的每一个变化,都要问一个为什么,都要试着给出定量的解释. 只有在对观测数据有了深刻理解,我们才能区分什么是正常变化、 什么是异常变化、 什么是干扰,才能进一步从异常中寻找是否有地震前兆异常,实现地震预报. 虽然用这样的做法取得成果似乎要缓慢一点,但通过科学的途径,一旦取得了正确的认识,获得的成果就会比较实在.

1 唐山水井的同震地温变化观测

唐山地下水观测井位于唐山市区大钊公园,一般称唐山井,又称唐山矿井或山西水2井. 该井构造上位于燕山褶断带与华北平原沉降带的结合部,北东向唐山断裂与北西西向隐伏断裂于此通过(张子广等, 19982002). 该井于1969年10月30日终孔,完钻井深286.6 m(图 1),第四纪冲积层厚度10.27 m,其下到150.84 m为石炭二迭纪砂岩,于150.84 m见含水层奥陶纪灰岩. 含水层厚度50 m左右,是该井孔的出水段,地下水埋藏类型为岩溶裂隙承压水,局部承压,封闭性不好. 由于逆断层造成地层重复,因此至240.63 m又见砂岩等石炭、 二迭纪地层. 钻孔完工后封孔至217 m,套管深度154 m,以下为裸孔. 该井位于城市地下水开采漏斗中心部位,水位埋深较大,且长趋势下降. 1983—2000年水位平均埋深51.7 m,现今约为62 m.

图 1 唐山井水文地质柱状图 (张子广等,2002)

自1981年1月开始采用SW40-1型水位计进行地下水位观测,可以记录到固体潮和有气压效应,水位也受到开采地下水的影响. 一般认为,由于唐山井正好打在断层上,而且该井的井-含水层系统的固有周期为21.9 s,与瑞利面波的振动周期相匹配,对地震波的响应比较敏感,能记录到全球范围内多数强震的地震波(张子广等, 19982002).

该井在125 m深处(约在水下63 m)安置有测量水温探头,测温分辨率可达0.000 1°C,从2001年9月开始观测. 自记录温度以来,该井对许多远震观测到同震水位变化的同时,也观测到了水温的同震变化,部分实例见表 1图 2. 一个显著的特点是,在所有记录到的地震中,水温均发生了降低. 降温过程一般比较迅速,降温持续时间一般与井水震波峰值出现时间有关,在初至波到达后开始降温,在地震波峰值到达时降温幅度达最大,而后慢慢恢复到原来的稳态温度值. 降温过程一般用10~20分钟,恢复一般用0.5~1小时,个别达数小时.

表 1 唐山井记录到的水位与水温同震变化的部分地震

图 2 唐山井记录到的部分地震(表1)的水位和水温的同震变化 (a)~(j) 分别为表1中的11次地震的记录, 地震编号在水位变化图中标出

为什么同震记录中水温总是降低?为什么一小时左右又可以恢复?下面本文就此问题进行探讨.

2 同震降温和震后恢复现象的解释

为解释这一现象,我们对降温机制提出两个设想,并进行了定量计算,以考察该设想的合理性. 计算中取水面为坐标原点,井水距水面深为155 m,套管距水面深度为92 m,温度计位于水下深度63 m.

第一个设想: 由于气、 水、 岩层、 钢质套管等热导率不同,造成水井稳态温度分布时井水温度略高于井壁温度,地震波到达时井壁含水层的水进出于井内,这种较冷的水混入较暖的水的混合过程,降低了井水的温度.

在不存在热源时,稳态热平衡方程为(梁昆淼,1998)

其中,T为温度,K为热传导系数. 在柱轴对称坐标系下该方程表示为

取离井管足够远处(r=5 m)为绝热边界,水表面和井底200 m深部为恒温边界条件,取上下温差为0.775℃.求解井内温度分布,观察它与随深度温度线性增加的均匀介质中的地温场的差异.钢套管、岩石、水和空气的热导率分别取为50,1.5,0.6和0.025 W/(m·K).求解结果表明,仅在井底和套管底部热导率差异对水温偏离均匀场较大,其余部分与均匀场差别极微小(图 3).模型虽然表明在套管底部和井孔底部,的确由于套管热导率高、岩石热导率居中、井水热导率低,这种热导率的差异造成的所谓"热折射",造成了套管底部和井底局部井中心水温高于相邻井壁温度,但是这种温度空间分布的差异有限,仅在这两个局部.而在远离这两个局部1 m以外的绝大部分深度内,水温与相邻井壁 处于相同温度.另外,这种温度的差异十分有限,尚不到0.001℃.这样局部微量的、仅冷 0.001℃的水混合于井水中,不可能造成井水千分之几到百分之一度的温度变化.

图 3 井孔内由于套管、 岩石、 水的热导率不同而造成的温度扰动 (a) 井孔底部附近等温线; (b) 套管底部附近等温线

第二个设想: 由于井水处于稳态时,水温一般随深度线性增加,上部的水温较低,下部的水温较高,在井水振荡过程中水受到扰动而产生温度变化. 控制井水温度的主要机制不再是热传导,而是井上部的较冷的水与井下部较暖的水的混合的水动力弥散过程. 这时,我们用涡旋热传递方程(热弥散方程)来研究其温度变化过程和特征( Pinsonet al, 2007):

湍流热扩散系数D为分子热扩散率D0与涡旋热扩散系数D*之和,后者可以比前者大数个数量级,与水流速度有关,一般D*=αivm,参量m一般在1~2之间(倪龙,马最良,2005). 由于该孔固有振荡频率在20 s左右,幅度可达几十厘米,因此孔内水体垂直运动速度幅值不妨估计为50 mm/s,时间平均也在30 mm/s以上. 在这样的水流速率下,湍流热扩散系数远远大于分子热扩散系数. 在水井深超过百米的条件下,其雷诺数达5 000以上,因此井孔内水运动可能已出现湍流,这时水流速对热弥散系数的影响将更为明显,弥散系数估计可以达到1 m2/s数量级. 在计算中将取不同数值试验计算,而其值最终确定将有待于在井孔内不同深度设置温度计,根据未来大地震时对水位和不同深度温度测量去计算.

在同震水位振荡时,上边界条件又可以考虑两种极端可能性:一种是假定水面来不及进行热交换,即取为绝热边界条件,这时下面温度较高的水加热表面的冷水,水面的温度达到升高的上限;另一种是假定水面温度为常数,即表层水如果趋于温度升高,其热量被迅速传入空气带走,水面温度始终维持常数,这是水面温度的下限.真实情况不会超过这两种极端情况的温度范围.下边界取为绝热边界条件;初始条件则均取为从水表面到井底线性增加,从0℃~0.775℃.

对于表层为绝热边界条件的第一种极端情况,又分为两种情况讨论: 第一种情况,如果热弥散系数在不同深度为常数,则解为图 4所示. 浅于二分之一水深的点,由于混入了下方较热的水,温度总是上升; 深于二分之一水深的点,由于混入了上方较冷的水,温度总是下降. 而温度测点位于距水面63 m深度,处于井水总深度的二分之一以上,似乎温度应该总是上升,而不是下降. 然而,考虑热弥散系数是速度(因而也是深度)的函数,则结果就会不同.

图 4 当上边界为绝热边界条件,D= 1 m2/s时,地震波扰动后0,5, 10,15分钟后的温度变化曲线

第二种情况,考虑热弥散系数由于水动力条件不同,在不同深度可以有不同的值: 在套管内部的水,由于质量守恒,震荡时速度相同; 而地震波到达时水被从没有套管的含水层位挤出或吸入时,鉴于质量守恒,在套管以下的水的垂直速度将随深度线性变化,即

而如果D随速度u线性变化,则按此模式,在接近孔底的一段距离内,水速很小,D就更小,因此热弥散系数很小,仅在距孔底60 m以上的范围内,活跃的弥散才会发生. 这样,计算出的温度同震变化如图 5a所示. 这时,地震波到达后5分钟,靠近水上表面的水体温度上升,而更深部的水温度下降,但分界线不再是在井深的二分之一处,而是位于较浅的40余米处; 10分钟后,靠近水上表面的温度上升段扩展到50 m深度; 15分钟后60 m以上均为升温,60 m以下温度仍然下降. 也就是说,最表浅的水温度一直上升,而稍深的水温度则先下降再上升,深度大于井深1/2的水体在振荡持续过程中温度一直下降. 其中观测点温度随时间的变化如图 5c所示,观测点温度下降幅度最大约为0.003℃. 如果热弥散系数在所有深度都变小1/2,计算结果见图 5bc所示. 其温度变化总的特征不变,但幅度和时间变化率会有些差别.

图 5 上表面为绝热边界条件的地震波扰动后的温度变化情况 (a)地震波扰动后5,10,15分钟后温度随深度变化曲线.热弥散系数随深度变化为:0~92 m深度时D=1 m2/s, 92~155 m深度时D线性衰减至零;(b)地震波扰动后不同时间温度随深度变化曲线.热弥散系数随深度变化为: 0~92 m深度时D=0.5 m2/s,92~155 m深度时D线性衰减至零;(c)测点(深度63 m)温度随时间的变化, 点线与图4定常热弥散率情况对应,虚线与图(a)的情况相对应,实线与图(b)的情况相对应

对于水的上表面为恒定温度的第二种极端情况,计算结果见图 6a所示. 这时水表层的温度不变或变化很小,所有深度的水温都下降. 深度越大,温度下降幅度越大,而且只要扰动持续,温度就随时间持续下降. 观测点温度随时间变化见图 6b所示,温度下降幅度可以超过0.01℃.

图 6 上表面为恒温边界条件的地震波扰动后的温度变化情况 (a)地震波扰动后5,10,15分钟后温度随深度变化曲线.热弥散系数随深度变化为:0~92 m深度时 D=1 m2/s,92~155 m深度时D线性衰减至零;(b)测点(深度63 m)温度随时间的变化

再对震后被扰动的温度恢复热平衡所需时间进行计算.假定井内水体温度受扰动而温度偏离平衡状态,用轴对称平面有限单元法计算恢复平衡所需要的时间.水和岩石的比热 分别取为4 186和860 J/(kg·K),密度分别取为1 000 kg/m3和2 536 kg/m3.井孔中心温度变化结果如图 7所示.可以看出,当温度扰动为0.001℃时,大约用1个多小时恢复;温度扰动为0.005℃时,则需数个小时恢复到原温度.该结果与观测的结果大体吻合,但略大于实际观测值.这与计算使用的是轴对称模型,而实际热传递是三维的有关.

图 7 初始温度为零度,井水温度存在10-3和 5×10-3℃扰动时的井水温度随时间的变化
3 讨论和结论

唐山井同震温度总是降低的现象,引起了不少观测者的关注. 鱼金子等(1997)谷元珠等(2003)都注意到了水动力学因素的重要性. 苏门答腊大地震后,对于“水位振荡而温度下降这样一个普遍存在的客观事实”(中国地震局监测预报司,2005),部分学者的“一种观点认为是地下水振荡过程中,水中的气体从溶液中分离、 溢出,……吸收周围介质的热量……,导致水温度降低”. “另一种解释是,井水含水层周边上层地下水随着振动效应的作用,加大了向下垂直运动的速率,低温水快速混入观测含水层,引起温度的快速下降”. 对这些设想的深入研究不能仅仅局限于定性的讨论,还需要定量化的研究. 笔者提出的震荡时水受扰动导致不同深度不同温度水的混合,从而造成同震水温变化的弥散模型,对此现象能够给予较好地定量化地解释.

人们熟知,热传递有多种形式. 对于水这一类的并非良导体的物质,热传导是一种缓慢的过程,热扰动显著影响到1 cm外需要十余分钟的时间,不能解释同震井水振荡时迅速的温度变化; 热对流是一种有效的热传递形式,但在井水振荡中,水虽然有向上向下的交替流动,但就整体来讲,长时间的平均速度和位移却是零,因此含平流项的热传递方程并不能解答此问题. 本文提出的主要机制是分子的弥散,如同清水中的一滴墨水会弥散开来,而且在静水中弥散慢,在流水或振荡的水中,特别是湍流发生、 存在涡旋时弥散快一样,分子动能不同的分子会弥散,而弥散系数与水的宏观速度密切相关,在静水中弥散系数很低,在同震水位振荡时弥散系数大大增加. 温度较高的一些高分子动能水分子弥散到冷的低分子动能的水中,以及温度较低处一些低分子动能水分子弥散到温度较高处,形成水温的变化,在一定的条件下,形成同震水温降低现象.

我们的模型计算讨论了所需的条件: 一种可能是如果井水表面由于与空气接触而便于迅速热交换,使井水水面维持常温,则降温的现象及其降温幅度不难得到解释. 这时井孔下部的较热的水因弥散作用渗入冷水而温度降低,接近上表面的较冷的水虽然因为弥散作用掺入了下面较热的水,但由于在与空气的界面上迅速将热传走,因此温度并不升高,整个井孔温度均有不同程度的下降. 另一种可能是如果井水表面不能迅速进行热交换,其最极端的情况是上表面处于绝热状态,这时井孔下部的较热的水将因渗入上部的冷水而变冷,井孔上部较冷的水将因渗入下部较热的水且不能从表面散出而变热. 如同上一节指出的,计算中还要考虑到振荡时套管内的水和更深部位与含水层相通的井孔中的水流速不同,因而弥散系数也不同,才能具体确定各个部位水温的变化特征. 实际情况应该介于这两种极端情况之间,虽然目前的观测资料尚不足以确定究竟哪种情况更接近实际,但无论哪种情况,都预期唐山井孔内水的振荡应该造成观测点温度同震下降. 该解释十分直接和简单,物理机制合理. 验证这一模型的最佳办法是在观测井孔内不同深度设置高精度温度观测探头,等待下一次大地震的同震变化记录,看看不同深度的部位观测到怎样的同震温度变化. 在更浅的部位观测温度变化将有助于确定哪种边界条件更合乎实际情况; 多个深度的观测将有助于我们对实际弥散系数的大小作出合理估计. 鉴于印尼苏门答腊大地震后我国许多台站都观测到同震地下水温变化,而且绝大多数均为降低,本文的模型解释可能具有普遍意义.

在地震观测,特别是在地震预报的研究中,我们不能对观测曲线知其然而不知其所以然,对每一个微小的“异常”变化,都应当力求明了其物理机制,获取其定量解释,以科学钻研精神和求实作风去探讨未知的奥秘. 在这种意义上,本模型的验证观测也是值得进行的,我们期待着能够看到这一实际检验观测结果.

参考文献
[1] 曹新来, 边庆凯. 2004. 华北地区强地震前地下水动态的重现性异常[J]. 地震学报, 26(增刊): 154-161.(1)
[2] 车用太, 刘喜兰, 姚宝树, 等. 2003. 首都圈地区井水温度的动态类型及其成因分析[J]. 地震地质, 25(3): 403-420. (1)
[3] 谷元珠, 车用太, 鱼金子, 等. 2003. 塔院井水温微动态研究[J]. 地震, 23(1): 102-108.(2)
[4] 梁昆淼. 1998. 数学物理方程[M]. 3版. 北京: 高等教育出版社: 529.(1)
[5] 刘伟, 胡久长, 顾申宜, 等. 2000. 海口ZK26井水温动态特征[J]. 华南地震, 20(1): 43-47. (1)
[6] 陆明勇, 黄辅琼, 刘善华, 等. 2005a. 地壳变形与地下水相互作用及其异常关系初探[J]. 地震, 25(1): 67-73.(1)
[7] 陆明勇, 牛安福, 陈兵, 等. 2005b. 地下水位短临异常演化特征及其与地震关系的研究[J]. 中国地震,21(2):269-279.(1)
[8] 毛德培, 王树明, 胡智文, 等. 2002. 大姚地热动态特征分析[J]. 地震研究, 25(1): 42-47.(1)
[9] 倪龙, 马最良. 2005. 热弥散对同井回灌地下水源热泵的影响[J]. 建筑热能通风空调, 24(4): 7-18.(1)
[10] 王道. 2002. 昆仑山口西8.1级地震的远程地下水动力学效应[J]. 内陆地震, 16(3): 213-223.(1)
[11] 杨竹转, 邓志辉, 赵云旭, 等. 2005. 云南思茅大寨井水位同震阶变的初步研究[J]. 地震学报, 27(5): 569-574.(1)
[12] 鱼金子, 车用太, 刘五洲. 1997. 井水温度微动态形成的水动力学机制研究[J]. 地震, 17(4): 389-396.(1)
[13] 张子广, 张素欣, 郑云贞, 等. 2002. 山西水2井和岳42井水位记震能力分析[J]. 西北地震学报, 24(3): 263-266.(3)
[14] 张子广, 万迪堃, 董守玉. 1998. 水震波与地震面波的对比研究及其应用[J]. 地震, 18(4): 399-404. (2)
[15] 中国地震局监测预报司. 2005. 2004年印度尼西亚苏门答腊8.7级大地震及其对中国大陆地区的影响[M]. 北京: 地震出版社: 418.(3)
[16] 庄光国. 1994. 台湾海峡7.3级地震水震波特征研究[J]. 地震, 19(2): 199-203.(1)
[17] Chadha R K, Pandey A P, Kuempel H J. 2003. Search for earthquake precursors in well water levels in a localized seismically active area of reservoir triggered earthquakes in India[J]. Geophys Res Lett, 30(7): 1 416.(1)
[18] Cooper J H H, Bredehoeft J D, Papadopulous I S, et al. 1965. The response of well-aquifer systems to seismic waves[J]. J Geophys Res, 70(16): 3 915-3 926.(1)
[19] Igarashi G, Wakita H, Sato T. 1992. Precursory and coseismic anomalies in well water levels observed for the February 2, 1992 Tokyo Bay earthquake[J]. Geophys Res Lett, 19(15): 1 583-1 586.(1)
[20] King C Y, Azuma S, Igarashi G, et al. 1999. Earthquake-related water-level changes at 16 closely clustered wells in Tono, central Japan[J]. J Geophys Res, 104(B6): 13 073-13 082.(1)
[21] Pinson F, Gregoire O, Quintard M, et al. 2007. Modeling of turbulent heat transfer and thermal dispersion for flows in flat plate heat exchangers[J]. International Journal of Heat and Mass Transfer(in Press).(1)
[22] Wang C Y, Dreger D S, Wang C H, et al. 2003. Field relations among coseismic ground motion, water level change and liquefaction for the 1999 Chi-Chi (MW=7.5) earthquake, Taiwan[J]. Geophys Res Lett, 30(17): 1 890.(1)
[23] Yechieli Y, Bein A. 2002. Response of groundwater systems in the Dead Sea Rift Valley to the Nuweiba earthquake: Changes in head, water chemistry, and near-surface effects[J]. J Geophys Res, 107(B12): 2 332.(1)