四川会理川-31 井水位变化与降雨量关系及异常识别
2021-03-15芮雪莲官致君林圣杰
芮雪莲,杨 耀,官致君,林圣杰
(四川省地震局,成都 610041)
0 引言
地下水位作为流体学科重要的物理观测量,是一个包含大气降水、气压、固体潮、地应力场等多种影响因素的复合参数[1]。其中,降水是地下水位变化的主要影响因素,它使得多数井孔所观测记录的水位呈现雨季上升、旱季下降且比较容易识别的年变规律,但当年降水量随着旱、涝年份的不均匀而出现明显差别时,就会成为地下水位异常识别的重要干扰因素[2]。车用太等人通过研究北京地区、福建汤坑与黑龙江尚志等典型观测井认为,降雨对深井水位的微动态影响并不是垂直入渗而是暴雨积水的荷载作用到含水层中去的结果[3-4]。为分析排除降雨对井水位的影响,学者们提出响应函数法[5-6]、卷积滤波法[7]、组合水箱模型法[8],这些方法在云南地区地下流体异常分析中已经得到了很好的应用[9-10]。此外,也有学者通过统计分析,研究了不同时期井水位与降雨量之间的关系及异常特征[2,11-13]。但目前针对四川境内的观测井相关方面的研究相对较少。
四川会理川-31 井位于南北地震带中南段,是国内地震监视跟踪工作重点区域。自开展观测以来,水位动态变化受降雨干扰明显,总体表现为雨季上升,旱季下降的特征。2019 年3 月20 日该井水位出现明显的上升,幅度达到0.728 m,同时该年雨季水位上升幅度也明显增大。为了判定两次水位上升异常是干扰造成还是由区域构造活动引起的地震前兆异常信息,本文在前人研究的基础上拟通过分析降雨量与水位变化幅度之间的相关性、井孔水化学组成及氢氧同位素特征,对两次水位上升的原因进行分析。
1 观测井概况及井水位异常特征
会理县处于南北地震带的中南段,地质构造复杂,地壳运动频繁,断裂发育。区内发育有安宁河深断裂带、磨盘山-昔格达断裂、宁会大断裂等。川-31 井位于会理县城内(102.25°E,26.65°N),处于安宁河断裂带南段与宁会断裂交界区域,距宁会断裂约3 km。宁会断裂沿宁南理力马河一线延展,走向NE,西南段被昔格达断裂带错断,东北端交汇于则木河大断裂,沿断裂带地震活动频繁,2008 年8 月30 日攀枝花-会理6.1 级地震距观测井仅60 km(图1)。
图1 川-31 井周边地质构造及2000 年以来地震分布图
观测井周边及附近区域地层为河湖相沉积,岩性主要是上白垩统雷打树组砂岩、下白垩统小坝组的砂岩、泥岩和上三叠至下侏罗统白果湾群砂岩、页岩,含煤。研究区内地下水资源丰富,纵横交错发育的断裂,破坏了含水层的连续性,同时也沟通了一些含水层间的水力联系,部分断层出现了充水及地下水的局部富集。区域内地下水的补给、径流、排泄条件及动态变化受气候、地貌及构造等因素的控制,除小型山间盆地及裂隙孔隙水外,一般以降雨补给为主。山地与河谷区的气象特征直接影响区域地下水的动态变化。
川-31 井深约260 m,钻孔揭露地层岩性为下白垩统小坝组砂岩、泥岩、泥质砂岩及砂泥岩互层,观测含水层深度为182~192 m。地下水类型为裂隙承压水,观测井位于山区,其地下水的动态变化直接受到区域气象特征的影响。
会理川-31 井水位于2019 年3 月20 日开始快速上升,并于4 月12 日开始出现回落,最大上升幅度达0.728 m。自观测以来,该井年变上升时间与当地雨季保持一致(每年5—10 月),但此次水位上升时间较往年相比,提前约2 个月。同时,2019 年6 月23—26 日水位也出现了明显上升。通常情况下,水位出现上升主要由以下原因引起[10,14]:①环境干扰:区域降水或者人为抽注水引起含水层内水量的变化;②区域构造应力作用下,含水层固体骨架发生变形从而导致孔隙压发生变化,进而引起水位发生变化。因此,为了进一步明确观测井水位发生变化的原因,收集相关降雨等辅助资料和进行井孔水化学特征分析是十分有必要的。
2 降雨量与水位变化关系分析
川-31 井所处地区属中亚热带半湿润气候区,据2009—2019 年数据统计,会理川-31 井所属地区年降雨量为791~1 503 mm,各年间的降雨量差异较大。一年之内,雨旱季分明,降雨量≥50 mm 的月份主要集中在5—10 月。
观测资料表明,川-31 井水位的变化受到该地区降雨的影响明显,年变动态为降雨补给型,水位变化主要表现为雨季上升,旱季下降的年变形态。由于川-31 井位于山间,考虑到降雨有效入渗及间接补给等其他因素,本文主要研究雨季降雨对水位的影响,而非年降雨量的影响。
2.1 水位变化与降水量相关性分析
2.1.1 雨季降水量与水位变化的关系
根据2009—2018 年每年水位变化幅度与降雨量数据分析,可认为水位变化与降雨量之间存在着密切的关系(图2)。且降雨对井水位变化的影响并不是短时间的,而是对此后年变化幅度均会产生影响。为进一步研究雨季降水量与年变幅度之间的关系,将各年雨季累计降雨量与年变幅度值做线性回归拟合二者之间的线性关系(图3),可以看出二者之间具有一定的相关性(相关系数R=0.60)。
图2 川-31 井井水位(月均值)与降雨量的同轴曲线
图3 井水位年变幅和雨季降雨总量的关系
2011 年会理地区降雨量为统计时段内(2008—2019 年)最低,全年降雨量仅793.8 mm,其中雨季(5—10 月)降雨量741.2 mm,占全年降雨量的93.4%,为2009—2018 年内唯一年降水量在900 mm 以下的年份。由此造成该变化周期内(2011 年5 月至2012 年5 月)低 水 位 值(4.945 m)较 起 始 水 位 值(4.527 m)低0.418 m,该变化周期内的低水位值为10 年(2009—2018 年)统计范围内的最低值。2012 年雨季降水升至1 023.7 mm,含水层水量得到补充,高水位恢复至往年水平,由此造成2012 年水位变化较其他时间明显增高。
2.1.2 月降水量与水位变化的关系
水位变化正常年份(2009—2018 年),会理川-31 井水位动态变化曲线形态为谷-峰-谷型,即井水位在春冬季节,该地区几乎无降水,井水位基本保持在较低水位,水位变化线较平稳,而在夏季由于该地区降水较多,水位受此影响出现明显上升,进入秋季之后,降雨明显减少,水位趋势下降逐渐保持平稳。本文对该地区雨季降雨量及其对应的月水位值变化进行一元线性拟合回归分析。水位月变幅值在拟合线1 倍方差线范围时,则认为该井水位月变化量与降水量符合该线性拟合关系(相关系数R=0.63)(图4)。
图4 井水位月变幅和月累计降雨量的关系
2012 年6 月降雨量为260.1 mm,但是该月水位变化幅度明显高于其他同等降雨量下的水位变化幅度值,这是由于2011 年降雨明显减少,年变水位低值到达历史低值,同时6 月降雨量较五月出现明显增加,导致本月水位出现明显上升造成水位变化幅值较大。此外通过数据分析表明,井水位变化不仅与累计降雨量有关,与降雨方式也存在密切的关系。2018 年10 月降雨量仅为62.8 mm,但是该月水位变化量幅度为81 cm,明显高于同等月降雨强度下的水位变化幅度,不在拟合范围内,这是由于9 月27 至10 月1 日连续降雨,且日均降雨量可达30 mm,导致水位升高,并在10 月2 日达到本月水位最高值,此后该月只有零星降雨,且日降雨量均小于5 mm,降雨对水位并没有形成影响。此期间内水位保持持续下降趋势。2017 年6 月、8 月和2019 年9 月水位月变幅值受降雨方式的影响超出拟合范围。
2.1.3 集中降雨量与水位变化的关系
川-31 井含水层岩性主要为紫红色砂岩,泥砂岩,透水性较好,同时由于其处于山地地区,地下水补给和排泄速度快。故该井对降雨响应快速,较强降雨后很快可见水位的上升变化,雨停后水位则快速回落,达成一个新的平衡。
由于降雨对水位的影响不仅与降雨量相关,同时还受到降雨形式、植被、地形、包气带岩性、地下水的埋深等众多因素的影响[15]。因此,本文在进行降雨量-水位月变化相关性拟合分析时,个别点偏离拟合线较远。为进一步研究降雨形式对水位的影响,明确二者之间的关系,本文对该井2009—2018年集中降雨时段水位的变化进行详细分析。
根据降雨及水位对降雨的响应特征,本文“集中降雨”指日降雨量在25 mm 以上或者连续多天降雨导致水位有明显上升变化时的降雨量。通过统计2009—2018 年每一次集中降雨时段的降雨量及对应的水位变幅,对其进行拟合发现,其拟合关系较好(相关系数R=0.84)(图5)。通过拟合水位变幅(y)受集中降雨量(x)影响的公式为
式(1)可以作为经验公式用于川-31 井降雨量与水位变化的关系分析,在今后的震情监视跟踪时用于分析水位变化的异常判定。但是由于降雨对水位的影响不仅与降雨量相关,同时还受到降雨强度、地质条件等其他因素的影响,且本次研究样本量有限,故此经验公式具有一定的适用范围。当集中降雨量明显超出这一范围时,需要进行重新分析。
图5 集中降雨与水位变幅的关系
2.2 利用降水量反演水位动态变化
王旭升等基于响应函数法描述降雨入渗补给的滞后特征并建立了观测井降雨-水位动态的组合水箱模型,该模型利用降水量反演水位变化特征,并通过模拟值与实测值对比去除降水量对水位的影响[8]。孙小龙等基于该模型,模拟了云南姚安井水位的变化曲线,并对该井的水位异常变化与降雨之间的定量关系进行了研究,得到了较为理想的应用效果[10]。本文基于前人的研究成果,利用孙小龙等人编写的Matlab 程序,更加客观性地研究分析了川-31 井水位变化与降水量之间的相关性特征。
首先,利用川-31 井水位2017—2018 年正常动态变化时期的降水与水位资料(取值均为月均值),拟合求取了模型计算所需的响应函数参数t0、β,动力学时间因子Th、Tjq、Tiq(Th反映的是断裂带贮存水的能力;Tjq反映的是降水直接渗入断裂带的能力;Tiq反映的是断裂带获得侧向补给的能力)及排泄基准面高度z0(可用断裂带周围的最低排泄点高程代替)[8],得出以上6 个参数的取值分别为z0=4.2,β=3.1,t0=0.1,Th=3.3,Tjq=1.3,Tiq=51.3。从井水位观测值与模拟值曲线对比图(图6)可以看出,在水位正常变化时期(2017—2018 年)模拟值与实测值吻合度较好,说明模拟效果理想,所得模拟参数可靠。模型的平均绝对误差为0.093 m。
图6 川-31 井水位观测值与模拟值对比及相对水位动态变化曲线
基于2017—2018 年模拟所得参数,对2019 年1—12 月水位进行模拟。为形象直观地说明水位实测值与模拟值之间的差值,本文引入相对水位的概念,即相对水位等于水位实测值与模拟值的差值。以相对水位±0.1 m 作为误差参考线(图6 中虚线)来判断水位是否异常。图6 中将相对水位超过误差线的事件分别标识为A、B、C、D、E 共5 个事件,A~D 事件均发生在雨季,造成该类异常的原因可能是由于本文所用模型在处理高频率降雨输入尚不完善造成,从而形成明显的雨季扰动。综合2.1.2 节中月降雨量与水位变化规律分析结果,可将A~D 事件视为由雨季扰动引起,不作为水位异常事件处理。
图6 中相对水位在2019 年2—9 月表现为负异常,表明实测水位值高于模拟水位值,即以E 事件为异常起始事件,后续此项异常持续至9 月。值得注意的是,该异常事件开始时间2 月,3—4 月会理地区月累计降雨量均小于10 mm,可认为无降雨扰动,其后异常动态无法用模型受到雨季扰动解释,故综合分析认为此类事件为区域构造活动引起的异常动态。5—9 月为当地多雨时段,降雨量明显增加,该段时间内异常产生的原因可认为是雨季扰动和区域构造活动叠加引起。
3 井水化学特征分析
地下水的化学组分和稳定同位素组成可以用于判定其水化类型、水-岩平衡反应特征及其组分来源[16-17]。利用Piper 图可以判定地下水的类型和地下水的混合作用;Giggenbach 图可判定地下水的水岩平衡状态,进一步揭示含水层的深浅部补给关系和地下水系统的开放与封闭程度[18],利用水样氢氧稳定同位素数据与大气降水线之间的关系研究地下水的运动过程[19]。本文对川-31 井水化学测试结果进行分析,以进一步了解该井地下水循环特征及补给来源。水样测试机构为中国地震局预测研究所,水化学组分测试结果见表1。
表1 川-31 井水化学成分及同位素分析结果
图7 川-31 井地下水水样测试结果分析图
Giggenbach 提出的Na-K-Mg 三角图可用来评价水-岩化学平衡状态,判断地下水循环深度等。从图7c 可知,测试结果落在Mg 端元附近,为未成熟水,表明该井地下水循环周期相对较快,埋藏较浅,水-岩之间尚未达到离子平衡状态,水岩作用仍在进行[18]。氢氧稳定同位素技术能够有效地识别地下水的来源与补给过程,可用于判定大气降水和地表水与地下水的补给关系[19]。川-31 井水样氢氧同位素测试结果位于中国西南地区大气降水线上(δD=7.54δ18O+4.84[21]),故认为该观测井地下水补给来源为大气降水(图7b)。
川-31 井水化学组分及稳定氢氧同位素分析表明,该井水循环深度较浅,井水来源主要为大气降水,与地表水水力联系密切,水位动态易受到降雨影响,故该井水位、水温出现异常时需分析是否存在降雨干扰。
4 结论
通过统计和数学物理模拟分析了降雨量与水位变化的相关性特征,并结合观测井地下水化学特征明确了川-31 井地下水循环特征及补给来源。
1)会理川-31 井水位变化与该地区降雨量之间关系密切。水位年变幅值与雨季累计降雨量之间存在一定的线性关系,2009—2019 年之间降雨量与水位变化幅值均拟合较好,符合线性相关的关系;雨季月降水量与水位月变幅之间呈线性相关,水位月变化不仅与降雨量有关,与降雨强度也密切相关。
2)川-31 井对降雨响应快速,集中降雨对井水位影响尤为明显,在一定范围内二者变化有较好的线性相关关系。通过对比模拟水位与实测水位,认为2019 年3 月水位明显上升为区域构造活动引起的异常事件,同年雨季水位与模拟水位值差别较大可能是降雨与前期构造活动叠加效应所致。
3)井孔水化学组分和氢氧同位素分析结果表明,会理川-31 井地下水以浅层溶解作用为主,水岩作用程度较弱;大气降水为地下水的主要补给来源。
4)综合上述分析,会理川-31 井2019 年3 月水位出现明显上升异常并非是降雨导致,而是与区域构造活动有关。同年雨季上升幅度较往年大,则可能是降雨与构造活动叠加作用造成。
本文在分析川-31 井降雨与水位变幅之间关系的过程中由于研究方法、数据等方面的限制,其相关系数不太理想,在今后的震情跟踪工作中需要丰富数据量以进一步研究完善。
致谢感谢审稿专家给本文提供宝贵的修改意见,孙小龙研究员为本文提供计算程序,凉山州防震减灾局王娟及会理县防震减灾局杨朝坤、罗文俊等人为本文提供了降雨量等其他相关数据,在此一并表示衷心的感谢!